Systematics and Phylogeography of The Non-Ethiopian Speckled-Pelage Brush-Furred Rats (Lophuromys Flavopunctatus Group) Inferred From Integrative Genetics and Morphometry


 Background: The speckled-pelage brush-furred rats (Lophuromys flavopunctatus group) has been difficult to define given conflicting genetic, morphological, and distributional records that combine to obscure meaningful accounts of its taxonomic diversity. In this study, we inferred the systematics, phylogeography, and evolutionary history of the L. flavopunctatus group using maximum likelihood and Bayesian phylogenetic inference, divergence times, historical biogeographic reconstruction, and morphometric discriminant tests. We compiled comprehensive datasets of three loci (two mitochondrial [mtDNA] and one nuclear) and two morphometric datasets (linear and geometric) from across the known range of the genus Lophuromys.Results: The mtDNA phylogeny supported the division of the genus Lophuromys into three primary groups with nearly equidistant pairwise differentiation: one group corresponding to the subgenus Kivumys (Kivumys group) and two groups corresponding to the subgenus Lophuromys (L. sikapusi group and L. flavopunctatus group). The L. flavopunctatus group comprised the speckled-pelage brush-furred Lophuromys endemic to Ethiopia (Ethiopian L. flavopunctatus members [ETHFLAVO]) and the non-Ethiopian ones (non-Ethiopian L. flavopunctatus members [NONETHFLAVO]) in deeply nested relationships. There were distinctly geographically structured mtDNA clades among the NONETHFLAVO, which were incongruous with the nuclear tree where several clades were unresolved. The morphometric datasets did not systematically assign samples to meaningful taxonomic units or agreed with the mtDNA clades. The divergence dating and ancestral range reconstructions showed the NONETHFLAVO colonized the current ranges over two independent dispersal events out of Ethiopia in the early Pleistocene.Conclusion: The phylogenetic associations and divergence times of the L. flavopunctatus group conform to demonstrated hypotheses surrounding the paleoclimatic and ecosystem refugium impacts on the evolutionary radiation of rodents dependent on stably humid conditions in the East Africa region. The overlap in craniodental variation between distinct mtDNA clades among the NONETHFLAVO suggests unraveling underlying ecomorphological drivers is key to reconciling taxonomically informative morphological characters. The genus Lophuromys requires a taxonomic reassessment based on extensive genomic evidence to elucidate the patterns and impacts of genetic isolation at clade contact zones.

. As such, they are an essential ecological component in these biodiversity hotspots, serving both as prey to raptors and small carnivores and as predators of invertebrates [11,37]. Moreover, their association with ecosystems characterized by stable annual precipitation makes them models for investigating how changing ecosystem-climate processes impact phylogeographic and speciation trends. Altogether, they demand stable taxonomic accounts to guarantee accurate appraisal of their diversity, ecological roles, and ecosystem functions.
The current distribution of the genus Lophuromys re ects relatively well-structured phylogeographic patterns. The L. sikapusi group spans a pantropical African range, from western Guinea to western Kenya, the Kivumys group is restricted to the Congo Basin east of the Congo River and the Albertine Rift, and the L. avopunctatus group is distributed primarily in eastern to central Africa [5,7,10,11,22]. Despite this remarkable geographic range, the spatiotemporal in uence of geographical features and climatic oscillations on evolutionary radiation in the genus Lophuromys remains largely unknown [7]. Studies using larger genomic datasets, like Komarova et al. [16], have uncovered complex reticulate evolution and recurrent mitochondrial introgression among the Ethiopian avopunctatus members, clearly illustrating that single-gene phylogenies, especially mitochondrial loci, should not serve as the exclusive basis for taxonomic assignment. For the non-Ethiopian avopunctatus members, even knowledge of mitochondrial DNA (mtDNA) diversity is limited. There is a necessity rst to test the extent to which the mtDNA re ects taxonomic units and biogeographical trends across its distribution, and then to contrast it with nuclear data.
In this study, we evaluated the taxonomic limits and biogeographic patterns in the genus Lophuromys using a comprehensive mtDNA (Cytochrome b; CYTB) dataset. We then focused on the non-Ethiopian avopunctatus members and complemented the CYTB alignment with Cytochrome c oxidase I (COI), and Interphotoreceptor retinol-binding protein (IRBP), and two morphometric datasets (geometric landmarks and linear measurements). The speci c aims were i) to elucidate the systematics of the NONETHFLAVO in the context of their position in the genus Lophuromys and ii) to investigate when and how the NONETHFLAVO diverged and dispersed to colonize their present ranges. Methods

Sampling
We compiled three datasets for the combined genetic and morphometric analyses. Sampling across the genus Lophuromys was possible for CYTB only and covered the currently known range of the genus (Fig. 1). Sampling for the full dataset (CYTB, COI, IRBP, and linear and geometric data) was possible for the non-Ethiopian L. avopunctatus members only, for which we sampled the currently known range, representing type localities (or their environs) of all species currently classi ed under or associated with the group (Fig. 1, Additional le 3 Fig. S1, Additional le 1 Table S1). The skulls are deposited at the Field Museum of Natural History, Chicago, USA (FMNH), Kunming Institute of Zoology, Kunming, China (KIZ), and National Museums of Kenya, Nairobi, Kenya (NMK).

