Integrative taxonomy using historical specimens provides evidence for a single species of bushbuck, Tragelaphus scriptus (Mammalia: Bovidae)

Abstract

The bushbuck (Tragelaphus scriptus) is an antelope widely distributed throughout sub-Saharan Africa that varies greatly in body size, horn length and coat patterning. Recent taxonomic work has suggested reclassification as two or more species. Molecular studies have shown a deep mitochondrial divergence between north-western (scriptus) and south-eastern (sylvaticus) bushbuck populations, and morphological assessments, based largely on coat patterning, have divided the bushbuck into as many as eight species and 32 subspecies. Molecular and morphological approaches have, to date, not been brought together in a single integrative study. Using 43 museum cranial specimens, we sequenced complete mitochondrial genomes and used three-dimensional geometric morphometrics to compare genetic and phenotypic divergence between the two main subclades of the T. scriptus species complex. We found a significant but weak relationship between phenotypic and genetic distances. Mean body size differed between males of each subclade, but morphological differences attributable to sexual dimorphism within subclades were greater than differences between subclades. Despite a mitochondrial divergence of nearly 3 Myr, cranial morphological differences between subclades were minor and could be accounted for largely by minor size differences between males (size–shape allometry). These results support the recognition of a single bushbuck species, as indicated by nuclear genetic data.

three-dimensional geometric morphometrics, cranial morphology, mitogenomics, museomics, species delimitation

INTRODUCTION

For millennia, morphological techniques have served as the primary methodological approach for testing taxonomic hypotheses. In only a few decades, increasingly inexpensive and accessible DNA sequencing techniques have enabled molecular systematics largely to supplant morphological analyses in taxonomic and phylogenetic research (at least on extant or recently extinct taxa). In particular, high-throughput sequencing and improvements in the successful isolation of ancient DNA material have created new opportunities for researchers working with molecular data from specimens in historical collections. However, the last decades have also seen the development of modern multivariate statistical methods in the analysis of morphology, with increases in computational power and the recent development and distribution of hardware and software for handling three-dimensional (3D) modelling and data. Three-dimensional geometric morphometric methods have reinstated quantitative morphological methods in modern zoological studies. Considering molecules and morphology, each of these methodological frameworks presents its own biases and limitations to the evolutionary biologist (Padial et al. 2010; Schlick-Steiner et al. 2010, Catalano et al. 2014). In order to overcome the difficulties of molecular or morphological research alone, researchers are now turning to integrative taxonomic methods.

Integrative taxonomy can provide systematics research with a more nuanced perspective of complex evolutionary histories (e.g. hybridization, the identification of cryptic species, reconciliation of gene discordance) than individual approaches alone (Chaplin et al. 2020). Technological and methodological advances have allowed new combinations of previously incompatible data types (e.g. shape data in geometric morphometrics with haplogroup data from mitogenomics; Gray et al. 2022). The production of 3D digital models and whole-genome sequencing (and the large quantities of data that these methods produce) have become easier to manage and analyse. Recent studies (e.g. Heckeberg 2020, Petzold and Hassanin 2020) have combined multiple lines of evidence successfully to investigate questions of species delimitation and patterns of evolutionary history.

In this study, we use integrative taxonomic methods to combine morphological and molecular data: 3D landmarks from surface scans and mitogenomes assembled via state-of-the-art molecular approaches, both from the same set of historical specimens. We address patterns of variation in different data types and recent evolutionary history of a widespread antelope species complex to resolve a long history of taxonomic conflict.

The bushbuck [Tragelaphus scriptus Pallas, 1766] is a medium-sized antelope that lives across most of sub-Saharan Africa and has adapted to a wide variety of habitats across its range. These adaptations have produced much variation in phenotypic characteristics, such as horn length and shape, body size, and coat patterning (Kingdon 1982). This phenotypic variation has motivated more than a century of morphological research on the taxonomy of the species complex. These studies have produced myriad classification systems suggesting as few as one species (Plumptre and Wronski 2013) and up to eight (Castelló 2016), or as many as 27 subspecies (Allen 1939).

Many of these morphological studies, regardless of the relative granularity of their taxonomic classifications, have observed a phenotypic split [defined either as separate species (Groves and Grubb 2011, Castelló 2016) or subspecies (Plumptre and Wronski 2013)] based on variations in coat patterning among populations in north-west and south-east Africa, with the presentation of two distinct phenotypes split roughly across the East African Rift Valley. Kingdon (1982) described both ‘harnessed’ (‘T. s. scriptus’) stock in the north-west and ‘sylvan’ (‘T. s. sylvaticus’) in the south-east, in addition to six other subtypes, but also noted the potential for much phenotypic variation even at a single location. Castelló (2016) elevates eight types of T. scriptus to species status: four are classified in the ‘scriptus’ group and the remainder in the ‘sylvaticus’ group. Drawing on a combination of data including coat patterning, linear cranial measurements and mitochondrial DNA (mtDNA) studies, Groves and Grubb (2011) describe eight species (three in a ‘T. scriptus cluster’; five in a ‘T. sylvaticus cluster’).

