Rapid ecological specialization despite constant population sizes

Background The bushbuck, Tragelaphus scriptus, is a widespread and ecologically diverse ungulate species complex within the spiral-horned antelopes. This species was recently found to consist of two genetically divergent but monophyletic lineages, which are paraphyletic at mitochondrial (mt)DNA owing to an ancient interspecific hybridization event. The Scriptus lineage (T. s. scriptus) inhabits the north-western half of the African continent while Sylvaticus (T. s. sylvaticus) is found in the south-eastern half. Here we test hypotheses of historical demography and adaptation in bushbuck using a higher-resolution framework, with four nuclear (MGF, PRKCI, SPTBN, and THY) and three new mitochondrial markers (cytochrome b, 12S rRNA, and 16S rRNA). Methods Genealogies were reconstructed for the mitochondrial and nuclear data sets, with the latter dated using fossil calibration points. We also inferred the demographic history of Scriptus and Sylvaticus using coalescent-based methods. To obtain an overview of the origins and ancestral colonisation routes of ancestral bushbuck sequences across geographic space, we conducted discrete Bayesian phylogeographic and statistical dispersal-vicariance analyses on our nuclear DNA data set. Results Both nuclear DNA and mtDNA support previous findings of two genetically divergent Sylvaticus and Scriptus lineages. The three mtDNA loci confirmed 15 of the previously defined haplogroups, including those with convergent phenotypes. However, the nuclear tree showed less phylogenetic resolution at the more derived parts of the genealogy, possibly due to incomplete lineage sorting of the slower evolving nuclear genome. The only exception to this was the montane Menelik’s bushbuck (Sylvaticus) of the Ethiopian highlands, which formed a monophyletic group at three of four nuclear DNA loci. We dated the coalescence of the two lineages to a common ancestor ∼2.54 million years ago. Both marker sets revealed similar demographic histories of constant population size over time. We show that the bushbuck likely originated in East Africa, with Scriptus dispersing to colonise suitable habitats west of the African Rift and Sylvaticus radiating from east of the Rift into southern Africa via a series of mainly vicariance events. Discussion Despite lower levels of genetic structure at nuclear loci, we confirmed the independent evolution of the Menelik’s bushbuck relative to the phenotypically similar montane bushbuck in East Africa, adding further weight to previous suggestions of convergent evolution within the bushbuck complex. Perhaps the most surprising result of our analysis was that both Scriptus and Sylvaticus populations remained relatively constant throughout the Pleistocene, which is remarkable given that this was a period of major climatic and tectonic change in Africa, and responsible for driving the evolution of much of the continent’s extant large mammalian diversity.