Genetic data
Total DNA was extracted from muscle or liver tissue preserved in absolute ethanol at -80°C using the sodium dodecyl sulfate method [83]. The DNA was PCRampli ed using gene-speci c primer pairs (Additional le 2 Table S2). The PCR reaction template comprised of 20 µl volumes (0.5 µl primer pairs, 10 µl PCR Master Mix, 8.5 µl water, and 0.5 µl DNA template); the cycling temperature, time settings, and primers were speci ed as shown in Additional le 2 Table S2.
The ampli ed product was sequenced in forward and reverse directions using the ABI Genetic Analyzer (Applied Biosystems), assembled in Geneious Prime® 2020.2.4 (https://www.geneious.com, Accessed September 2020), and aligned in Aliview v.1.26 [84] using MUSCLE [85]. After dropping duplicates and sequences with a high ratio of gaps/ambiguous bases, we retained 803 CYTB sequences, of which 316 were newly generated, and the rest downloaded from GenBank [86] and the African Mammalia database [87] (Additional le 1 Table S1). From the new CYTB sequences, we subsampled from the unique haplotypes and extracted 138 COI and 100 IRBP sequences, which were aligned separately and concatenated in SequenceMatrix [88]. The alignment of concatenated loci was available for the non-Ethiopian L. avopunctatus members only and comprised 91 sequences, 3088 bp long (1140 bp CYTB, 717 bp COI, and 1231 bp IRBP), after matching similar sample identi cations. We con rmed that there were no premature stop codons, indels, or heterozygous bases in MEGA X v.10.1.8 [89] and resolved heterozygous bases in the IRBP alignment using PHASE [90] in DnaSP v.6 [91]. All the unique new sequences were submitted to GenBank (accession numbers MW464441 -MW464606).

Morphometric data
Morphometric variation among the non-Ethiopian L. avopunctatus members was inferred using a linear dataset of 725 skulls [310 , 363 , and 23 unsexed specimens] and a geometric dataset of 635 two-dimensional cranial images [278 , 338 , and 19 unsexed specimens] (Additional le 1 Table S1). The samples were age-classi ed based on the stage and pattern of M 3 wear into three age classes: young adults -fully erupted M 3 but very little to no visible wear, adults -medium wear on M 3 , and old adults -medium to extensive M 3 wear. Consequently, the geometric dataset comprised 29% young adults, 40% adults, and 31% old adults, while the linear dataset comprised 28% young adults, 41% adults, and 31% old adults. The samples' assignment to separate sex and age categories was used to illustrate the potential impacts of ontogenic variation and sexual dimorphism, which can confound differences between taxonomic units [92,93]. We used TPSUtil v.1.74 and TPSDig2 v.2.30 [94] to digitize 37 landmarks on the 2-dimensional skull images (Additional le 3 Fig. S2) and processed the resulting dataset in MorphoJ v.1.07a using Generalized Procrustes Analysis (GPA). The GPA untangles shape and size to produce centroid size (CS) and Procrustes coordinates. We performed a regression analysis with CS as an independent variable and the Procrustes coordinates as a dependent variable and used the resulting regression residuals as shape variables, free of CS variation, for consequent analyses. For the linear craniodental variation analysis, we used the same measurements and extraction techniques as in Onditi et al. [72].  The mitochondrial phylogeny of the genus Lophuromys was reconstructed from an alignment of 241 CYTB sequences (1140 base pairs long, 711 distinct  patterns, 443 parsimony-informative, 88 singleton sites, and 609 invariant sites), which included single longest sequences of each haplotype identi ed in the initial 803 sequences. The alignment represented all the species currently recognized in the genus Lophuromys [5,22], except L. medicaudatus and L. eisentrauti, and L. dieterleni for which we could not obtain representative new material or publicly available sequences. Sequences of Acomys ignitus, Deomys ferrugineus, and Uranomys ruddi, downloaded from GenBank, were used as outgroups. We used maximum likelihood (ML) and Bayesian inference (BI) methods for phylogenetic reconstructions, based on a GTR+F+G4 model of nucleotide substitution, which was identi ed as the best-tting under the Bayesian information criterion (BIC) in ModelFinder [95]. The ML analysis was performed using IQ-TREE v.1.6.12 [96] in PhyloSuite v.1.2.2 [97] using 100,000 ultrafast bootstrap replicates [98] to estimate branch support (BS). The BI analysis was performed in MrBayes v.3.2.7a [99] with two independent runs involving 10 million generations each, sampled every 1000 th run, using the reversible-jump Markov chain Monte Carlo (MCMC) [100] to estimate posterior probability support [PP]. The BI results were visualized in Tracer v.1.7.1 [101] to diagnose convergence using the effective sample size values (ESS), with values >200 considered adequate. The majority-rule consensus tree was annotated after discarding 25% as burn-in. The resulting trees from the ML and BI analyses were graphically edited in FigTree v.1.4.4 [102].

Species delimitation and genetic diversity analysis
Initial principal component analysis (PCA) tests of the morphometric datasets showed both the linear and geometric craniodental characters did not cluster samples consistent with current taxonomic units a priori. Therefore, we used the CYTB dataset to delimit operational taxonomic units (OTUs) representing valid species. We used species delimitation methods that can reliably identify common species units without prior assignment of samples to taxonomic units and implemented both tree-based and distance-based algorithms. For tree-based species delimitation, we used the branch-cutting method (BCUT, Mikula [103]), the multi-rate Poisson Tree Processes algorithm (mPTP, Kapli et al. [104]), and the single threshold general mixed Yule coalescent model (GYMC, Fujisawa and Barraclough [105]). We used the genus-wide ML tree as input in BCUT and mPTP analyses and a time-calibrated tree reconstructed in BEAST2 v.2.6.3 [106] for GMYC. The BCUT and GMYC analyses were performed in R v.4.0.3 [107] using functions provided by the author for the former and the splits package [108] for the latter. The mPTP analysis was implemented using the command-line options with four MCMC runs of 500 million generations, each sampled every 50,000 runs with a 10% burn-in, with convergence con rmed from a visual inspection of the combined likelihood plot. Finally, the distancebased delimitation was performed using the Automated Barcode Gap Discovery method (ABGD, Puillandre et al. [109]). The same genus-wide CYTB alignment for the phylogenetic reconstructions was used as input. The analysis was run in the ABGD web server (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html, Accessed 10 November 2020) using the K80 Kimura measure of distance, 0.001 -0.1 prior bounds for intraspeci c divergence, and a 0.75 relative gap width.
We used haplotype networks to inspect further the genealogical relationships between the delimited OTUs. Haplotype networks visualize genetic relationships among haplotypes, and because they do not force branching schemes, they may re ect evolutionary relationships better than the phylogenetic trees [110]. The haplotype networks were reconstructed using haplotypes generated in DnaSP and visualized in PopART v.1.7 [111] based on the Median Joining Network algorithm [112]. The genetic divergence within and between the delimited OTUs was explored using various indices of genetic diversity estimated in DnaSP, including the number of haplotypes (h), the probability that randomly selected haplotype pairs are different -haplotype diversity (Hd), the mean number of nucleotide differences (θ), and θ per site between randomly selected sequence pairs -nucleotide diversity (π). The genetic distances between and within the resolved OTUs/clades were estimated in MEGAX based on the number of nucleotide differences per site averaged between sequence pairs (uncorrected pdistances).

Estimation of divergence times
The divergence between main clades in the genus Lophuromys was inferred using the genus-wide CYTB alignment based on the coalescent-based approach in BEAST2. We applied secondary calibrations of the most recent common ancestor (MRCA) since Lophuromys has no fossil record. Two secondary calibration points were speci ed; the divergence between the L. sikapusi and L. avopunctatus groups, which was estimated by Aghova et al. [113] to ca. 3.71 million years ago [Mya] (con dence: 2.66-5.05) and the root node of the subfamily Deomyinae (having included sequences of Acomys ignitus, Deomys ferrugineus, Uranomys ruddi as outgroups). According to Aghova et al. [113], diversi cation within Deomyinae commenced ca. 13.8 Mya (95% highest posterior density interval [HPDI]: 12.04-16.01). We used lognormal priors with a mean of 1.31 and standard deviation (SD) of 0.1 (median 3.71 Mya) for the divergence between the L. sikapusi group and L. avopunctatus group and a mean of 2.628 and SD of 0.06 (median 13.8 Mya) for the MRCA of the Deomyinae subfamily members. Because the three genes, CYTB, COI, and IRBP, were available for the non-Ethiopian L. avopunctatus members only, we estimated a species tree of the group using the StarBEAST2 package [114] of BEAST2. The time-calibration was based on the divergence between the L. sikapusi and L. avopunctatus groups as speci ed above, following the inclusion of an L. sikapusi sequence as an outgroup. Two separate and unlinked substitution, clock and tree models corresponding to the mitochondrial (CYTB+COI) and nuclear (IRBP) loci were set, tted with uncorrelated lognormal clock and Yule speciation models. The time-calibrated phylogeny of the genus Lophuromys and the non-Ethiopian L. avopunctatus species tree were implemented with two MCMC runs, each 100 million generations-long, sampled every ten thousand runs were performed. The sampling convergence was assessed in Tracer; all the parameters had ESS values >400. The runs, including tree and log les, were combined in LogCombiner after discarding 10% as burn-in. The trees were summarized in TreeAnnotator and graphically edited in FigTree.

Biogeographical analysis
We reconstructed species ancestral ranges in RASP v.4.2 [115,116], based on the dispersal-extinction cladogenesis (DEC) model [117] which was selected with the BioGeoBEARS R package [118] best-tting to our dataset. The DEC model uses a species tree (with branch lengths scaled to evolutionary divergence times) and the geographical areas where the species (tree tips) occur to estimate ancestral ranges. The input tree was reconstructed from a reduced (single sequences from each GMYC-delimited OTU) time-calibrated genus-wide CYTB tree based on the same secondary calibrations as above. The major biogeographic ecoregions were de ned according to Dinerstein et al. [119] [https://ecoregions2017.appspot.com/ -Accessed 5 th November 2020], with slight modi cations. A total of six ecoregions were used; Albertine Rift montane forests, Guinea-Congo forests, East African montane forests, Eastern Arc forests, Ethiopian montane forests, and Southern Rift Montane forests. Because neither the ETHFLAVO nor NONETHFLAVO are monophyletic and range overlap exists between the Kivumys group, L. sikapusi group, and NONETHFLAVO, dispersal was allowed between all the ecoregions.

Morphometric analyses
The linear variables were initially transformed by natural logarithms to enhance their multivariate normality. The presence of outliers in the geometric dataset was checked for in MorphoJ for each species group, while in the linear dataset, we used Tukey's 1.5*IQR rule using a custom R script (http://goo.gl/UUyEzD, Accessed 1 st October 2020). From the combined 725 skulls for linear morphometry, <2% outliers existed for any of the 14 measurements across species groups, therefore, they were simply replaced with the respective group mean for each measurement. In the geometric dataset, a single sample identi ed as an outlier and was simply excluded from consequent analyses. The craniodental differences between groups were explored using discriminant function analysis (DA) in IBM SPSS Statistics v25 based on the within-group covariance matrices. In the DA, each group was assumed to have equal prior probabilities, so that cases were equally assignable to any group regardless of sample size. To test how classi cation accuracy compared to random assignment, we used the leave-one-out cross-validation model, where a discriminant function classi es cases based on all other cases except itself. Discriminant analysis is preferable when delimitating interspeci c morphological differences due to its ability to estimate the combination of characters that best distinguish groups [93,120].
Statistical signi cance of between-clade differences was estimated with the multilevel pairwise comparison using permutational MANOVA (multivariate analysis of variance) in the pairwiseAdonis R package [121]. We also used dendrograms of group mean clusters following MANOVA (performed using the manovacluster MATLAB function [www.mathworks.com/help/stats/manovacluster.html?s_tid=srchtitle, Accessed 1 st October 2020]) to visualize the multivariate craniodental relationships between clades. The three main species groups in the genus Lophuromys were well-supported (BS and PP >0.95) and occupied relatively speci c geographic areas (Fig. 1). The L. avopunctatus group was distributed primarily in highland regions of east and east-central Africa, with ETHFLAVO1 and ETHFLAVO2 being endemic to Ethiopia and the NONETHFLAVO1 and NONETHFLAVO2 spanning a broader range over the Eastern Afromontane Highlands south of Ethiopia (Fig. 1). The L. sikapusi group traversed the Guinea-Congo forest belt, with a primarily west to central Africa range (Fig. 1). The Kivumys group distribution, on the other hand, was restricted to the Albertine Rift and adjacent lowlands to the east of the Congo River (L. luteogaster), where it overlapped ranges with the L. sikapusi group and the NONETHFLAVO (Fig. 1).
Overall, the between-group genetic distances (uncorrected p-distance) was highest in Kivumys group versus L. sikapusi group (16.9%), the Kivumys group versus L. avopunctatus group was comparably distant at 16.4%, and the L. sikapusi group versus L. avopunctatus group were relatively less differentiated (11.5%).

Mitochondrial phylogeny of the L. avopunctatus group
Within the L. avopunctatus group, 12 main clades were resolved from the non-Ethiopian samples (NONETHFLAVO); three in the rst subgroup -NONETHFLAVO1 -and nine in the second subgroup -NONETHFLAVO2 (Fig. 2). Among the Ethiopian samples (ETHFLAVO), the number of clades in the two subgroups differed between ML and BI trees, with either one (L. chrysopus) or two (L. chrysopus and L. simensis 1) in the rst subgroup -ETHFLAVO1 -and the rest (9-10) in the second subgroup -ETHFLAVO2, (Fig. 2, Additional le 3 Fig. S3). The ETHFLAVO1 and ETHFLAVO2 subgroups were separated by an 8.74% genetic p-distance, slightly higher than the 5.8% p-distance that separated NONETHFLAVO1 and NONETHFLAVO2. Over-all, the p-distances between clades corresponding to the NONETHFLAVO (NONETHFLAVO1 and NONETHFLAVO2) were consistently higher than within clades, and were comparable, albeit averagely lower, to the p-distances between the ETHFLAVO clades -ETHFLAVO1 and ETHFLAVO2 ( Table 2).
The second NONETHFLAVO subgroup, NONETHFLAVO2, contained nine distinct clades -L. machangui, L. sabuni, L. makundii, L. dudui, L. rita, L. cf. cinereus, L. laticeps, L. stanleyi, and L. zena ( Fig. 2 and Fig. 3). The phylogenetic relationships and divergence times between the clades are shown in Fig. 2 and Fig. 3, their respective geographic ranges in Fig. 1 and Additional le 3 Fig. S1, and their genetic and evolutionary diversity in Table 1 and Table 2. Table 1 Estimates of evolutionary divergence between and within clades in the L. avopunctatus group. The number of base substitutions per site from averaging over all sequence pairs between groups are shown for within-clade (shaded diagonal) and between-clade (upper matrix: actual average site differences, lower matrix: uncorrected p-distances) comparisons estimated in MEGA X [1]. The non-Ethiopian L. avopunctatus members are highlighted in bold fonts. 2.3. Concatenated mitochondrial and nuclear phylogeny of the non-Ethiopian L. avopunctatus members The concatenated mitochondrial tree of the NONETHFLAVO was generally congruent to the genus-wide CYTB topology, with minor differences in sister relationships between clades (Additional le 3 Fig. S4). The NONETHFLAVO1 subgroup was distinct from NONETHFLAVO2, each separately monophyletic (Additional le 3 Fig. S4). In the NONETHFLAVO1 subgroup, notable differences with the CYTB topology included the paraphyly of L. zena + L. stanleyi clade, which contrasted the monophyly in the CYTB tree. The L. zena and L. stanleyi clades were also positioned at the root of NONETHFLAVO2, unlike in the CYTB tree. Except for L. laticeps and L. rita, which were not successfully sequenced for COI and therefore not included in the concatenated mitochondrial analysis, the rest of the NONETHFLAVO clades maintained congruent topologies as in the CYTB tree (Additional le 3 Fig. S4). On the other hand, the nuclear (IRBP) phylogeny did not correspond to the CYTB or concatenated mitochondrial tree, with most clades included in polytomies (Additional le 3 Fig. S5). In the NONETHFLAVO1 subgroup, L. aquilus merged with verhageni in monophyly while the L. kilonzoi samples remained monophyletic but not sister to the L. aquilus + L. verhageni (Additional le 3 Fig. S5).

Mitochondrial species delimitation, genetic distances, and networks
Each of the four delimitation methods produced an incongruous number and topology of splitting OTUs based on the genus-wide CYTB trees (Fig. 2). The mPTP identi ed 25 OTUs which differed from the 39 identi ed by BCUT, 46 identi ed by ABGD, and 56 identi ed by GMYC (Fig. 2). The L. sikapusi and L. chrysopus clades were consistently split into at least three OTUs across the methods, except in mPTP (Fig. 2). Several other clades were split as multiple OTUs by at least one of the delimitation methods, including those of NONETHFLAVO2 [L. zena, L. stanleyi, L. dudui, and L. machangui] (Fig. 2). Based on currently recognized species in literature and the haplotype networks, we resolved the 25-56 delimited OTUs to represent 33 clades. Of these, there were two clades in the Kivumys group, eight in the L. sikapusi group, and 23 in the L. avopunctatus group [12 clades corresponding to the NONETHFLAVO and 11 clades corresponding to the ETHFLAVO (Fig. 2)].
The evolutionary diversity (uncorrected p-distance) within the NONETHFLAVO clades (0.2% -1.6%) was systematically lower than between-clade diversity (2.39% -6.14%) ( Table 2). The haplotype networks of the L. avopunctatus group (combined ETHFLAVO and NONETHFLAVO) depicted composite genealogical relationships between clades that were not apparent in the phylogenetic trees, but altogether suggested a common evolutionary origin (Fig. 4).
The NONETHFLAVO clades with lower π and θ, such as L. rita, L. verhageni, and L. aquilus, had fewer haplotypes (Fig. 4, Table 2). Also, the more broadly sampled clades such as L. zena, L. stanleyi, and L. machangui had more haplotypes and higher haplotype diversity than those from single/fewer localities such as L. aquilus and L. verhageni (Table 3). These broadly sampled clades also revealed that haplotype networks were only slightly in uenced by sampling coverage, such that, within a clade, different localities were not uniquely systematically clustered (Additional le 3 Fig. S6).  (Fig. 3). Internal divergences within the L. avopunctatus group were more recent than in the Kivumys group and sikapusi group; with the oldest lineage, L. chrysopus, appearing ca.  5. The divergence times of the NONETHFLAVO based on the species tree of the three genes (CYTB + COI + IRBP) was generally more recent than the genus-wide CYTB estimates but maintained the relationships between clades (Additional le 3 Fig. S7). Because the time-calibrated phylogeny using the three genes was possible for the NONETHFLAVO only, we relied on the CYTB divergence dates for all divergence references.

Historical biogeography of the genus Lophuromys
Divergence within the genus Lophuromys likely originated in the Guinea-Congo/Albertine Rift forests, from where several dispersal events (28 dispersals versus six vicariance events) led to the colonization of current ranges (Fig. 5). These dispersals mostly occurred within ecoregions (mainly in the Guinea-Congo and Ethiopian Highlands forests) than between ecoregions (Fig. 5). The divergence in the Kivumys group likely originated in the same area as the genus, while in the Lophuromys subgenus, the Guinea-Congo forests formed the ancestral range, after which the L. sikapusi group remained in the Guinea-Congo forests while the L. avopunctatus group dispersed to the Ethiopian Highlands. From the Ethiopian Highlands, the NONETHFLAVO species colonized current ranges over two southward dispersal events (Fig. 5). The rst dispersal was by the NONETHFLAVO1 ancestor to the East African montane and Eastern Arc forests after which vicariance caused consequent divergences (Fig. 5). The NONETHFLAVO2 ancestor later dispersed to the Albertine Rift forests, from where both dispersal and vicariance events resulted in the colonization of the Congolian forests, East African montane forests, Eastern Arc forests, and the Southern Rift Montane forests (Fig. 5).

Morphometric analysis of the non-Ethiopian L. avopunctatus members
Morphometric analyses were conducted on combined datasets (age and sex) because sexual dimorphism and ontogenic variation did not impact discriminant tests between the CYTB clades. Overall, L. dudui had the smallest skull (based on CI and CS), while the L. aquilus skulls were the largest (Fig. 6, Additional le 2 Table S3). The combined-clades' morphospace following linear (DA LIN ) and geometric (DA GEO ) discriminant tests overlapped randomly with no evident systematic pattern delimiting the CYTB clades. Therefore, we partitioned the datasets into two groups around the two phylogenetic subgroups; NONETHFLAVO1 (L. aquilus, L. verhageni, L. kilonzoi) and NONETHFLAVO2 (L. sabuni, L. makundii, and L. machangui, L. stanleyi, L. dudui, L. laticeps, L. cf. cinereus, and L. zena).
The DA LIN and DA GEO classi cation results were consistent over-all, however, most clades were more correctly classi ed (more distinguishable) by DA GEO , especially in NONETHFLAVO2 (Fig. 6, Additional le 3 Fig. S8). Between-clade differences between the two skull datasets were not unidirectional, with DA GEO achieving higher correct classi cation than DA LIN in NONETHFLAVO1 but not in NONETHFLAVO2 (Table 3). In NONETHFLAVO1, all three clades were distinct, with DA correctly classifying >85% of each clade into the respective given group ( Table 3). The L. aquilus and L. verhageni were the most correctly classi ed by either DA LIN or DA GEO (Table 3, Fig. 6). The L. verhageni skulls were smaller than the adjacent L. aquilus or L. kilonzoi (Additional le 2 Table S2). The L. verhageni and L. aquilus skulls were more closely related to each other more than to L. kilonzoi (Fig. 6, Additional le 3 Fig. S8, Table 3).
In NONETHFLAVO2, the morphospace of L. zena and L. stanleyi markedly overlapped in DA LIN and DA GEO , between themselves and with several other clades, mainly L. cf. cinereus, L. dudui, and L. laticeps (Fig. 6, Table 3, Additional le 3 Fig. S8). The L. laticeps skulls were highly indistinguishable from other clades, being least correctly classi ed in the NONETHFLAVO2 subgroup and the combined pool of all clades.
The range-restricted clades, such as L. verhageni, L. aquilus, and L. makundii, were less ambiguously delimited and highly correctly classi ed by DA LIN and DA GEO (Fig. 6, Table 3, Additional le 3 Fig. S8). In contrast, more broadly sampled clades such as L. zena and L. stanleyi were less distinctly discriminated against from other clades (Fig. 6, Table 3). Statistical differences between clades based on multilevel pairwise comparison were signi cant, even for clades with overlapping morphospace (Additional le 2 Table S4).

Phylogenetic relationships within the genus Lophuromys
The deeper phylogenetic relations in the genus Lophuromys, including the validity of the subgeneric divisions (Lophuromys and Kivumys) and older lineages (Kivumys group and L. sikapusi group), and their respective monophyly have been relatively uncontested in recent checklists [5,6,22]. In contrast, species accounts in the 'speckled pelage' groups, the L. avopunctatus group, combining the Ethiopian endemics (herein as Ethiopian L. avopunctatus members -ETHFLAVO) and the non-Ethiopian ones (herein as non-Ethiopian L. avopunctatus members -NONETHFLAVO), have changed rapidly recently. In consensus, our genus-wide phylogenetic inference based on the CYTB gene supports the deep divergence of the genus Lophuromys into three distinct deeply-diverged groups that correspond with the widely recognized species groupings; i) Kivumys group (Kivumys subgenus), ii) L. sikapusi group (Lophuromys subgenus), and iii) L. avopunctatus group (Lophuromys subgenus). The two Lophuromys groups -L. sikapusi group and L. avopunctatus groupare separated by a much lower mtDNA divergence (p-distance) compared to the almost equidistant p-distance separating them from the Kivumys group. Based on the CYTB divergence, the Kivumys group appear distinct enough to be classi ed in a distinct genus, with the L. sikapusi group and L. avopunctatus group also distinct enough to be included in separate subgeneric groups. However, testing these conjectures hinges on broader nuclear loci sampling and morphological characterizations across the genus. Overall, these ndings broadly agree with the current classi cation of species into the Kivumys group [6]. However, we importantly highlight that classifying species into the L. sikapusi and L. avopunctatus groups based on morphological characterization is essentially ambiguous unless backed up by genetic evidence.
Within the Kivumys group (subgenus Kivumys), high CYTB differentiation (13.39% p-distance) between the two species represented in our dataset -L. woosnami and L. luteogaster -clearly delimitate them as distinct lineages. Together with L. medicaudatus, whose CYTB sequences were not included in the study, all three species in the Kivumys subgenus have been recorded from overlapping ranges, i.e., in the northeastern and eastern DRC forests and bordering montane forests of the Albertine Rift, with L. woosnami extending into western Burundi, Rwanda, and Uganda [5,6,22,35,38]. A thorough investigation of niche partitioning and other ecomorphological strategies inherent in gene ow and adaptive genetic divergence within the Kivumys subgenus is necessary to clear up their evolutionary history. Such a study would also illuminate the precise nature and limits of their ranges (whether sympatric, syntopic, or parapatric).
sikapusi group and forms a sister relationship with L. sikapusi (separated by 11.42% p-distance), representing a potentially undescribed species. While describing rodents of western and southwestern Guinea, Denys et al. [39] considered their Lophuromys samples from Mankountan, Boffa prefecture, Guinea, as a tentative new species in the sikapusi group due to its unique head and body length, pending morphometric and genetic evidence for a taxonomic description. Our CYTB evidence suggests that the Lophuromys from the Conakry region of Guinea might belong to their undescribed taxa. Still, without more the morphometric evidence to con rm its morphological relationships in the sikapusi group and range limits, we retain a similar provisional species status for L. sp.1.
The utilization of pelage coloration to resolve the systematic grouping of species is rather debatable in the genus Lophuromys. For instance, the inclusion of L. dieterleni (Mt Oku, Cameroon) and L. eisentrauti (Mt Lefo, Cameroon) in the L. avopunctatus group by Musser and Carleton [6] (citing Verheyen et al. [25]) is quite doubtful. Our CYTB tree shows all unspeckled-pelage taxa clusters in the L. sikapusi group (L. angolensis, L. ansorgei, L. huttereri, L. nudicaudus, L. rahmi, L. roseveari, L. sikapusi, and L. sp.1), well distinct from the speckled-pelage L. avopunctatus group. From an ecomorphological outlook, the craniodental relationship between L. dieterleni and L. eisentrauti and the L. avopunctatus group [25,26] might simply be signals of convergent adaptive responses to local environments [40,41], which are taxonomically uninformative without genetic evidence. Then again, the genetic and craniodental a nity of the unspeckled L. pseudosikapusi to the Ethiopian endemics [13] confounds further the overall phylogenetic relationships within and between the speckledpelage (L. avopunctatus group) and unspeckled-pelage (L. sikapusi group) species. From our ndings, assigning clades to either the L. sikapusi group or the L. avopunctatus group is unambiguous based on genetic relationships, and future genetic studies are likely to resolve the L. dieterleni and L. eisentrauti membership.
The assignment of ETHFLAVO clades to corresponding species is a nontrivial task due to the recently clari ed taxonomic accounts of the Ethiopian Lophuromys [15][16][17]. For example, the pairs of highly divergent CYTB clades of L. simensis (L. simensis 1 and L. simensis 2) and L. melanonyx (L. melanonyx 1 and L. melanonyx 2) comprise the multiple haplogroups within the same species due to past mtDNA introgression events. Such introgressions have also been con rmed in L. brunneus, of which we only sampled one haplogroup, summing up to the 12 mtDNA Lophuromys lineages endemic to Ethiopia. Nuclear genomic data support these 12 lineages represent nine species (L. chrysopus, L. melanonyx, L. simensis, L. avopunctatus, L. brunneus, L. pseudosikapusi, L. menageshae, L. chercherensis, and L. brevicaudus) which differ by karyotypes, morphology, and preferred elevation, i.e., types of ecosystems [13,15,16]. Interestingly, mitochondrial introgression, apparently common in the avopunctatus group, was not detected in the rest of the genus Lophuromys, although it should be noted that nuclear genetic data are relatively scarce outside the Ethiopian avopunctatus members.