Interestingly, recent molecular studies of the bushbuck have described a similar north-west–south-east continental division between mitochondrial sequences. Indeed, two subclades of bushbuck haplogroups are divergent in mitochondrial phylogenies of the genus in that they sort as non-monophyletic clades, separated from one another by other tragelaphine species (Willows-Munro et al. 2005, Moodley et al. 2009). These mitochondrial subclades are also roughly divided by the Rift Valley into north-western and south-eastern populations (Moodley and Bruford 2007, Hassanin et al. 2018). In this study, we refer to the north-western mitochondrial subclade as the scriptus subclade and the south-eastern as the sylvaticus subclade (after Moodley and Bruford 2007).

Investigating this mitochondrial division, Hassanin et al. (2018) hypothesized that the mtDNA of the scriptus subclade was acquired by introgression after hybridization in the Early Pleistocene between T. scriptus and an extinct species closely related to Tragelaphus angasii Angas, 1849Angas 1848. More recently, a study combining mitochondrial and nuclear sequence data (Rakotoarivelo et al. 2019a) showed that this rift in the mitochondrial data is also present in the nuclear phylogeny, although to a much lesser extent. Notably, however, the nuclear data in that study support the monophyly of T. scriptus, despite the differing number of chromosomes between the two subclades [cited by Hassanin et al. (2018) as evidence for two species]. For Rakotoarivelo et al. (2019a), the hybridization event hypothesized by Hassanin et al. (2018) reconciles the discordance between mitochondrial and nuclear phylogenies and supports the hypothesis of a single species, with many zones of ongoing secondary contact between the two subclades.

Despite the apparent biogeographical correlation between phenotypic and genetic variation in the bushbuck, no study to date has integrated morphological and molecular data from the same specimens. Prior morphological taxonomic research on the bushbuck has been restricted to comparison of coat patterning and a few linear measurements (Groves and Grubb 2011); molecular studies of both mitochondrial and nuclear DNA have used only short segments. Here, we combine morphological and molecular data (3D cranial surface landmark data and whole mitochondrial genome sequences) drawn from the same set of specimens to determine whether patterns of variation in cranial morphology reflect the deep divergence observed in the mitochondrial subclades. We took advantage of improvements in ancient DNA extraction techniques and state-of-the-art molecular methods to enrich and sequence whole mitochondrial genomes from the same historical specimens from which we collected morphological data. Leveraging historical collections to study patterns of both genetic and phenotypic diversity using the same specimens, we highlight the importance of historical museum collections as repositories of biodiversity across space and time. This publication also serves as a case study for contemporary integrative taxonomic work, combining current morphological and molecular methods to clarify the complex evolutionary history and taxonomic status of the bushbuck. Using these integrative methods, we tested for correlation between phenotypic and genetic variation and determined that the two mitochondrial subclades of bushbuck cannot be distinguished by their cranial morphology.

MATERIALS AND METHODS

Geometric morphometrics

Specimens

We performed 3D geometric morphometric analyses on 43 adult bushbuck crania (27 males and 16 females) from the zoological collections of the Museum für Naturkunde, Berlin (Table 1). Thirty-nine of the specimens had location information (Fig. 1), broadly representing the modern range of the bushbuck. Twenty-two specimens had collection dates, which ranged from 1897 to 1938.

Table 1.

Bushbuck (Tragelaphus scriptus) cranial specimens used in the analyses. Specimens without GenBank accession numbers were not sequenced or did not yield sufficient mitochondrial DNA; NA subcade refers to these specimens. Subclade taxonomy was determined by molecular data (see Fig. 4). DRC refers to the Democratic Republic of the Congo.