INTRODUCTION
The bushbuck (Tragelaphus scriptus) is a well-known, highly diverse species complex of spiral-horned antelopes, inhabiting most of sub-Saharan Africa (Moodley et al., 2009;Hassanin et al., 2012). This species complex is unique, being the most widespread and ecologically diverse of any bovid species and occurring in approximately 73% of the total land area of sub-Saharan Africa. Across this vast and heterogeneous region, bushbuck can be found in most habitat types (Moodley & Bruford, 2007) from forested to xeric zones and ranging in altitude from sea-level to 4,000 m.
Phenotypic diversity among bushbuck populations is unprecedented among the bovids. Although the number of subspecies described varies, up to twenty four have been recognised by a single author (Lydekker, 1914;Allen, 1939). The complex can be subdivided into two divergent morphological groups which inhabit the western and northern (T. s. scriptus group) and eastern and southern (T. s. sylvaticus group) parts of the species range ( Fig. 1), hereafter Scriptus and Sylvaticus for ease of reference. Scriptus is smaller and less dimorphic, but it possesses a heavily striped white harness-like pattern, whereas most populations of the larger Sylvaticus have little to no striping at all. Although known to favour areas of thick cover wherever they occur, bushbuck do not inhabit the dense rainforest of the Congo basin, preferring the mosaic landscapes at its fringe. The two groups are therefore separated in the west and south by the Lower Congo valley and the Congo basin respectively, but in eastern Africa Scriptus and Sylvaticus come into secondary contact from the northern end of the Albertine rift along the Imatong and Didinga Mountains of South Sudan following the rift into the Ethiopian Highlands (white arrows, Fig. 1). Within this zone of contact, the phenotypic integrity of each form may be maintained through habitat preference; the Scriptus form inhabits the low lands while the large, dark and heavy-coated Sylvaticus montane bushbuck inhabits the high altitude forests, although evidence of gene flow has been observed (Moodley & Bruford, 2007).
Scriptus and Sylvaticus can also be separated genetically. Initial mitochondrial (mt)DNA studies divided the bushbuck into Scriptus and Sylvaticus, but with either lineage more closely related to other Tragelaphus species than to each other (Moodley & Bruford, 2007;Moodley et al., 2009). This mtDNA paraphyly prompted some authors to regard Scriptus and Sylvaticus as independent species (Moodley et al., 2009;Hassanin et al., 2012;Hassanin et al., 2018), possibly evolving through convergent evolution (Wronski & Moodley, 2009). However, a recent analysis of nuclear DNA among the spiral-horned antelopes showed that Scriptus and Sylvaticus, although genetically divergent, are reciprocally monophyletic and so the bushbuck may henceforth be considered a single species. Paraphyletic Scriptus Figure 1 The land cover of Africa reconstructed from remotely sensed data (redrawn from Mayaux et al., 2004). The geographical distribution of sampling localities included in the present study are shown on the map. Taxa are plotted as dots and designated either blue for Scriptus or black for Sylvaticus. Samples are numbered according to Table 1 and Sylvaticus mtDNA lineages thus arose through an ancient interspecific hybridization event (Rakotoarivelo et al., in press). Across its range, the bushbuck was highly structured into 23 phylogenetically distinct haplogroups (Scriptus 8; Sylvaticus 15), each with differing levels of ecological specialization. Among the more specialized haplogroups, the montane (T. s. meneliki, T. s. powelli, T. s. barkeri and T. s. delamerei), and more xeric-adapted (T. s. decula, T. s. dodingae, T. s. fasciatus1, T. s. fasciatus2 and T. s. roualeyni) appear to have evolved more than once through convergent evolution (Moodley & Bruford, 2007). Much of the mtDNA variation in the complex is structured according to ecoregion (Olson et al., 2001), suggesting local ecological conditions as a driver for the evolution of specialization. Ecological conditions are in turn driven by a combination of local geology and an oscillating Pleistocene paleoclimate (Vrba, 1995;Bobe & Behrensmeyer, 2004;Hernandez Fernández & Vrba, 2005). However, where the species evolved and its subsequent routes of colonization and diversification are still unknown. Despite the research potential of this system, only mtDNA data have been generated for this species to date. Not only is the mitochondrial genome a single locus, it is also maternally inherited so mtDNA structure may not be representative of nuclear DNA structure in species with sex biases in dispersal/philopatry. Genetic drift is also more effective in sorting non-segregating mtDNA lineages as their effective population size is approximately four times smaller than segregating nuclear DNA. Therefore, whether the nuclear genome is structured similarly, or even whether Scriptus and Sylvaticus constitute different nuclear lineages, is unknown. Furthermore, demographic analyses that may evidence population responses to paleo-environmental conditions and a spatially-informed phylogeographic analysis of origins and colonisation routes have never been carried out.
To test the hypotheses of variation, structure and potential adaptation purported by previous mtDNA work, we sequenced representative bushbuck from across the species range using a higher-resolution multilocus framework of four nuclear introns, complemented by three further mtDNA markers. We further reconstructed both the demographic and phylogeographic histories of the bushbuck complex using this new data set to shed further light on the evolution of this species.