Species divergence and biogeography
The nested phylogenetic relationships between the ETHFLAVO and NONETHFLAVO conform generally to evolutionary processes speculated previously [13,14]. While our ndings support the Ethiopian highlands as the cradle of the speckled-pelage Lophuromys, the precise nature of their evolutionary radiation, including processes characterizing the observed differentiation between clades remains a matter for speculation, mainly owing to the strong effect of mtDNA on the inferred phylogeny. In any case, it is currently not possible to ascertain whether long-distance dispersal and/or montane-forest bridges promoted the divergence and dispersal of non-Ethiopian avopunctatus members out of the Ethiopian Highlands. Dispersal along a north-south axis, i.e., out of Ethiopia to southern Afromontane Highlands is relatively like that of other montane-adapted rodents [42][43][44] and attributed to montane forest expansion during Pliocene-Pleistocene interglacials.
The timing of the NONETHFLAVO and NONETHFLAVO out-of-Ethiopia dispersals, albeit based on a single mitochondrial locus, coincide with the repeated expansion and contraction/isolation of montane forests and their faunal assemblages during the humid intervals of Pleistocene glacial-interglacial cycles [17,[45][46][47][48]. Within the L. avopunctatus group, these events likely connected the southern Ethiopian Highlands with Albertine Rift montane forests, and Kenyan and Tanzanian Highlands across the currently arid Turkana depression [43,45,[49][50][51]. The L. avopunctatus group is primarily restricted to humid/wet habitats which are currently con ned to montane areas in East Africa. These species could only have dispersed when the East Africa Highlands were connected with similarly suitable habitats. The rst out-of-Ethiopia dispersal by the NONETHFLAVO1 ancestor and consequent range retention in the northern EAMs concur with their prolonged stability that preceded the formation of most of the Kenya Highlands, Tanzanian Highlands, and Albertine Rift montane forests. The split and dispersal of the L. aquilus + L. verhageni clade from L. kilonzoi, the consequent split of L. aquilus from L. verhageni, and the appearance of several clades in the NONETHFLAVO2 subgroup, all happened in the mid-late Pleistocene. This coincides with wet climate periods that made it possible to cross currently dry valleys such as those isolating Mt. Kilimanjaro and Mt. Meru and the Turkana depression in northern Kenya and southern Ethiopia [44,52]. The absence of genetic evidence of this rst dispersal in East African highlands such as the Kenyan Highlands, suggests these mountains served as 'steppingstones'. The ' rst colonizers' -NONETHFLAVO1 -were presumably replaced by the more successful 'second colonizers' -NONETHFLAVO2, such as L. zena in the Kenyan highlands. Whether or not the second colonizers hybridized with the rst ones remains unclear from mtDNA, and genomic analyses should be applied to investigate this possibility.
The Albertine Rift Valley is a crucial biogeographical feature in the radiation of the NONETHFLAVO and is likely an active barrier to gene ow on either side.
The four distinct clades whose ranges are separated by the Albertine Rift (L. dudui, L. stanleyi, L. cf. cinereus, and L. laticeps) suggest that they are not able to cross and have not experienced gene ow since their divergence. While L. stanleyi occurs widely eastward of the Rwenzori Mountains, it does not extend west of the mountains, whereas the range of L. dudui begins in Virunga National Park, and only extends westwards. The Albertine Rift might have been a barrier to L. stanleyi's westward dispersal and L. dudui's eastward dispersal. This hypothesis is also consistent with the occurrence of L. laticeps on the eastern and L. cf. cinereus on the opposite western side of the Albertine Rift around Lake Kivu and Lake Tanganyika, with either presumptively unable to cross. Notably, L. sabuni, which is the only clade whose occurrence span both anks of the Albertine Rift, appear to have dispersed between the Rukwa Rift and Lake Tanganyika and then southwards to Chishimba Falls (Northern Zambia), where it was recently recorded (Sabuni et al. [53]. Other forest rodents have ranges that span the Albertine Rift, unlike observed here for L. dudui, L. stanleyi, L. cf. cinereus, and L. laticeps. For instance, the Malacomys longipes [54] and the Praomys jacksoni [55] occur on both sides of the Albertine Rift Valley.

Morphological variation within the non-Ethiopian avopunctatus members
Most species in the NONETHFLAVO have overlapping craniodental characters in morphospace, making our large dataset of linear measurements and geometric landmarks unreliable as the exclusive evidence to infer species limits. For instance, the range of skull morphology of L. stanleyi and L. zena (both linear and geometric) signi cantly resembles the skull forms of all other clades in the NONETHFLAVO, except L. verhageni and L. aquilus, which unambiguously cluster and have the least overlap with any other species in the group. The L. stanleyi and L. zena clades exemplify a typical systematic problem in the NONETHFLAVO, where morphological evidence cannot classify samples to meaningful species units using taxonomically informative characters. Accounting for phenetic variation in the NONETHFLAVO, beyond their common ancestry, requires more comprehensive genomic analyses to disentangle the underlying ecomorphological processes among species occurring in similar habitats. Without such genomic evidence, the taxonomic accounts of several clades are best not considered reliably resolved when based on linear or geometric morphometrics only.
Divergence dates and biogeographic patterns in the NONETHFLAVO suggest that the drivers of craniodental variation t multiple non-exclusive hypotheses associated with the correlation of ecomorphological divergence with speciation [56]. The relatively recent divergence of most of the clades suggests ecologically-mediated adaptive evolution might not be predominant speciation drivers between congeners in sympatry [57][58][59]. Except for L. aquilus and L. verhageni which are restricted to single mountain ecosystems, all the NONETHFLAVO species appear to have non-specialized niches as they are not restricted to high montane habitats. They are, thus, more likely to exhibit non-specialized morphological traits that are taxonomically uninformative [60]. The treatment of the L. avopunctatus group by Verheyen et al. [12] and Verheyen et al. [14] highlights the use of craniodental and external morphology data to recognize populations as unique species with minimal use of genetic data. However, inter/intraspeci c taxonomic delimitation among rodents often have fewer taxonomically informative stable morphological states, possibly due to nonadaptive and or rapid adaptive radiations [61]. These in uences might hinder a replicable de nition of taxon-speci c phenotypic traits [62][63][64], leading to the subjective interpretation of valid species. While our geometric morphometry appears generally more sensitive at detecting variabilities between clades compared to linear measurements, just like in other cases [65], over-all, both datasets produced virtually similar results.

Taxonomic assessment of the non-Ethiopian avopunctatus members
While most of the OTUs recovered in the L. avopunctatus group represent species currently named, the between-clade genetic distances and CYTB incongruence with morphometric and IRBP gene results raise more taxonomic questions than resolution. For instance, only a few mutations at CYTB separate L. aquilus from L. verhageni (2.8% p-distance) and L. stanleyi from L. zena (2.9% p-distance), which is among the closest between-species CYTB divergences in the L. avopunctatus group. While such low sequence divergence between these sister clades indicates a recent separation of gene pools (at least at mtDNA), it nonetheless, raises concerns about the species' taxonomic validity, suggesting the need to synonymize them in future taxonomic revisions, without more genetic support, especially since no clear diagnostic morphological differences delimit them. Moreover, the IRBP failure to delimitate several distinct mtDNA clades in the NONETHFLAVO might relate to its slow mutation rate which makes it unable to resolve deeper and or short branches among rodents [66,67]. Nevertheless, future taxonomic reassessments of the genus Lophuromys should utilize more comprehensive genomic analysis (such as multiple nuclear loci) which are likely to be more informative in delimitating phylogenetic relationships. Such genomic evidence would also quantify the level of distinctiveness between close relatives that are allopatric such as L. aquilus and L. verhageni and the level/absence of gene ow between parapatric ones such as L. stanleyi and L. zena as is prevalent in the ETHFLAVO [16].
The genetic diversity within lineages such as L. sikapusi and L. chrysopus, for instance, showed the delimited OTUs within them were similarly distinguishable based on CYTB comparably to several clades in the NONETHFLAVO. It appears that taxonomic classi cations of Lophuromys species that is not based on extensive nuclear evidence, should yet be regarded inconclusive (i.e., within the Kivumys group, L. sikapusi group, and NONETHFLAVO).
The NONETHFLAVO1 subgroup -L. aquilus, L. verhageni, and L. kilonzoi The samples from Mt. Kilimanjaro, Mt. Meru, and northeastern Tanzanian Eastern Arc Mountains form distinct monophyletic lineages, representing species currently recognized as valid. L. aquilus was described by True [23] from Mt. Kilimanjaro and con rmed by Verheyen et al. [14] to be the only Lophuromys along the entire elevation gradient. Lophuromys verhageni was described by Verheyen et al. [12] as an endemic of Mt Meru, while L. kilonzoi was described by Verheyen et al. [14] from the Magamba, East Usambara. Perhaps because of fewer informative sites in shorter sequences, L. aquilus, L. verhageni, and L. kilonzoi had a different phylogenetic topology in Verheyen et al. [14]. Our expanded CYTB sampling supports the three species are minimally differentiated, forming a sister clade to one of the haplogroups of L. simensis. The current CYTB phylogeny, therefore, provides a clearer picture of the phylogenetic relationship between L. aquilus, L. verhageni, and L. kilonzoi and their position in the genus. The historical biogeographical reconstruction suggested the colonization and divergence of L. aquilus and L. verhageni resulted from vicariance events that coincide with the Pleistocene climatic oscillations which fragmented humid montane forests in East Africa [45,54,55,68,69]. The savannas separating their current ranges were substantially stable even across glacial cycles in the late Pleistocene [69]. The occurrence of L. aquilus and L. verhageni is also consistent with the endemism of Crocidura newmarki on Mt. Meru [70] and Myosorex zinki on Mt Kilimanjaro [71]. The divergence and dispersal of the NONETHFLAVO1 ancestor coincided with temporary biogeographical contacts between the Ethiopian Highlands and other Afromontane forests in the early Pleistocene [69]. After initial colonization, montane forests were again fragmented by climatic oscillations, in the process facilitating allopatric speciation. Similarly, the patchy distribution of Praomys delectorum across the Eastern Arc Mountains and Southern Montane Forests was probably driven by comparable vicariance events [68]. The higher genetic diversity within L. kilonzoi suggests that it remained in the ancestral range after diverging from the ancestor of L. aquilus + L. verhageni. Similar divergence and diversity patterns were observed for the forest-dependent Praomys delectorum [68], where the MRCA of populations from the Eastern Arc Mountains predated those from Mt. Kilimanjaro and Mt. Meru and correlated with genetic diversity.
The northern part of NONETHFLAVO: L. stanleyi and L. zena The sister clades from the Kenyan and Ugandan highlands, L. zena and L. stanleyi, comprise the northern part of the non-Ethiopian avopunctatus members. Our sample coverage of these two clades was the most comprehensive to date and substantially extend their known ranges. Lophuromys zena [31], thought to be endemic to the higher elevations of central Kenya [12,14], occurs in all the stably humid ecosystems in Kenya, including Loita Hills forests in the southeast of their range to western (Mt. Elgon, Cherangani Hills, and Kakamega Forest) and southwestern (Victoria Basin -Yala Swamp) Kenya. The distribution of L. zena overlaps with L. stanleyi and L. ansorgei in the Kakamega Forest and with L. ansorgei in Yala Swamp. This distribution supports the conclusions of Onditi et al. [72] that L. zena could have been much more widespread in Kenya than previously known [12,14]. The range of L. stanleyi is also much more extensive than previously described. Sabuni et al. [53] extended the range of L. stanleyi (delimited by mtDNA) into northwestern Tanzania beyond its Mt Rwenzori type locality [12], where it was thought to be restricted. Here we provide evidence that L. stanleyi occurs through much of Uganda, spanning southeastern South Sudan and northeastern Uganda forests eastwards to the Kakamega Forest in Kenya (its eastern limit) and south into northern Rwanda.
The 'Karamoja/Uganda gap' [47,73] was not a barrier to the dispersal of L. stanleyi through Uganda to connect the Kenya Highlands and Albertine Rift montane forests, as was the case for the forest-dependent Hylomyscus [47,74]. Generally, the sister relationship of L. zena and L. stanleyi (minimal CYTB divergence) reinforce biogeographic a nities between the Albertine Rift montane forests and the Kenya Highlands [47][48][49]. Furthermore, the occurrence of L. zena and L. stanleyi in both lowland and highland forests suggest a phylogeographic pattern shaped also by an opportunistic ecological strategy, unlike true forest-specialists such as the Hylomyscus denniae and Sylvisorex granti groups that are restricted to high-elevation forests [47,48]. The biogeographies of L. zena and L. stanleyi mirror similar patterns as the more widespread Praomys jacksoni which colonized both montane and lowland forests. However, the absence of a taxonomic structure based on IRBP further reinforces the need to apply genomic analyses, especially in zones of secondary contacts, such as in Kakamega, to shed light on the level of their reproductive isolation and the taxonomic validity.
The southern part of NONETHFLAVO: L. machangui and L. sabuni These two signi cantly supported sister clades correspond to L. machangui and L. sabuni, both described by Verheyen et al. [14] from Mount Rungwe and the Mbizi Mountains (U pa Plateau), respectively. Their sister relationship and late Pleistocene divergence coincide with the split of L. verhageni and L. aquilus, attributable to the late Pleistocene expansion of moist forests that likely enabled them to disperse to the current ranges, whose suitability was later restricted to highland forests. Overall, the distribution of L. machangui and L. sabuni reveal biogeographical trends that both coincide and contrast with other small mammals in the region, suggesting that other taxon-speci c functional traits, such as dispersal ability and habitat speci city versus generality also in uenced their evolutionary radiation. For instance, the distribution of L. machangui suggests the Makambako Gap has not barred its dispersal, similar for other small mammals including Myosorex kihaulei [75], but has barred the dispersal of Praomys delectorum [68] and Otomys lacustris [76]. Within the range of L. sabuni, Kerbis Peterhans et al. [73] recently described two species in the genus Hylomyscus, Hylomyscus stanleyi from Mbizi Forest Reserve and Hylomyscus mpungamachagorum from Mahale National Park, suggesting that the so-called Karema Gap was a barrier to the dispersal of these Hylomyscus species but not to the dispersal of L. sabuni. Overall, the close craniodental and genetic a nity between L. sabuni and L. machangui to each other in comparison with members from other clades in the NONETHFLAVO2 subgroup suggests they have experienced somewhat similar ecological selection resulting in similar ecomorphological characteristics [77]. The craniodental and genetic a nities between L. sabuni and L. machangui also concurs with the oral and faunal a nity between the Southern highlands of the northern end of Lake Malawi and the Mbizi Forest, attributed mostly to the absence of a substantial biogeographical barrier between them. More studies are needed to delineate genetic differentiation across the range of L. machangui and L. sabuni, and detail how isolation by distance and geographical features have impacted dispersal.
The western part of NONETHFLAVO: L. dudui and L. rita The L. dudui clade comprised samples from the northeastern DRC montane highlands of the Albertine Rift -Rwenzori Mountains, westwards to the Kisangani -Bomane -Yaenero areas. This distribution leaves a ca. 480 km sample gap between the eastern limits (Epulu -Tshiabirimu -Ituri) and western limits near Baliko, Boende on the left bank of Congo River. The inclusion of samples from both sides of the Congo River in the L. dudui clade modi es the original description as well as consequent accounts of L. dudui, where it has been thought to be restricted between the right bank of the Congo River and the western foothills of the Albertine Rift mountains [12,14]. The current range of L. dudui resembles that of Praomys mutoni and Praomys jacksoni [55] both of which occupy lowland forests on both banks of Congo River in the Kisangani region [55,78]. Morphologically, L. dudui is easily diagnosable from the nearby NONETH2 members due to its distinctly small skull (Fig. 6, Additional le 2 Table S3), consistent with previous ndings [12,14]. The L. dudui range overlaps with that of L. rita, which was assigned to samples spread over an expansive area in the Congo Basin, spanning southwestern DRC (Kinshasa) to the northeast (Kisangani, left bank of Congo River) and southwards to northwestern Zambia. Although we are unable to make skull comparisons with the holotype, this clade forms a well-de ned mtDNA lineage, probably representing the L. rita described by Dollman [31] from south of Lake Tanganyika in NE Zambia (Mporokoso) and Lufupa River, Katanga, DRC. Despite its expansive range, L. rita appears bound to the central Congo basin by the Congo and Lualaba Rivers, which have likely limited its dispersal, like Praomys minor in the central Congo Basin [55]. Our geographic sampling of L. rita is notably sparse relative to its distribution and more surveys are necessary to resolve the full range and genetic diversity within the L. rita clade. Importantly, a formal taxonomic reassessment is required to validate the morphological relationship of the L. rita clade with the holotype and topotypes.

L. makundii
Specimens attributed to the monophyletic L. makundii derive from the Mount Hanang (type locality) northwards over the Lake Manyara and Ngorongoro crater to Mt Kitumbeine. Several 'unsuitable' dry corridors which currently isolate L. makundii from Eastern Arc Montane forests, Albertine Rift Mountains, Kenyan Highlands, and even the nearby Mount Kilimanjaro and Mount Meru seem to have impacted its dispersal after the initial colonization event. However, the occurrence of Crocidura montis, Crocidura hildergadea, Otomys angoniensis, Grammomys dolichurus/macmillani, Graphiurus murinus, and Praomys taitae in similar habitats as L. makundii in the north-central Tanzania region [53,79] suggest that its biogeographical a liation to other Eastern Afromontane forests in the region is recent. The relatively isolated range of L. makundii likely imposed a more rigid barrier to genetic exchange with other lineages after divergence [80], and might explain why it is the only other clade in the NONETHFLAVO, besides L. kilonzoi, that retains monophyly in the IRBP tree. Still, the minimum divergence time and possible sister relationship to either L. dudui or L. laticeps show that it is more closely a liated to the Albertine Rift than the NONETHFLAVO1 range. As such, L. makundii probably colonized its current range when moist forests connected the currently isolated volcanic mountains during the late Pleistocene climate uctuations.
L. cf. cinereus and L. laticeps The L. cf. cinereus samples overlap the Kahuzi-Biega National Park locality from where Dieterlen and Gelmroth [28] described L. cinereus. Following the initial proposal by Dieterlen [11] that the external and craniodental distinctness used by Dieterlen and Gelmroth [28] to describe L. cinereus were, in fact, morphotypes of L. laticeps, there has since been no formal taxonomic reassessment of its validity, with the foundational references [6,12] maintaining its synonymy to L. laticeps. Our mtDNA, nuclear (IRBP), and craniodental tests showed similar differences between the L. cf. cinereus and L. laticeps clades, comparable to the distances within and between other NONETHFLAVO2 clades, including the sister clade, L. rita. The L. cf. cinereus skulls overlapped most with L. laticeps, L. dudui, and L. stanleyi, consistent with the earlier rationale for its synonymy [6,12]. A formal taxonomic revision of L. cf. cinereus, is needed to validate and update its distribution, and genetic and phenetic relationship to other non-Ethiopian L. avopunctatus members. Such a revision would importantly update the occurrence of L. cinereus (herein as L. cf. cinereus), which was perceived restricted to the type locality [28,81], to extend from Kahuzi-Biega NP to the Itombwe Massif and southwards ca. 300 kilometers to Mt. Kabobo -on the western shore of Lake Tanganyika. Thomas and Wroughton [29] considered L. laticeps as a morphologically unique lineage among its close relatives allied to L. aquilus [23] due to a broader lower braincase and shorter palatal foramina. In agreement, our L. Laticeps skulls had the broadest BBC and one of the shortest PPL in NONETHFLAVO2 dataset (Additional le 2 Table  S3). The L. laticeps clade is also genetically well-differentiated, comparably to close relatives -L. cf. cinereus, L. stanleyi, and L. laticeps. There is a need to formally reassess the taxonomy of L. cf. cinereus and L. laticeps, to clarify and update their distinctness from other lineages in the NONETHFLAVO, not to be synonymized under L. aquilus.
L. margarettae, L. rubecula, and L. major No genetic OTUs could be matched to L. margarettae, L. rubecula, or L. major, despite sampling from their respective ranges -Mathews Range, Mount Elgon, and proximity of Ubangi River. L. margarettae was described by Heller [82]  Similarly, L. rubecula described by Dollman [27] is another species we are unable to con rm without new material. Our Mt. Elgon samples cluster genetically and craniodentally with L. zena. However, we lacked samples from other parts of the Mt Elgon ecosystem, without which we cannot dismiss L. rubecula's occurrence or its validity. Future surveys of Mt Elgon should employ elevational strati ed sampling transects on the Kenya and Uganda sides to substantiate the occurrence limits (or the absence thereof) of L. rubecula.
Finally, we were also unable to verify the validity of L. major, which was described by Thomas and Wroughton [29] from the Bwanda area, Ubangi River, DRC. The ranges of the presupposed nearest congeners -L. dudui and L. rita are considerably south of its type locality and without new material from the area, we cannot verify the validity of L. major or approximate its relationship to other species in the Lophuromys genus.

Conclusion
Despite being one of the most widely occurring and abundant rodents in east-central and East African montane and lowland rain forest habitats, the taxonomy and biogeography of the non-Ethiopian L. avopunctatus members (NONETHFLAVO) remain poorly understood. Our utilization of the CYTB gene to reconstruct the phylogenetic relationships of the genus Lophuromys and combined mitochondrial and nuclear genes and morphometrics (geometric and linear characters) to analyze the systematics of the NONETHFLAVO revealed a substantial number of novel ndings. Most of the species recognized previously based on morphology are supported as well geographically structured mtDNA lineages but lack stable informative craniodental characters capable of reliably assigning samples to putative species units a priori. The NONETHFLAVO colonized its current range over two independent dispersal events out of Ethiopia in the early Pleistocene, with the two resulting subgroups remaining respectively monophyletic but nested in the Ethiopian avopunctatus members. Despite our study providing a comprehensive scenario for the evolution, phylogeography, and genetic differentiation of the NONETHFLAVO, a formal mitochondrial-nuclear phylogenetic incompatibilities, as done recently for the Ethiopian L. avopunctatus members [16]. Ultimately, such a comprehensive genomic phylogenetic approach, even in the absence of craniodental data, is likely to reliably delimit the unique population pools corresponding to valid species and the resolution of species groups. Currently, NONETHFLAVO occurs in ecosystems with stable annual precipitation which are susceptible to changing climates, meaning that ongoing climate uctuations, warming climates, and continued habitat fragmentation are likely to fragment further their distributions, resulting in substantial range shifts and or loss of habitat. This would probably drive divergent ecomorphological adaptations between sympatric and parapatric close relatives, with taxonomic implications that are essential from a conservation point of view. The collection and handling of animals adhered to the wildlife research regulations of the respective countries. New eldwork in Kenya conducted by joint teams from the National Museums of Kenya and Kunming Institute of Zoology was permitted by Kenya Wildlife Service and Kenya Forest Service, whose o ces also provided key logistical support. Other new material was obtained from The Field Museum of Natural History (Chicago, USA) and we are grateful to the curators for facilitating access to study their collections.

Consent for publication
Not applicable Availability of data and materials All data generated or analyzed during this study are included in this article and its supplementary information les.

Competing interests
The authors declare that they have no competing interests