ID	Specimen number	Sex	Country	Year	GenBank accession no.	Subclade
1	ZMB Mam 106581	Male	–	–	–	NA
2	ZMB Mam 106591	Male	Cameroon	–	ON474903	Scriptus
3	ZMB Mam 106592	Male	Angola	1903	–	NA
4	ZMB Mam 106597	Male	Cameroon	–	ON474904	Scriptus
5	ZMB Mam 106604	Male	Ghana	–	ON474905	Scriptus
6	ZMB Mam 106608	Male	Mozambique	–	–	NA
7	ZMB Mam 106609	Male	–	–	ON474906	Sylvaticus
8	ZMB Mam 106620	Male	Angola	1902	ON474907	Sylvaticus
9	ZMB Mam 106625	Male	Cameroon	1909	–	NA
10	ZMB Mam 106626	Male	South Africa	–	ON474908	Sylvaticus
11	ZMB Mam 106659	Male	–	1909	–	NA
12	ZMB Mam 106686	Male	Cameroon	1909	–	NA
13	ZMB Mam 106691	Male	Tanzania	1909	ON474909	Scriptus
14	ZMB Mam 107602	Female	Equatorial Guinea	1909	ON474910	Scriptus
15	ZMB Mam 107662	Female	Cameroon	1912	ON474911	Scriptus
16	ZMB Mam 107664	Female	DRC	–	ON474912	Sylvaticus
17	ZMB Mam 107665	Female	Tanzania	1897	ON474913	Sylvaticus
18	ZMB Mam 107667	Female	Ethiopia	1900	ON474914	Sylvaticus
19	ZMB Mam 107680	Female	Cameroon	–	ON474915	Scriptus
20	ZMB Mam 107683	Female	Angola	–	ON474916	Scriptus
21	ZMB Mam 107722	Female	Angola	–	ON474918	Scriptus
22	ZMB Mam 107732	Female	Tanzania	1909	ON474919	Sylvaticus
23	ZMB Mam 107735	Female	Cameroon	1910	ON474920	Scriptus
24	ZMB Mam 107740	Female	Cameroon	1904	ON474921	Scriptus
25	ZMB Mam 107751	Male	Cameroon	–	ON474922	Scriptus
26	ZMB Mam 107796	Female	South Sudan	–	–	NA
27	ZMB Mam 107812	Female	Cameroon	–	–	NA
28	ZMB Mam 16105	Female	Togo	1912	ON474890	Scriptus
29	ZMB Mam 2169	Male	South Africa	–	–	NA
30	ZMB Mam 31590	Male	Tanzania	–	–	NA
31	ZMB Mam 32456	Male	Ethiopia	–	ON474891	Sylvaticus
32	ZMB Mam 37177	Male	DRC	–	ON474892	Sylvaticus
33	ZMB Mam 37178	Male	DRC	–	ON474893	Sylvaticus
34	ZMB Mam 40389	Male	Angola	1924	ON474894	Scriptus
35	ZMB Mam 40543	Male	Ethiopia	–	ON474895	Sylvaticus
36	ZMB Mam 41189	Female	Senegal	1929	ON474896	Scriptus
37	ZMB Mam 41729	Male	Kenya	1925	ON474897	Sylvaticus
38	ZMB Mam 45680	Male	South Africa	1933	ON474898	Sylvaticus
39	ZMB Mam 48268	Male	Nigeria	1910	ON474899	Scriptus
40	ZMB Mam 4931	Male	–	–	ON474889	Scriptus
41	ZMB Mam 60992	Male	Tanzania	1938	ON474900	Sylvaticus
42	ZMB Mam 60995	Male	Tanzania	1938	ON474901	Sylvaticus
43	ZMB Mam 90802	Female	Tanzania	1913	ON474902	Sylvaticus

Figure 1.

Map of specimens with location information used in this study (39 individuals; four additional specimens did not have location information and are not shown here). Mitochondrial subclade designations (scriptus and sylvaticus) are based on the molecular phylogeny, and are designated NA (grey) where this could not be determined. Location coordinates are based on locality information on historical specimens and were estimated as closely as possible. The modern range of the bushbuck is shown in grey (IUCN SSC Antelope Specialist Group 2016). See Supporting Information, Fig. S2 for a map labelled with specimen numbers and countries.

Digitization, landmarking and principal component analysis

We created 3D surface models of all crania using a 3D surface scanner (Artec Spider) and Artec Studio Professional software (v.13). We scanned all external surfaces of the crania, then aligned and registered the partial scans, created a fast fusion to check for missing surfaces (repeating the scanning, alignment and registration steps when necessary), removed outliers and non-cranial elements, created a sharp fusion, applied texture, simplified the fusion to one million vertices, and exported the model in ply format with vertex colour. We used default parameters except where otherwise specified.

To characterize cranial shape, we used a set of 53 fixed 3D landmarks placed on the surface of each cranial model at anatomically homologous points chosen to capture morphological variation in bovid crania (Bibi and Tyler 2022). The Supporting Information (Table S1) details the landmark numbers and anatomical positions, which are illustrated in Figure 2. We performed all landmarking and geometric morphometric analyses in RStudio (v.1.4.1103; RStudio Team 2019) and the R statistical environment (v.4.0.2; R Core Team 2020). Functions refer to the geomorph package (v.4.0.0; Adams and Otárola-Castillo 2013) except where otherwise stated. We placed landmarks on the centred surface models using the digit.fixed() function, then Procrustes-transformed the raw landmark coordinates using both the gpagen() function, in order to extract centroid size, and the function bilat.symmetry(), to account for bilateral symmetry reflected in the left and right landmarks (Supporting Information, Table S1). Finally, we performed a principal component (PC) analysis with the function gm.prcomp() on the symmetric component of shape variation from this latter set of transformed landmarks.

Figure 2.

Surface model of Tragelaphus scriptus cranium in dorsal, lateral and ventral views, with points representing 53 digitized landmarks used to characterize cranial shape (22 blue landmarks mirrored on left and right sides; 9 orange medial landmarks). See the Supporting Information (Table S1) for anatomical descriptions of the locations of landmarks.

In order to visualize the effects on cranial shape attributable to size and PCs, we selected one specimen close to the mean shape (ZMB Mam 41189) and warped it to the mean shape calculated with the function warpRefMesh(). We then warped this mean shape mesh into models representing the extremes of centroid size in our sample set, in addition to the minima and maxima of the first three PCs (Supporting Information, Fig. S1).