Taxon sampling
A total of 27 bushbuck individuals (excluding outgroups) were included in this study. Samples sourced previously by Moodley & Bruford (2007) were re-extracted and representatives of all 23 mtDNA haplogroups were selected ( Fig. 1; Table 1). As outgroups, we used both the distantly related Bos taurus as well as the most closely related lesser kudu (Tragelaphus imberbis) to root trees in several of the phylogenetic analyses.

DNA sequencing
Four nuclear intron DNA markers (MGF-mast cell growth factor, PRKCI-protein-kinase CI, B-spectrin non-erythrocytic 1-SPTBN, and THY-thyrotropin) were amplified and sequenced in the 27 individuals above using previously published primers and methodology (Matthee et al., 2001). Additionally mtDNA sequences were amplified and sequenced from three mtDNA cytochrome b (Cyt b), 12S rRNA, and 16S rRNA (for mtDNA PCR and primer details see Arnason, Gullberg & Widegren, 1993;Simonsen, Siegismund & Arctander, 1998). In order for downstream comparison of summary statistics, the same individuals were sequenced for each locus. Sequences from each gene were first aligned using ClustalW (Thompson, Higgins & Gibson, 1994) as implemented in BioEdit (Hall, 1999), using default settings and thereafter manually to optimize homology. All heterozygous sites in the nuclear DNA were coded using the appropriate IUB code. Model selection for the best fitting substitution model for each gene was conducted in jModelTest (Posada, 2008;Darriba et al., 2012) under the Bayesian information criterion, which was preferred over the Akaike information criterion, to guard against over parameterization by averaging the likelihood over all included parameters.

Analysis of genetic diversity and positive selection
The number of variable sites, number of parsimony informative sites and nucleotide frequencies were estimated for both mtDNA and nuclear DNA separately in MEGA 7 (Kumar, Stecher & Tamura, 2016). Further, for each locus we calculated standard diversity statistics in DnaSP 5.0 (Librado & Rozas, 2009). These included: the number of polymorphic sites (s), number of haplotypes, haplotype diversity (Hd), nucleotide diversity (π), and average number of pairwise differences per sequence (k). Summary statistics were also calculated for the total data and for each major clade inferred form phylogenetic analyses.
We used several analyses to test each of our seven loci for neutrality. The McDonald and Kreitman test (MKT) was used to detect signatures of selection and measure the amount of adaptive evolution within a species at the molecular level. Under this test, a neutrality index (NI) quantifies the direction of departure from neutrality, comparing the ratio of non-synonymous to synonymous variation between species (Dn/Ds) with the ratio of non-synonymous to synonymous variation within species (Pn/Ps). NI was calculated using the Standard and Generalized McDonald-Kreitman Test (MKT; Egea, Casillas & Barbadilla, 2008) website. Because silent mutations are neutral, a neutrality index lower than 1 (i.e., NI <1) indicates an excess of non-silent divergence, which occurs when positive selection is at work in the population. When positive selection is acting on the species, natural selection favors a specific phenotype over other phenotypes, and the favored phenotype begins to go to fixation in the species as the allele frequency for that phenotype increases (Biswas & Akey, 2006). Furthermore, we used the coalescent parameters Tajima's D (Tajima, 1989) and Fu 's Fs (Fu, 1997) to test for departures from the neutral theory and these were calculated in DnaSP v5.

Phylogenetic analyses
Phylogenetic reconstruction was performed using both maximum likelihood (ML) and Bayesian approaches using the software Garli 2.0 (Zwickl, 2006) andBEAST v2.4.5 (Bouckaert et al., 2014) respectively. The total data matrix was partitioned by gene, with the parameters of nucleotide substitution models (12S rRNA −HKY + I + G, 16S −HKY, Cyt b-HKY + I, MGF −TIM1 + I, PRKCI −HKY, SPTBN −HKY, THY −TIM1ef + I) and unlinked across partitions. Each ML analysis was initiated from a random starting tree, with nodal support assessed using 1,000 bootstrap replicates. A 50% majority rule consensus tree was constructed using the CONSENSE program in the PHYLIP package (Felsenstein, 2005). Using BEAST, five independent runs of 1 billion generations each were performed; each run consisted of four Monte Carlo Markov chains (MCMC), with topologies sampled every 100,000 generations. The program Tracer 1.6 (Rambaut et al., 2014) was used to determine that the effective sample size (ESS) had reached >200 for all parameters. In each simulation the first 20% of generations were discarded as burn-in. Genealogies were also reconstructed for the nuclear and mitochondrial data sets and for each gene independently using the same MCMC parameters.

Molecular dating
We dated our nuclear phylogeny, since the mtDNA of bushbuck are paraphyletically related (Moodley et al., 2009), and so mitochondrial branch lengths may be upwardly biased. Multiple fossil calibration points were used to scale nodal depth estimation. We calibrated the bushbuck divergence based on the earliest appearance of T. scriptus s.l. in the fossil record in Kenya (Leake & Harris, 2003) and Ethiopia (Kalb et al., 1982) as early as 3.9 Mya and a minimum age of constraint of 2.58 Mya as suggested by Hassanin & Douzery (1999). An exponential distribution was used with a 2.5% probability quantile set at the age of the fossil with hard bound at the youngest bound and a soft maximum bound, beyond which it is unlikely that the divergence actually occurred. Our last calibration point constrained the evolution of the tribe Tragelaphinii 5.72 Mya (95% probability, 4.7-6.7 Mya;Deino et al., 2002). In the latter case, a normal distribution was used allowing for the actual node age to be equally younger or older than the fossil record. Phylogenetic relationships and divergence times were estimated using an uncorrelated relaxed lognormal Bayesian molecular clock approach in BEAST v. 2.4.5 software (Bouckaert et al., 2014). A Yule speciation process was applied to the tree inference through the MCMC (Markov chain Monte Carlo) with a random starting tree. All other parameters were the same as in previous analysis.

Inferring historical demography
In addition to Tajima's D and Fu's Fs, which may be used to infer demography in neutrally evolving loci, demographic changes in both clades were also inferred from the observed mismatch distribution for each of the populations, calculating the raggedness index (R2) according to the population expansion model in DnaSP (Librado & Rozas, 2009). This measure quantifies the smoothness of the observed mismatch distribution, with lower raggedness characterizing a population that experienced a sudden expansion, whereas higher raggedness values suggest stationary or bottlenecked populations (Harpending et al., 1993;Harpending, 1994). Lastly, changes in effective population size were inferred using Bayesian Skyline Plots (BSP: Drummond et al., 2005). These plots utilize the coalescent properties of gene trees to plot population size changes over time, however, the inferred population sizes could potentially be biased downwards (population decline) if the sample set is significantly genetically structured (Ho & Shapiro, 2011;Heller, Chikhi & Siegismund, 2013). To account for biases due to genetic structure, we divided the data into Scriptus and Sylvaticus groups and reconstructed their demographic histories separately using BEAST (Bouckaert et al., 2014). In order to incorporate stochastic differences between gene genealogies in the estimation of population parameters, we constructed multi-locus Extended Bayesian Skyline Plots (EBSP;Heled & Drummond, 2008) for each clade. In addition, EBSP estimates posterior probabilities for the number of population size change events. A mitochondrial divergence rate of 0.056 per million years was used (Arbogast & Slowinski, 1998) as well as appropriate inheritance scalars were used to account for potential difference in effective population size between mtDNA and nuclear DNA. The lengths of the MCMC chains were set to 1 billion to achieve effective sample sizes (ESS) and proper mixing of Markov chains.