We used ordinary least squares ANOVA with randomized permutation to examine the relationships between cranial shape and three variables: mitochondrial subclade (scriptus vs. sylvaticus), sex and size. We tested the comparison shape ~ subclade for significant morphological differences between the two bushbuck subclades. Specimens were assigned to a subclade by the mitochondrial phylogeny (see Results). We used the comparison shape ~ sex to examine the significance of sexual dimorphism and compare it to differences between subclades (sexual determination in bushbuck crania is relatively straightforward, because females are hornless). Finally, we used the comparison shape ~ size to assess the significance of size allometry in the generation of morphological differences between subclades and to determine whether the two subclades follow similar allometric trajectories (slopes). Here, we used the natural log of centroid size as a proxy for body mass. For models without clade as a variable, we included all T. scriptus specimens (N = 42); otherwise, only specimens for which we obtained mitochondrial sequence information were included (N = 33). All comparisons used the procD.lm() function in geomorph (Adams and Otárola-Castillo 2013).

We visualized the relationship between phylogeny and morphospace using the phylomorphospace() function in the phytools package (v.0.7-70; Revell 2012), including only those bushbuck specimens for which we obtained both morphological and sequence data (N = 33). We used the function physignal() (multivariate version of the K-statistic; Blomberg et al. 2003, Adams 2014) to estimate the degree of phylogenetic signal present in the Procrustes-transformed shape coordinates of specimens for which we had sequence data.

Mitogenome reconstruction

Sampling

Cardini (2020) discussed the impact of autocorrelation in museum collections (specimens often come from the same locality and can often be closely related individuals), which can make them poor representatives of within-species variation. To counteract this, we have included specimens from many sampling locations across the species range, and the results of the phylogenetic analysis indicate that very few individuals in this study have an extremely close relationship (e.g. the individuals from Cameroon in Figs 4, 7).

Figure 3.

Plots of size allometry examined via ln-transformed centroid size and principal component 1 (PC1) (A); or regression of shape (B). Note the strong influence of size on shape. The key is the same as that shown in Figure 1; labels refer to specimens in Table 1.

Figure 4.

Maximum likelihood phylogeny of all sequenced bushbuck specimens with other tragelaphines. Node numbers indicate approximate likelihood ratio test (aLRT) support values; only nodes with > 95% support are shown. Nodes with < 80% support were collapsed. Numbers refer to specimens in Table 1 or GenBank reference specimens in Table 2. Dots are coloured by country and correspond to those in Figure 7.

Laboratory procedures

To determine the most productive, least destructive skeletal element source of DNA for this set of specimens, we first prepared ~30 mg of ground bone material sampled from five different locations on the cranium of one specimen (ZMB Mam 106608), then extracted the DNA using a QIAamp DNA Investigator Kit following the protocol described by Rohland and Hofreiter (2007) and modified by Dabney et al. (2013). After extraction, we determined the degree of DNA degradation using an Agilent D1000 DNA ScreenTape assay and TapeStation Analysis Software v.3.1.1. We measured DNA concentration with a Qubit 2.0 Fluorometer and dsDNA HS Assay Kit (Thermo Fisher Scientific). Only the sample from the nasal turbinate bone produced any detectable concentration of DNA (1.06 ng/µL). For the remaining specimens, we prepared 30–400 mg tissue samples from the nasal turbinate bones and tissue, extracted the DNA, and checked the degradation and concentration following the protocol outlined above (Supporting Information, Table S2).

We prepared multiplexed DNA libraries for paired-end Illumina sequencing using a NEXTFLEX Rapid DNA-Seq Kit 2.0 and Bundle with NEXTflex-HT Barcodes 1–96 (Bioo Scientific/PerkinElmer). We followed the manufacturer’s instructions for ‘Option1 – no size selection’, with the following exceptions: DNA was not sheared before library preparation; Clean-Up Bead Volume was set to 45 µL instead of 35 µL after Adapter Ligation, to avoid losing small target-fragments; and Resuspension Buffer was decreased to 25 µL to concentrate amplified libraries. Finally, library size, quality and quantity were determined as above. To obtain a sufficient quantity for subsequent enrichment, several libraries were re-amplified.

After this, we enriched the libraries for targeted sequencing using myBaits (Arbor Biosciences, USA). We designed two different biotinylated RNA baits, complementary to sequence targets for sylvaticus and scriptus subclades (Kit A and B), using existing mitogenomes from GenBank (Table 2). Assuming a high amount of non-target templates, we followed the recommendations for ‘Enriching low-copy, degraded and/or contaminated DNA Libraries’, using the Arbor Bioscience myBaits protocol for ‘Dual Round Enrichment’ and the myBaits manual v.4. In total, we enriched 44 individuals in pools of two to four libraries per enrichment reaction in the following groups: 16 individuals/libraries with Kit A (sylvaticus); 17 individuals/libraries with Kit B (scriptus); 6 ‘Unknown’ individuals/libraries, each sample with both kits; and 3 ‘Unknown’ individuals/libraries with Kit A (sylvaticus). Enrichment pools were checked for quality and quantity as above.

Table 2.

Ten additional Tragelaphus mitogenomic sequences (nine from GenBank and one from this study) included in the molecular analyses. These include two reference sequences (Tragelaphus scriptus) used in the bioinformatics workflow and eight sequences representing the remaining tragelaphine species.