Bayesian phylogeographic reconstruction
We attempted to reconstruct the phylogeographic history of two major clades of the bushbuck complex using our nuclear DNA data set. To do this, we employed the spatial diffusion approach under a Bayesian discrete phylogeographic framework in BEAST 1. 8.4 (Lemey et al., 2009;Drummond et al., 2012). Five independent runs of 1 billion generations each were performed; each run consisted of four Monte Carlo Markov chains (MCMC), with topologies sampled every 100,000 generations. We used three geographical states corresponding to the continental regions where both lineages are present: west (W), east (E), and south (S). These phylogeographic analyses were run under a constant-size coalescent model, with molecular clock parameterised as described above and with a random starting tree as tree model. Bayesian Stochastic Search Variable Selection (BSSVS) was used to identify those rates (colonization routes) that were frequently invoked to explain the diffusion process (Lemey et al., 2009). The maximum clade credibility (MCC) tree was computed and annotated using the BEAST module TreeAnnotator v1.8.4 (Drummond et al., 2012). We then used SpreaD3 v0.9.6 (Bielejec et al., 2016; https://github.com/phylogeography/SpreaD3) to analyze and visualize the spatial diffusion incorporated in our Bayesian phylogeographic reconstruction. This was done by mapping the location-annotated MCC tree with the 95% highest posterior density (HPD) of node locations which was then export as a keyhole markup language (KML) file for animation of the spatial diffusion in virtual globe software. The final results were overlaid onto a base map of Africa. We also ran a further ancestral reconstruction analyses on our nuclear DNA data set using a statistical dispersal-vicariance model (S-DIVA) in the package RASP 4.0 (Reconstruct Ancestral State in Phylogenies; Yu, Harris & He, 2010;Yu et al., 2015). In this analysis we used the same West, East and South geographical regions and loaded the posterior sample of 10,001 trees previously produced in BEAST v. 2.4.5 software (Bouckaert et al., 2014). Only the most likely reconstruction was considered for each node.

Genetic variation and its relationship to taxonomy and biogeography
To test whether nuclear DNA supported the hypothesis that ecology has driven genetic diversification in this complex (Moodley & Bruford, 2007), we tested the fit of a comprehensive biogeographic model (Olson et al., 2001) to the nuclear DNA data, relative to that of taxonomic and geographic models using a multiple regression on genetic distance matrices (MRM), implemented in DISTLM (Anderson, 2004). MRM involves a multiple regression of a response matrix on any number of explanatory matrices, where each matrix contains distances or similarities. Pair-wise genetic distances of nuclear DNA data between all 27 samples was used as the response matrix. The MRM method also allows the use of covariables to assess a models conditional effect on that of explanatory matrices. We defined the basic units for the taxonomy model relative to the proposed phenotypic classification of the bushbuck based on the combined classifications of Grubb-Best (Best, 1962;Grubb, 1985) used in Moodley & Bruford (2007) and also the recently published scheme of Groves & Grubb (2011). It should be note that the regression tests employed here test the taxonomic partitions in the data, and not whether these partitions comprise species-, subspecies-or population-level entities. A matrix of geographic coordinates (latitude and longitude) was included as a covariable to assess the possible the effect of isolation-by-distance (IBD) on the model being tested. In a wide-ranging species, IBD Notes. S, number of polymorphic sites; π , nucleotide diversity; h, number of haplotypes; Hd, haplotype diversity; k, average number of nucleotide differences; S/k, expansion coefficient. Statistically significant results were indicated by asterisks: *P < 0.05, **P < 0.01. may significantly influence genetic structure due to the geographic distance separating the widely distributed sampling locations. MRM method allows the quantification of this effect, conditional on that of biogeography and taxonomy.

RESULTS
This study generated a total DNA sequence alignment of 4,676 bp, of which ingroup taxa accounted for 353 segregating sites. Nuclear introns were less diverse (2,596 bp, 26 segregating sites) than mitochondrial genes (2,080 bp, 353 segregating sites, see Table 2). All DNA sequences were found to be evolving neutrally (MKT: χ 2 P > 0.1).

Structure and divergence
Phylogenetic analyses of mitochondrial ( Fig. 2A) and nuclear (Fig. 2B) multilocus alignments yielded highly concordant ML topologies. Both marker sets recovered two well supported Scriptus and Sylvaticus lineages, although the level of phylogenetic resolution was much higher for mtDNA, recovering the general topology originally observed by Moodley & Bruford (2007), despite much smaller sample sizes. By contrast nuclear introns identified the lineage of the Kidepo bushbuck (T. s. dodingae) as well as a Nile-Abyssinian (T. s. bor-T. s. decula) bushbuck clade within Scriptus. The Sylvaticus clade was also less

Demographic analyses
We found both Fu's Fs and Tajima's D indices to be slightly negative among nuclear and mitochondrial loci, for both Scriptus and Sylvaticus (Table 3). However, only locus SPTBN1 returned statistically significant indices, allowing a rejection of the neutrality/constant population size null hypothesis at the species level. Furthermore, the frequencies of pair-wise differences within each population were also consistent with a null hypothesis of constant population size, with non-significant raggedness indices (R2) for all mismatch distributions (Table 3)

Bayesian phylogeographic reconstruction
We used both discrete Bayesian phylogeography and statistical dispersal-vicariance approaches to reconstruct patterns of spatial dispersal and the ancestral location for the origin of the species complex. Within Scriptus, both analyses separated a well-supported T. s. dodingae-T. s. decula clade in the east, from bushbuck inhabiting regions across the Nile and further west (including the Nile bushbuck, T. s. bor, Figs. 5A and 6A). Sylvaticus also comprised significant phylogeographic structuring, with Menelik's bushbuck most divergent, and other more derived lineages separated into eastern and southern groups (Fig. 5B) or into different groups with histories of either dispersal or vicariance (Fig. 6B). Both approaches identified East Africa as the most likely ancestral location for the origin of the bushbuck radiation. From this origin, Bayesian phylogeography invoked dispersal events in a westward direction for Scriptus and in a southward direction for Sylvaticus, both events occurring on either side of the Congo basin. On the other hand, S-DIVA analysis allowed for the possibility of vicariance, rather than dispersal, as an explanation for nuclear DNA spatial branching patterns. According to this analysis, the initial split from an ancestral Scriptus was a westward dispersal of T. s. bor into central Africa, followed by vicariance that separated T. s. decula from T. s. dodingae, and a subsequent secondary dispersal into West Africa. In contrast, the initial stages of the Sylvaticus radiation into southern Africa are all characterized by vicariance events, with dispersal only invoked for more derived lineages around the Great Lakes region.

Ecological adaptation
MRM analysis revealed that biogeography explained a significant 95% of the nuclear genetic variation within the species complex (Table 4). Taxonomic designation and geographic distance accounted respectively for 88% and 26% of the variation, with the more recent taxonomy of Groves & Grubb (2011) outperforming previously used schemes. Under the conditional influence of isolation by distance, both biogeographic and taxonomic models account for 41% and 65% of the genetic variation respectively.

Patterns of genetic diversity
Nuclear genetic diversity was moderate across the species complex. However, mtDNA diversity was exceptionally high, with only a handful of studies reflecting similar levels (Arctander, Johansen & Coutellec-Vret, 1999;Smitz et al., 2013). That mtDNA diversity was higher than nuclear DNA is expected given the differences in mutation rates between the two sets of loci (Nei & Kumar, 2000), however, particularly high mtDNA diversity in bushbuck also reflect a Pliocene mtDNA introgression event (Hassanin et al., 2018), where Scriptus obtained a nyala-like mitochondrial genome. Within each lineage, the higher diversity of Sylvaticus at both nuclear and mtDNA reflects a slightly earlier coalescence time relative to Scriptus (Fig. 3).

Origins, divergence and secondary contact
Fossil records from the mid-Pliocene (approximately 3.9 Mya) of proto-bushbuck are known from several sites in eastern and southern Africa. T. scriptus remains were recovered in Ethiopia (Kalb et al., 1982) and Kenya (Harris et al., 1988;Leake & Harris, 2003). We observed a more recent diversification of Sylvaticus and Scriptus lineages. Since these fossils predate the estimated divergence within the bushbuck, they suggest a possible ancestral origin from north-east Africa. This is indeed the inference from our Bayesian phylogeography and S-DIVA reconstructions, supporting an origin for the species in East Africa. Until the late Pliocene, east Africa was densely forested habitat (Partridge, Wood & DeMenocal, 1995;Reed, 1997), supporting the idea that ancestral bushbuck were both forest dwelling and used its peculiar harnessed striping pattern as an adaptation for camouflage in closed habitats (Moodley & Bruford, 2007). There is some evidence that striping patterns among other bovids are also associated with living in forest habitat (Stoner, Caro & Graham, 2003). The past 3-2 Mya has seen a major paleoclimatic shift that led to the expansion of grassland habitats in Africa, consequently inducing a drastic change in ungulate community structure, specifically in north-east Africa (Bobe & Behrensmeyer, 2004;Hernandez Fernández & Vrba, 2006;Trauth et al., 2007). This also coincided with major geomorphological processes along the Gregory and Albertine Rifts (Vrba, 1995;Reed, 1997). The combination of paleoclimatic shifts and tectonic uplift have shaped the phylogeography of terrestrial African vertebrates (Flagstad et al., 2001;Trauth et al., 2007;Lorenzen et al., 2010;Voelker, Outlaw & Bowie, 2010;Faulkes et al., 2011;Barlow et al., 2013;Jacobs et al., 2013). The Scriptus-Sylvaticus divergence can also be traced back to this time, and their extant distributions on either side of the Rift Valley (Fig. 1) suggest vicariance of the two lineages, on the basis of the major tectonic uplift events along the East African Rift system. Statistical dispersal-vicariance analysis also suggest that a dispersal-only view of branching events within Scriptus and Sylvaticus from their East African origin may be too simplistic. The evolutionary history of Sylvaticus is predominated by vicariance events, which may help explain why phenotypic diversity is higher in this lineage compared to Scriptus, where a history of mainly dispersal was invoked.  Since divergence, Scriptus and Sylvaticus appear to have remained geographically isolated, however, gene flow between the two cannot be discounted. Although mitochondrial and nuclear multilocus haplotypes were not shared between Scriptus and Sylvaticus, the most common allele at nuclear genes PRKC1, SPTBN1 and THY and the mitochondrial 12S rRNA were shared among samples of both lineages. Shared alleles may indicate polymorphisms that were present in an ancestral bushbuck population, but they may also indicate post-divergence gene flow between Scriptus and Sylvaticus. A further analysis with whole genome sequences may yet shed further light on the role of introgression in the evolution of this species complex.

A stable Pleistocene demographic history
Both bushbuck lineages appear to have been demographically stable through the mid to late Pleistocene (Table 3, Fig. 4), despite most of the diversity within each lineage having evolved during this time. This is a surprising result, as the Pleistocene is known for its dramatic climatic fluctuations. Ungulate population sizes are inherently linked with climate change over evolutionary timescales (Lorenzen et al., 2011), and the distributions of herbivores would presumably have shifted in accordance with vegetation change. In sub-Saharan Africa, Pleistocene population expansions of large mammals such as the kob (Birungi & Arctander, 2000), Jackson's hartebeest (Flagstad et al., 2001), Cape buffalo (Van Hooft, Groen & Prins, 2002;Smitz et al., 2013), hippopotamus (Okello et al., 2005) and lion (Barnett et al., 2014) tend to corroborate this view. Pleistocene demographic contractions such as that of the brown hyaena (Westbury et al., 2018) occur less commonly among African mammals, with most declines taking place during the Holocene, as observe for drill baboons (Ting et al., 2012) and white rhinoceros (Moodley et al., 2018). Yet, during the same period, bushbuck mitochondrial and nuclear DNA shows little evidence of demographic change since the Scriptus-Sylvaticus divergence. It is possible that bushbuck, being highly adaptable and ubiquitous generalists, are less demographically affected by climatic fluctuations, and that evolutionary change occurs more through vicariance than population size changes.

Rapid ecological specialization
Demographic stability also appears to be at odds with high levels of variation observed both morphologically and genetically. The extant genetic diversity in both Sylvaticus and Scriptus was generated in the late Pleistocene, <1 Mya, but with most divergences occurring within the last 0.5 Mya. Much of this diversity is reflected in mitochondrial DNA ( Fig. 2A), and has been described previously (Moodley & Bruford, 2007). Although, fewer divergence events were identified with nuclear intron sequences, a large proportion of the nuclear sequence diversity could be attributed to biogeography, even when conditioned on geography (Table 4). This lends strong support to the hypothesis that local ecology has helped shape the structure of genetic diversity in this species.
By dating our nuclear tree we were also able to estimate a more realistic timeframe for the onset of divergence events in the species complex, compared to the mtDNA-based timeframes reported by Moodley & Bruford (2007). Within Sylvaticus, Menelik's bushbuck (T. s. meneliki) was first to diverge into cooler habitats of the Ethiopian massif. Larger size, a darker and thicker coat are typical of several mammalian montane forms (egs. Red squirrel, Paraxerus palliates; Saola, Pseudoryx nghetinhensis). Bergman's rule predicts an increase in size among colder-adapted species (Bergmann, 1847;Freckleton, Harvey & Pagel, 2003;Clauss et al., 2013), whereas darker and thicker coats help in thermoregulation (Caro, 2005;Clusella-Trullas et al., 2008;Mills & Hes, 1997;Amy & Kunz, 2012). The early differentiation of montane Menelik's bushbuck, and the more recent evolution of other montane bushbuck (eg. T. s. barkeri, T. s. delamerei) strengthens evidence for the independent convergence of the montane phenotype among Sylvaticus bushbuck.
The Somali bushbuck (T. s. fasciatus) is also large in size and is able to survive deep into the xeric interior of the Horn of Africa along the watercourses of the Wabi Shebelle and the Juba River. This population comprises two paraphyletic mtDNA lineages ( Fig. 2A) and independent nuclear lineages (Fig. 2B), suggesting the bushbuck colonized the Somali arid zone through two migration or range expansion events of different coastal bushbuck populations from the south.
Within Scriptus, the Nile-Abyssinian bushbuck (T. s. bor-T. s. decula) clade diverged into the more open, drier habitats of the mosaic region on the fringes of the Sahel. This is reflected in phenotype, as most Scriptus populations are strikingly patterned with the typical bushbuck ''harness'', striping is reduced in those Scriptus populations in more open habitats such as T. s. bor, T. s. decula and T. s. dodingae. There is also a suggestion of reduced patterning among Sylvaticus bushbuck. Although much less strikingly coloured, individuals in some Sylvaticus populations such as the Chobe bushbuck (T. s. ornatus) and the Ituri bushbuck (T. s. dianae) may be more heavily patterned with vertical and horizontal stripes and spots. However, such individuals become rarer in populations to the south where habitats are drier and more open. A similar loss of patterning occurs across the north-south range of the plains zebra, which is also suggested to be in response to open drier environments (Rau, 1978;Leonard et al., 2005).

CONCLUSIONS
In the present study, we sequenced mitochondrial and nuclear DNA 27 individuals representing the range of genetically distinct haplogroups previously described within the bushbuck complex. Phylogenetic congruence was observed between mitochondrial and nuclear markers, both identifying two lineages that diverged in the late Pliocene (Scriptus and Sylvaticus), with further diversification into more specialised groupings during the Pleistocene. Although climatic upheaval during the Pleistocene may have promoted one of the more astonishing arrays of phenotypic diversity among mammals in Africa, we do not observe evidence that these changes were effected by decreases in population size (genetic drift). The strong association between genetic diversity and ecological region suggests that the exceptional diversity within the bushbuck complex may have been driven, at least in part, by parapatric speciation.