Taxon	GenBank ID/ Specimen No.	Source
T. scriptus	JN632707.1 (sylvaticus reference)	Hassanin et al. (2012)
T. scriptus	JN632705.1 (scriptus reference)	Hassanin et al. (2012)
T. buxtoni	NC_038064.1	Li et al. (2018)
T. oryx	NC_020750.1	Hassanin et al. (2012)
T. strepsiceros	NC_020752.1	Hassanin et al. (2012)
T. eurycerus	NC_020749.1	Hassanin et al. (2012)
T. spekii	NC_020620.1	Hassanin et al. (2012)
T. angasii	NC_020748.1	Hassanin et al. (2012)
T. imberbis	NC_020619.1	Hassanin et al. (2012)
T. spekii	ON474917 (ZMB Mam 107720)	This study

Sequencing was carried out at the Berlin Center for Genomics in Biodiversity Research (BeGenDiv/DE-Berlin). After test sequencing on an Illumina MiSeq, sequences were run in two batches on an Illumina NextSeq using paired-end reactions, the first with 26 enriched libraries and the second with 22 enriched libraries and one control, representing a total of 34 specimens: 33 bushbucks (T. scriptus) and one sitatunga (Tragelaphus spekii Speke, 1863).

Read cleaning and mitochondrial reconstruction

After sequencing, we cleaned the raw reads using a custom-designed workflow for museum specimens (see https://github.com/MozesBlom/NGSdata_tools). This workflow (see Supporting Information, Fig. S4) consists of: (i) removal of duplicates using Super deduper (v.1.0.0; Petersen et al. 2015); (ii) adapter trimming using Trimmomatic (v.0.36; Bolger et al. 2014); (iii) merging of overlapping read pairs with PEAR (v.0.9.10; Zhang et al. 2013); (iv) quality trimming with Trimmomatic (v.0.36; Bolger et al. 2014); and (v) removal of low-complexity reads using a Python script developed in house. Before and after cleaning, we inspected the quality of the reads visually using FastQC (v.0.11.5; Andrews 2010).

Previous studies have reported a deep mitochondrial split between T. scriptus and T. sylvaticus Sparrman, 1780 (see Introduction), and reference selection can have an impact on downstream evolutionary analyses (Prasad et al. 2022). To remove potential biases arising from unequal phylogenetic distance between focal specimens and the reference, we mapped each individual twice and used two reference genomes from GenBank (Table 2); one belonging to each of the two mitochondrial subclades (JN632705.1 for the scriptus subclade; JN632707.1 for sylvaticus). By mapping each individual against both references in parallel and comparing the obtained datasets, we were able to identify differing sites and thus mitigate reference bias.

We mapped the polished reads to each reference using the BWA-MEM algorithm (Li 2013) and used SAMtools (v.0.1.19; Li et al. 2009) as follows: (i) to convert the resultant SAM files into BAM format; (ii) to remove reads that did not pair properly with the reference; (iii) to sort and index the BAM files; and (iv) to merge the library-specific BAM files by individual. Next, we used FreeBayes (v.1.0.2; Garrison and Marth 2012) to call variants while considering only bases with a phred quality score > 20 from reads with a mapping quality score > 20. Subsequently, we filtered the variants with VCFtools (v.0.1.15; Danecek et al. 2011) and vcflib (Garrison et al. 2022) as follows: (i) by eliminating individual genotypes with a depth less than eight reads; (ii) by decomposing complex variants; and (iii) by removing indels. Finally, we created a consensus sequence for each individual and each reference by using BCFtools (v.1.3.1; Danecek et al. 2021). During variant calling, we ran FreeBayes deliberately under a fixed ploidy of two, requiring at least two reads supporting the alternative allele. However, mtDNA is haploid, and consistent support for alternative alleles could be the result of contamination. Therefore, we incorporated the reference allele in heterozygous sites only when calling the draft consensus sequences using BCFtools.

In order to reduce the risk of incorporating non-endogenous sequence data from other museum specimens, we carried out an additional, unconventional filtering procedure. This step consisted of a second mapping and variant calling iteration. We repeated the mapping and variant calling pipeline using the same settings as described above but used the consensus sequences called for each of the two original reference genomes instead. Ideally, reads mapped against their own consensus sequence should align perfectly, and spurious genetic variants could indicate DNA contamination or mapping errors. Therefore, we masked all positions at which alternative alleles persisted. After this second iteration, we also masked sites with a read depth below eight to avoid the incorporation of reference alleles in low-coverage regions. Finally, we compared the two consensus sequences for each individual, masked sites where they differed and used the resulting individual consensus sequence for downstream phylogenetic analyses.

In addition to the 33 bushbuck specimens successfully sequenced in this study, the two reference sequences, and the T. spekii sequence (ZMB Mam 107720), we also included in our phylogenetic analyses eight complete mitogenome sequences (from GenBank) representing all other tragelaphine species (Table 2). We aligned these sequences using the FFT-NS-2, progressive method in MAFFT (v.7.453; Katoh et al. 2002, Katoh and Frith 2012) and checked the alignment by eye in Mesquite (v.3.61; Maddison and Maddison 2019).

We then determined the best-fitting site model of gene evolution for this alignment (N = 49) using Bayesian information criterion in ModelFinder (Kalyaanamoorthy et al. 2017). We created a maximum likelihood (ML) tree using the desktop version of IQ-TREE (v.1.6.12; Nguyen et al. 2015) and calculated branch support using the Shimodaira–Hasegawa approximate likelihood ratio test (SH-aLRT; Guindon et al. 2010). We also produced a Bayesian tree in BEAST (v.2.6.7; Bouckaert et al. 2019; see Supporting Information, Fig. S3), which corroborated the results of the ML tree (Supporting Information, Fig. S4). To visualize the relationships among haplotypes, we used MEGA11 (v.11.0.8; Tamura et al. 2021) to remove all sites with ambiguous bases or gaps in one or more of the sequences in the alignment, then constructed a median-joining network (Bandelt et al. 1999) in PopART (v.1.7; Leigh and Bryant 2015). Given that most sequences were of similar length, sequences with many gaps were removed from the alignment before removing these sites and constructing the median-joining network. To calculate phylogenetic signal in the morphological data, we used the function physignal() to calculate a multivariate K-statistic and performed a Mantel test between the molecular and morphological distance matrices.

We used the following R packages for data manipulation and visualization: abind (v.1.4-5; Plate and Heiberger 2016), ade4 (Dray et al. 2007), ggpubr (v.0.4.0; Kassambara 2017), ggrepel (v.0.9.1; Slowikowski 2018), ips (v.0.0.11; Heibl et al. 2019), RColorBrewer (v.1.1-2; Neuwirth 2014), scales (v.1.1.1; Wickham and Seidel 2020), sf (v.0.9-8; Pebesma 2018) and tidyverse (v.1.3.1; Wickham et al. 2019).

RESULTS

Geometric morphometrics

Subclade, sex and size all have a significant effect on shape, although size and, especially, sex explain higher proportions of shape than does subclade (Table 3). After accounting for sex and size (shape ~ sex + size + subclade), differences in shape between scriptus and sylvaticus subclades are significant but small (Table 3). Although size and subclade both have a significant effect on shape, the interaction between size and subclade (shape ~ size * subclade) is non-significant (N = 33, P = 0.556, Z = −0.1414) in explaining shape, indicating similar size allometry in both scriptus and sylvaticus subclades. The interaction between size and sex in explaining shape (shape ~ size * sex) is also non-significant (P = 0.062, Z = 1.6077), indicating similar size allometry between males and females (despite the sexually dimorphic presence/absence of horns). The relationship between ln-transformed centroid size and the regression of size on shape is shown in Figure 3; a pairwise test of slopes confirms no significant differences between the two subclades (P = 0.0978).

Table 3.

Results of the ANOVA of sex, size and subclade on shape. Models that include the subclade variable were calculated using the subset of bushbuck specimens for which we obtained sequence data (N = 33); all other models included all bushbuck specimens (N = 43). The table includes significant results only; see the Results section for discussion of non-significant results.

Shape ~	Residual d.f.	d.f.	RSS	SS	MS	Rsq	F	Z	Pr(>F)
Size	41	1	0.096	0.015	0.015	0.139	6.593	3.726	0.001
Sex	41	1	0.092	0.020	0.020	0.177	8.843	4.586	0.001
Size + sex	40		0.088	0.024	0.015	0.214	5.437	4.555	0.001
Size		1		0.015	0.015	0.139	7.047	3.793	0.001
Sex		1		0.008	0.008	0.075	3.826	3.977	0.001
Subclade	31	1	0.076	0.009	0.009	0.107	3.710	2.861	0.002
Size + sex + subclade	29		0.058	0.027	0.009	0.318	4.508	4.219	0.001
Size		1		0.013	0.013	0.149	6.346	3.625	0.001
Sex		1		0.008	0.008	0.090	3.842	3.488	0.001
Subclade		1		0.007	0.007	0.078	3.335	3.265	0.001
Size + subclade	30		0.065	0.020	0.010	0.234	4.580	3.840	0.001
Size		1		0.013	0.013	0.149	5.844	3.542	0.001
Subclade		1		0.007	0.007	0.085	3.316	3.035	0.001

Abbreviations: MS, mean squares; Pr(>F) = P-value for F-value; Rsq, coefficient of shape variation attributed to each covariate; RSS, residual sum of squares; SS, sum of squares; Z, effect size (based on F-values).

Figure 5.

A, plot of all specimens along the first two principal components (PCs) of shape. B, cranial outlines transformed via thin-plate spline from the mean shape to the extremes of PCs 1 and 2. Note the separation of subclades along PC1, which largely reflects sexual dimorphism and size allometry (see Fig. 3). The key is the same as that shown in Figure 1; labels refer to specimens in Table 1.

Differences in size between scriptus and sylvaticus subclades are non-significant (P = 0.243, Z = 0.7417). Considering males and females separately, differences in size (size ~ subclade) are also not significant (males: N = 19, P = 0.202, Z = 0.8839; and females: N = 14, P = 0.307, Z = 0.5835). Interestingly, however, the relationship between shape and subclade (shape ~ subclade) is small but significant (Table 3). When sexes are analysed separately, this is significant only for males (P = 0.004, Z = 2.9857), not for females (P = 0.06, Z = 1.6139).

Principal component analysis

In Figure 5, we plot all specimens along the first two PCs of shape, which describe 24% and 10% of the variation, respectively. The remaining 41 components each describe < 9% of the variation.

Although the morphological differences between all extremes of shape and size are subtle, the largest skull (Centroid Size Max) resembles the minimum PC1 skull, and the smallest skull (Centroid Size Min) resembles the maximum PC1 skull (Fig. 3A). This correlation between PC1 and size is unsurprising, because the first PC of geometric morphometric analyses often loosely reflects size–shape allometry.

Principal component 1 describes 24.4% of the total shape variation and relates to both sexual dimorphism and size, with differences in cranial shape relating to the length and robustness of the cranium. The differences in shape variation described by the PC1 are similar to those of the models of minimum and maximum centroid size (Supporting Information, Fig. S1). Principal component 2 describes 10.2% of the total variation in shape and is related to the length of the snout and roundness of the braincase. Principal component 3 describes 8.1% of the total variation in shape and is related to the medial curvature of the ventral skull.

Phylogenetic inference

Our molecular phylogeny agrees with the deep divergence between scriptus and sylvaticus subclades found by other mtDNA studies (Hassanin et al. 2018). We also find a high level of structure in the sylvaticus subclade compared with scriptus, which might be attributable to the smaller sample size of scriptus individuals in our study or might be because of collection bias (Cardini 2020). However, this could also represent a real pattern of higher genetic structure and diversity in the sylvaticus group (additional data on the genetic structure of the sample and the results of the molecular experiments are given in the Supporting Information, Supplemental material; Table S3).

Phylogenetic signal

Sylvaticus and scriptus clades can be distinguished readily along PC1 (primarily representing size allometry), but within each subclade there is no clearly discernible relationship between phylogeny and morphology (Fig. 6). There is no significant evidence for phylogenetic signal in the Procrustes-transformed shape data (K ≈ 0, P = 0.6355, Z = −0.3522). However, we do find a small but significant positive correlation between molecular and morphological distance (r = 0.226, P = 0.0013; for additional distance comparisons subdivided by sex and clade, see Supporting Information, Fig. S5; Supplemental material).

Figure 6.

Phylomorphospace shows poor correspondence between mitogenomic relationships and morphological differences. Axes refer to the first and second principal components (PC1 and PC2) of the geometric morphometric analysis. The key is the same as that shown in Figure 1; labels refer to specimens in Table 1. Note that this figure includes only specimens for which mitochondrial DNA was extracted.

DISCUSSION

Morphology

Although we find significant differences in shape between the scriptus and sylvaticus subclades, these are only between males and not females. The observed differences can also be accounted for by minor (non-significant) size differences, particularly given that males of the sylvaticus subclade are larger, on average, and attain a larger maximum size (with high overlap). The amount of cranial shape variation explained by differences between subclades is far less than that explained by sexual dimorphism or size differences (Table 3), indicating minimal morphological divergence. Furthermore, although phylogenetic relationships differ along PC1, which largely reflects size–shape allometry (Figs 3, 5, 6), and although molecular and morphological distances are positively correlated, the effect is small, and there is no phylogenetic signal in total cranial shape.

Changes with increasing size (and more characteristic of sylvaticus) include a deeper face, a taller skull, and higher and more laterally set orbits (Supporting Information, Fig. S1). However, whether such differences could be detected consistently by eye at the individual level is not at all clear. The small magnitude of these differences and their high range of overlap indicate that the two bushbuck subclades do not exhibit strong morphological differentiation, at least as far their cranial shape is concerned. Given that the specimens used in this study represent nearly the full continental range of the species, this low amount of difference in cranial shape between the two subclades, which are separated both geographically and temporally and which exhibit remarkable biogeographical variation in coat patterning (Kingdon 1982), is striking. Despite a deep mitochondrial divergence, our geometric morphometric analyses found no significant evidence for morphological divergence between the scriptus and sylvaticus subclades.

Mitogenomics

Our mitochondrial phylogeny and haplotype network (Figs 4, 7) agree with previous mitochondrial phylogenies of bushbuck, confirming the presence of two deeply divergent and geographically distinct subclades within the species.

Figure 7.

Median-joining network of 28 sequences from bushbuck specimens, with missing and ambiguous sites removed. [Three sequences with many gaps (ZMB Mam 32456, ZMB Mam 40389 and ZMB Mam 107683) were removed before the analysis.] Branches are labelled with mutation distances; unlabelled branches have a distance of one. Note the structure, the long branch separating clades and the Tanzania/Nigeria haplogroup. One reference sequence was included from both clades (Table 2).

Using mitochondrial fragment data, Hassanin et al. (2018) dated the most recent common ancestor of sylvaticus to 2.61 [95% CI: 2.1–3.1] Mya and the most recent common ancestor of scriptus to 2.31 [1.7–3.0] Mya (largely in agreement with Moodley and Bruford 2007). Their tree based on the concatenation of nuclear genes corroborated this older origin of the sylvaticus clade: 0.78 [1.1–0.5] Mya for the most recent common ancestor of sylvaticus and 0.39 [0.6–0.2] Mya for that of scriptus. Rakotoarivelo et al. (2019b) reported similar findings and hypothesized that the bushbuck probably originated in East Africa, with scriptus dispersing north-west and sylvaticus dispersing south-east via a series of vicariance events. They described 23 phylogenetically distinct haplogroups (eight in scriptus; 15 in sylvaticus), showing that much of this structural variation aligns with habitat. This older age of origination in sylvaticus and the diversity of ecoregions in south-eastern Africa (owing, in part, to the greater span of latitude that sylvaticus inhabits) explain the higher level of phylogenetic structure in the sylvaticus subclade, something also observed in our findings (Fig. 8).

Figure 8.

Distribution of bushbuck crania according to size (log centroid size; A) and along the first three principal components of shape (B–D). Differences in size are non-significant between subclades and when assessed separately for males and females, despite the attainment of a larger maximum size in sylvaticus males. Differences in shape between subclades are significant, but only for males, not females.

Additionally, the advantages of our approach, which targeted morphology and DNA sequences from the same museum individuals, is highlighted by our finding of at least three specimens with mismatching geography and subclade. Based on their provenance, ZMB Mam 106620 from Angola (34 and 21 in Fig. 1) and ZMB Mam 106691 from central Tanzania (41 in Fig. 1) would have been assumed to belong to the sylvaticus subclade but appear in the scriptus group in Figure 4. However, their mitochondrial sequences place them securely within the scriptus clade. It is most likely that the location information associated with these specimens is misleading. Although it is possible that this incongruence indicates that dispersal and gene flow across the sylvaticus–scriptus divide can cross substantial distances, it is extremely unlikely that this is the case for specimen ZMB Mam 106691 (the locality given for this specimen is Dodoma, the capital of Tanzania, which might merely indicate that the specimen was acquired there, perhaps at a market). Interestingly, female bushbucks often spend their lives near their birthplace, whereas young male bushbucks tend to leave their home range at maturity (Kingdon 1982), which could promote the rapid formation of phylogeographical structure at the mitogenomic level, but not at the nuclear level.

CONCLUSION

Until recently, the incorporation of phylogenomic and morphological data had been restricted computationally and methodologically. Here, we took advantage of advances in integrative methods and incorporated 3D shape data with mitochondrial phylogenetics to investigate patterns of variation in mtDNA and cranial morphology in the bushbuck species complex. The lack of phylogenetic signal in cranial morphology suggests that the bushbuck should remain classified as a single species, in agreement with nuclear DNA analyses (Hassanin et al. 2018, Rakotoarivelo et al. 2019a, b). The success of this approach in a species delimitation problem supports the use and development of integrative methods in taxonomic studies and the further incorporation of morphological, phylogenetic, biogeographical and ecological datasets.

SUPPLEMENTARY DATA

Supplementary data is available at Zoological Journal of the Linnean Society online.

Table S1. Identification numbers and anatomical descriptions of the 53 landmarks used in the geometric morphometric analyses.

Table S2. Results of lab enrichment procedures.

Table S3. Coverage and DNA read statistics.

Figure S1. Models representing the extremes of centroid size in our sample set, as well as the minima and maxima of the first three principal components.

Figure S2. Map of specimens with specimen number and country of origin.

Figure S3. Bayesian phylogeny of all sequenced bushbuck specimens with other tragelaphines.

Figure S4. Visualisation of pipeline used for read cleaning.

Figure S5. Pairwise comparisons between molecular and morphological distance for all specimens with sequence data.

CONFLICT OF INTEREST

None declared.

ACKNOWLEDGEMENTS

We are grateful for the editorial staff of the Journal, who worked hard to facilitate the publication process during a time of transition; two anonymous referees for their time, suggestions, and comments, which contributed greatly to the development of the manuscript; and for support provided by the Abgeordnetenhaus of Berlin and the Museum für Naturkunde's Integrative Taxonomy funds. We also thank E. Amson, I. Baird, D. Breyer, C. Funk, N. Heckeberg, E. Hempel, I. Lazagabaster, K. Mahlow, F. Mayer, S. Mbedi, A. Petzold, T. Ramm, K. Spitzer and D. Willborn.

AUTHOR CONTRIBUTIONS

C.N.B. and F.B. designed the study and secured funding. C.N.B. scanned and landmarked all crania, with input from F.B. I.W. and C.N.B. isolated the DNA, and I.W. developed methods for DNA extraction and enrichment and prepared the libraries for sequencing. M.E. and M.P.K.B. conducted bioinformatics work and generated the sequence alignments. M.E., I.W., M.P.K.B. and F.B. provided support with the analyses and methodology. C.N.B. conducted all morphological and phylogenetic analyses with input from F.B. C.N.B. led the writing of the manuscript, and all authors contributed to revision of the manuscript.

FUNDING

This work was supported by the Taxonomy Fund of the Museum für Naturkunde; and the Studienstiftung des Abgeordnetenhauses von Berlin to C.N.B. The Berlin Center for Genomics in Biodiversity Research team was partially funded by the German Federal Ministry of Education and Research (BMBF, Förderkennzeichen 033W034A).

DATA AVAILABILITY

DNA sequences are available through GenBank. Three-dimensional surface models and an R script to reproduce all results are provided as Supplementary materials (https://doi.naturkundemuseum.berlin/data/10.7479/abew-e365).