Materials and Methods

Molecular data: Nine nuclear genes were selected based on the availability of sequence data for a broad selection of bovids: spectrin beta nonerythrocytic 1 (SPTBN1, table S1.1), kappa casein exon 4 (CSN3, table S1.2), protein kinase C iota (PRKCI, table S1.3), amelogenin X-linked (AMELX, table S1.4), thyroid stimulating hormone beta (TSHB, table S1.5), stem cell factor (SCF, table S1.6), signal transducer and activator of transcription 5A (STAT5A, table S1.7), thyroglobulin (TG, table S1.8), and tumor necrosis factor alpha (TNF, table S1.9; all supplemental tables available in the online supplement,  https://doi.org/10.5531/sd.sp.49). There were 107 bovid species with sequence data available for at least one of these genes on GenBank. I also obtained sequences for at least one of these genes for domestic pigs (Sus scrofa), Reeves' muntjac (Muntiacus reevesi), red deer (Cervus elaphus), mule deer (Odocoileus hemionus), caribou (Rangifer tarandus), and roe deer (Capreolus capreolus) to serve as outgroups. When a species had multiple sequences for a single gene available on GenBank, all sequences were aligned using the local alignment MAFFT algorithm, L-INS-I (Katoh et al., 2002). The majority consensus sequence for that species was then computed, ignoring gaps introduced in the alignment, in Geneious 10.2.3 ( https://geneious.com). BLAST searches of consensus sequences were performed against the GenBank nonredundant nucleotide database (nt) to ensure that they still produced close matches to the original genes. The gene sequence for each species was aligned using the L-INS-I algorithm in MAFFT, and after separate alignment of each gene, the alignments were then combined into a single nuclear data matrix with 7,372 molecular characters.

Mitochondrial genomes are available for many ruminant species, and thus they were an important source of molecular characters for the bovids in this study. Complete mitochondrial genomes for 123 species, including all six outgroup species from the nuclear DNA matrix (table S1.10 in the online supplement,  https://doi.org/10.5531/sd.sp.49) were retrieved. Ancient DNA mitochondrial genomes were included for the helmeted muskox (Bootherium bombifrons), the aurochs (Bos primigenius), and the steppe bison (Bison priscus). Because mitochondrial genomes are circular, but MAFFT can treat sequences only as linear during alignment, nine species' sequences were edited based on their GenBank annotations to start with tRNA-Phe, the starting point for all other sequences: Bison bison, B. priscus, Bos grunniens, Bos indicus, Bos mutus, Bos javanicus, Bos primigenius, Bos taurus, and Syncerus caffer. The global alignment algorithm of MAFFT, G-INS-I, was then applied to produce a mitochondrial genome alignment with 18,069 characters.

Morphological data: The morphological data comprised 196 craniodental and postcranial characters coded for 88 taxa analyzed across four studies (Gentry, 1992; Thomas, 1994; Bärmann, 2013, 2014). These data include 13 fossil species from nine of the tribes of Bovidae (table S2 in the online supplement,  https://doi.org/10.5531/sd.sp.49). Gentry (1992) analyzed fossil taxa separately from the extant matrix (with the exception of three extant species that appear in both matrices), and had several characters coded only for fossils; for this study this fossil matrix was combined with the extant matrix based on the character correspondences provided in the paper (Gentry, 1992). No character descriptions were available for characters 33–36 in the fossil matrix, thus they were not included in this analysis. Thomas (1994) included essentially the same set of characters from Gentry's (1992) matrix, with the addition of two characters. He also added the saola (Pseudoryx nghetinhensis), anoa (Bubalus depressicornis), and Sumatran serow (Capricornis sumatraensis), as well as two nonbovid outgoups, the Indian muntjac (Muntiacus muntjak) and the Siberian musk deer (Moschus moschiferus). Musk deer characters from Thomas' (1994) matrix were not included here, as the species was not represented in any other data partition in this study. Additional craniodental characters were incorporated from Bärmann's studies of bovids (2013; Morphobank Project 352) and the tribe Antilopini (2014). Behavioral and soft tissue characters from these analyses were also not included in the present study, because they are features not commonly preserved in the fossil record and thus would be unlikely candidates to provide support for the placement of fossil taxa in phylogenies. When character descriptions in the different matrices unambiguously referred to the same feature, the characters were combined, reordering codings or combining relevant states as needed. Character descriptions from each of these sources are included in appendix 1 with explanations of how states were modified to facilitate the matrix combination where appropriate. Thirty-four of these morphological characters were also encoded for four cervids based on Marcot's (2004) ruminant matrix to serve as an outgroup comparison for morphological character optimizations. Muntiacus reevesi and Odocoileus hemionus, two cervids for which much of the molecular data was available, were encoded as chimeras based on the congenerics Muntiacus muntjak (see Thomas, 1994; Marcot, 2004) and Odocoileus virginianus (see Marcot, 2004). A table describing which species are represented in each data partition is located in the online supplemental information 2, and all character encodings are available in the phylip matrix file (supplemental file S3) and table S4.6 (available online,  https://doi.org/10.5531/sd.sp.49).

Maximum likelihood tree search: To assess how different data partitions might support the final total evidence analysis, separate analyses of each matrix (morphological, nuclear DNA, mitochondrial DNA), as well as a combined nuclear and mitochondrial DNA matrix, were conducted. Each molecular dataset was partitioned into coding and noncoding regions, using by-codon partitions for each protein-coding gene. The best fitting model of molecular evolution available in RAxML, GTR with gamma rate heterogeneity (GTR+G) or GTR+G with a proportion of invariant sites (GTR+I+G), was then assessed using the greedy search algorithm and linked branch lengths in PartitionFinder (Lanfear et al., 2017) on the CIPRES science gateway. GTR+G was the best-fitting model for the nuclear data analysis, and GTR+I+G was the best-fitting model for the mtDNA, combined molecular, and total evidence analyses. Before running each tree inference in RAxML, the best-scoring rearrangement setting was assessed following the recommendations of the RAxML v.8 user manual (Stamatakis, 2014, 2016). To find the topology with the best likelihood score, 200 inferences were then run for each separate partition. One-thousand-replicate nonparametric bootstrap branch supports were computed and applied to the topology with the best likelihood score (Stamatakis, 2014, 2016). For the morphological matrix, the multistate gamma MK model of evolution with an ascertainment bias correction (Lewis, 2001) was used, and the same procedures for determining a rearrangement setting, obtaining the best topology, and calculating bootstrap support values were followed. RAxML does not accept polymorphic character states, and so any polymorphic characters in the morphological matrix were converted to missing before analysis.

Next, both morphological and molecular data were combined into a single total evidence matrix to estimate the total evidence phylogeny of Bovidae. This final matrix included 156 species. Analysis of the best rearrangement setting, identification of the most likely topology, and calculation of bootstrap support values were performed following the RAxML manual (Stamatakis, 2014, 2016) as was done for the analyses of the separate partitions. The same partitioning scheme for molecular data and model of molecular evolution from the analysis of the combined nuclear and mitochondrial matrix (GTR+I+G) was used with the addition of the MK model with ascertainment bias correction (Lewis, 2001) for the morphological characters. Bootstrap support values were calculated and applied to the most likely topology following the same procedures as in all other analyses. The final total evidence matrix in phylip format can be found in file S3 (available at  https://doi.org/10.5531/sd.sp.49). To assess the effects of missing morphological data on the ability to estimate fossil relationships, a total evidence phylogeny was also computed for the 88 species represented in the morphological matrix using the same models of molecular and morphological evolution as the full total evidence analysis.

Bayesian tree search: A time-calibrated Bayesian phylogeny was computed to examine potential differences in topology introduced by the choice of phylogenetic estimation method. Rather than combine partitions under the scheme assessed for the maximum likelihood analysis, all partitions were allowed to have separately-estimated models of molecular evolution, averaging across all reversible models available in the bModelTest package for BEAST2; all default priors were used for these models (Bouckaert and Drummond, 2017). When a single model is not strongly favored over others, averaging over all models rather than selecting a single model for each partition can account for the uncertainty (Bouckaert and Drummond, 2017).

Age ranges from Paleobiology Database (paleobiodb.org) provided the calibration times for all fossil species except Sivatragus bohlini, Tragelaphus sp., and Turcoceros grangeri (table 1). S. bohlini is Pliocene in age, according to Gentry (1992), thus I used the Pliocene stratigraphic age range (5.333–2.58 million years) from the International Chronostratigraphic Chart v 2020/03 (Cohen et al., 2013). For Tragelaphus sp., I used an age of 5.72 million years, based on Bibi (2013), and for T. grangeri, from the Tunggur Formation (Geraads, 2003), I used a range of 13–11.8 million years based on the age range of the formation as estimated from paleomagnetic correlations (Xiaoming et al., 2003). I used the fossilized birth-death process to model speciation, estimating all parameters except ρ (rho), the probability that an extant species would be sampled, which I fixed to 1 (Heath et al., 2014; Barido-Sottani et al., 2018). As the starting age for estimating the origin of Bovidae, I used a uniform distribution of 58.6 ±3.2 million years, based on the uniform dating model age with hard fossil calibration points estimated in a large-scale mitochondrial phylogeny of Artiodactyla (Hassanin et al., 2012).

Ancestral character estimation: Finally, to determine which characters map as synapomorphies for Bovidae or for clades within Bovidae, marginal ancestral state reconstructions were performed for each morphological character using RAxML-NG (Kozlov et al., 2019). Character 118 was excluded from the ancestral state reconstructions because it had either the same state encoded for all taxa or was polymorphic and marked as missing for analysis in RAxML, and thus could not be optimized in RAxML-NG. State changes were estimated on the most likely total evidence topology with the multistate MK model and correcting for ascertainment bias as in the morphological partition tree search. I then examined the best-supported character states from these analyses for state transitions at the nodes for Bovidae, the subfamilies Bovinae and Antilopinae, and the tribes as identified previous large-scale artiodactyl phylogeny using mitochondrial DNA (Hassanin et al., 2012). Because the total evidence analysis placed some fossil taxa in clades outside the tribes to which they most likely belong, I performed a second ancestral state reconstruction with these misplaced fossil species pruned from the phylogeny. I compared both ancestral state reconstructions to the characters coded for the fossil species to assess whether the synapomorphies would support their inclusion in those tribes, even if they were not correctly placed during the phylogenetic analysis.

TABLE 1.

Bayesian analysis fossil calibrations.

Topology

The most likely topology from the maximum likelihood total evidence analysis recovered monophyletic tribes of extant bovids; however, support for most tribes and internal nodes was below 75% (fig. 1). Hippotragini, Reduncini, Antilopini, Aepycerotini, and Bovini all had bootstrap support well below 50%. Caprini, Neotragini, and Tragelaphini had bootstrap supports of 55%, 56%, and 59%, respectively. The node describing the bifurcation between the single representative of Oreotragini in this analysis and the Cephalophini received similarly moderate support (53%). Cephalophini (71%) and Alcelaphini (82%) had relatively strong bootstrap support values. Most of the nodes describing relationships between the tribes and other deeper relationships (i.e., nodes for the subfamilies Bovinae and Antilopinae) also had low bootstrap support, while the highest levels of support were generally for nodes toward the tips. Only four of the 13 fossil species in the morphological data matrix aligned with their most likely tribes (table 2). Sivacobus palaeindicus, Miotragocerus amalthea, and Damalacra acalla placed as sister taxa to their most likely tribes. Sivatragus bohlini, although placed outside the tribe where it would be expected (Hippotragini), had the highest support for a fossil species based only on morphological data (62%), while the three species represented by ancient DNA mitochondrial genomes had strong support (84% Ovibos moschatus and Bootherium bombifrons, 99% for Bison priscus). The total evidence analysis of the 88 species for which morphological data were available resulted in the same relationships for fossil species as the full, 156 species analysis, thus only the phylogenies with more complete taxon sampling were considered for character optimizations.

FIGURE 1.

Total evidence maximum likelihood topology of Bovidae. Node labels are bootstrap supports over 50 and extinct species are labeled Ancient DNA or Fossil to denote how they appear in the matrix. Bootstrap values for tribes are displayed regardless of support level.

Several genera were not monophyletic on the total evidence phylogeny. Bison and Bos were not reciprocally monophyletic. The American bison (Bison bison) and the steppe bison (B. priscus) were each other's closest relatives, but linked to the sister taxa wild and domestic yaks (Bos mutus and Bos grunniens). Bootstrap support for this relationship was 92%. The gaur (Bos gaurus) was the sister taxon to this four-species clade, with 71% bootstrap support. The European bison (Bison bonasus) was sister taxon to a clade containing most of the rest of the species of Bos (except for the gayal, Bos frontalis), with 65% bootstrap support. In the Antilopinae, Sylvicapra grimmia nested within the genus Cephalophus, with moderate support (76%). Addax nasomaculatus as sister taxon to Oryx dammah had strong support (96%), but only in the analyses that included mitochondrial genomes; A. nasomaculatus only had mitochondrial DNA and morphological data in these analyses. Naemorhedus swinhoei nested within Capricornis, instead of with the other species from Naemorhedus, with strong support (98%). Nilgiritragus hylocrius was sister taxon to Ovis aries and Ovis orientalis, disrupting the monophyly of Ovis; however, this position was essentially unsupported at 32% bootstrap support. The total evidence analysis also placed Hemitragus jemlahicus within Capra, with similarly low support (33%). These genera also were not monophyletic on the most likely nuclear-DNA-only topology (fig. 2), with the exception of Oryx, which had data for fewer species and no data for A. nasomaculatus in the matrix.

The analyses of separate data partitions had different patterns of support, usually across the whole tree. The most likely topology for the morphological data showed very little support above the genus level; most internal nodes had no bootstrap support (fig. 3). This morphological topology resolved few of the tribes as monophyletic. The species of Reduncini were in a monophyletic group, with low bootstrap support, and the Cephalophini were monophyletic with moderate support. The Alcelaphini and Hippotragini were combined into a single group of species with low support, and the Caprini nested within the Antilopini. Fossils were generally misplaced with respect to other species of their likely tribes. Kobus subdolus maintained a sister taxon relationship with Kobus leche, and Sivacobus palaeindicus nested within the genus Kobus, albeit with low bootstrap support, rather than outside Reduncini + Pelea as it did in the total evidence phylogeny. Ugandax gautieri did not form a clade with other species of Bovini, and Aepyceros shungurae was not sister taxon to Aepyceros melampus. Turcoceros grangeri was sister taxon to Antilope cervicapra at the base of a clade containing species of both Antilopini and Caprini, closer to its likely tribe (Caprini), than it was on the total evidence topology. Miotragocerus amalthea in the morphological analysis was sister taxon to Tragelaphus scriptus; there was neither a monophyletic Tragelaphini nor a monophyletic Boselaphini on this topology.

TABLE 2.

Fossil species and tribes. Comparison of the most likely tribal affinity for each fossil species to its placement in the maximum likelihood total evidence phylogeny.

In the separate molecular data analyses, tree structure was fairly similar, and most tribes were monophyletic. All tribes of the subfamily Bovinae were monophyletic in the molecular data partitions. However, the nuclear data (fig. 2) resolved Pseudoryx nghetinhensis, the saola, as the sister taxon to Bubalina (the subtribe of water buffalo and African buffalo), while the mitochondrial data (fig. 4) placed it as sister taxon to all other species of Bovini. Nuclear data did not support a monophyletic Cephalophus, while mitochondrial DNA did, but placed Sylvicapra grimmia within the genus. Similarly, the nuclear DNA data did not support monophyly for Madoqua, Gazella, Nanger (Nanger granti, in particular, was grouped with Reduncini) in Antilopini, Damaliscus in Alcelaphini, nor Capra, Rupicapra, and Ovis in Caprini. Neither the nuclear nor the mitochondrial data recovered a monophyletic Bison. The mitochondrial phylogeny placed A. nasomaculatus within the genus Oryx, as in the total evidence analysis, but this placement could not be compared with the nuclear data analysis because A. nasomaculatus had no nuclear data represented in that matrix. Similar to the total evidence analysis, Naemorhedus swinhoei was located in the genus Capricornis in the mitochondrial DNA analysis. Relationships between tribes were generally similar in both the nuclear and mitochondrial analyses, except for the position of Oreotragus. It was sister taxon to Caprini + Hippotragini + Alcelaphini in the nuclear analysis, but sister taxon to Cephalophini in the mitochondrial analysis. The analysis of both molecular partitions as a single supermatrix (fig. 5) resulted in a topology that was nearly identical to the mitochondrial topology, differing only in the addition of species not represented in the mitochondrial matrix. This analysis placed Kobus vardonii as sister taxon to Pelea + Reduncini, with 94% bootstrap support.

In the dated Bayesian total evidence analysis (fig. 6), the major relationships between the tribes were unchanged from the maximum likelihood analysis. In most cases, fossil species were located well outside of their expected positions, with eight placing outside crown Bovidae entirely. Of the species that were located within extant tribes, M. amalthea, U. gautieri, and S. palaeindicus placed within their expected tribes, and D. acalla was located in Hippotragini instead of Alcelaphini, as in the maximum likelihood tree. S. bohlini was placed as the sister taxon to Neotragini in this topology, far outside the expected relationship to Hippotragini. As with the maximum likelihood phylogeny, support values were low, especially for the deepest relationships on the tree. Because this Bayesian analysis resulted in particularly unreasonable placements for the fossil taxa but otherwise produced the same topology for the extant tribes as estimated in the maximum likelihood analysis, it was not used for character optimizations.

Synapomorphies

Ancestral state reconstructions on both the total evidence topology (referred to as the full topology) and the topology with fossils located outside their likely tribes pruned (referred to as the reduced topology) produced a large number of characters that mapped as synapomorphies for tribes, subfamilies, and Bovidae (tables S4.1–4.6 in the supplement, available at  https://doi.org/10.5531/sd.sp.49). For the purpose of identifying characters that can link fossils to their most likely tribes within Bovidae, fossil comparisons were made for any character that met four criteria: the character mapped as a synapomorphy for that tribe (i.e., characters with a reconstructed state for the common ancestral node of a tribe that differed from the node immediately preceding it), it was coded for more than one species in the tribe, it had the same character state as most of those species, and it differed from the character state reconstructed for Cervidae (table 3). Not considered were characters that mapped as synapomorphies for Oreotragus oreotragus, as it was the only member of its tribe represented in this matrix, and thus its synapomorphies were indistinguishable from autapomorphies for that species. A larger number of characters mapped as synapomorphies on the reduced topology, but many of these characters were based on a single taxon or were not representative enough of the majority of the tribe to aid in the identification of phylogenetic affinities of the fossils. The majority of characters that met the criteria for inclusion mapped as synapomorphies on both the full and reduced topology (table S4.5).

Aepyceros shungurae had three characters in common with the proposed synapomorphies from Aepycerotini: long horn cores (char. 1, state 2), presence of a postcornual fossa (char. 12, state 1), and the presence of a transverse metaconid ridge on p3 (char. 55, state 1) (Gentry, 1992; Thomas, 1994). Gazella lydekkeri similarly had only one character that matched the synapomorphies of Antilopini: the position of the external auditory meatus' opening in lateral view is at or above the mandibular fossa (char. 184, state 1) (Bärmann, 2013, 2014). Two other characters could support inclusion of G. lydekkeri in Antilopini, a posterolingual cristid fusing the metaconid and entoconid on p4 (char. 133, state 1) and lack of contact between the nasal and lacrimal bones (char. 151, state 0), but these characters were also synapomorphies for Cephalophini and Hippotragini, respectively. Pachytragus crassicornis and Plesiaddax depereti both shared the presence of a flattened lateral surface of the horn core (char. 112, state 0) (Gentry, 1992) and mediolaterally compressed horn cores as a synapomorphy of Caprini (char. 196, state 1) (Gentry, 1992; Thomas, 1994; Bärmann, 2013, 2014). P. depereti also shared an indistinct paraconid on p3 and p4 (char. 56, state 0) (Gentry, 1992; Thomas, 1994; Bärmann, 2013, 2014) with Caprini, but this character also mapped as a synapomorphy for Cephalophini and Neotragini. Two characters matched between T. grangeri and the synapomorphies for Caprini: a roughly triangular-shaped tympanic bulla in lateral view (char. 179, state 1) and laterally directed opening of the external auditory meatus in ventral view (char. 183, state 0) (Bärmann, 2013, 2014). Two synapomorphies for Caprini that matched with T. grangeri, a third lobe of m3 consisting only of a hypoconulid (char. 137, state 2) and horn cores that diverge in frontal view (char. 193, state 1) (Bärmann, 2013, 2014), also mapped as synapomorphies for Boselaphini and Reduncini, respectively. Ugandax gautieri matched the synapomorphies for Bovini with just the presence of horn cores that are notably inclined in lateral view (char. 6, state 1) (Gentry, 1992; Thomas, 1994). Kobus subdolus and Sivacobus palaeindicus both have the reduncine synapomorphy of a braincase that is approximately the same width along the anteroposterior axis (char. 109, state 1) (Gentry, 1992). No characters mapped as a synapomorphy and met the inclusion criteria for Tragelaphini, and no other fossil species shared synapomorphies from either topology with their tribes.

FIGURE 2.

Nuclear DNA maximum likelihood topology. Node labels are bootstrap supports over 50.

Characters that mapped as synapomorphies for Bovidae on the full topology fit the criteria for inclusion better than those from the reduced topology. Many character reconstructions from the reduced topology that mapped as synapomorphies were based on few coded species and did not match with observed character states, and so were not considered in these fossil comparisons. In addition, the synapomorphic characters for Bovidae on the full topology were all coded for the cervid outgroups, lending more credence to their value for identifying fossils as bovids. Seven characters mapped as synapomorphies for Bovidae on the full topology, including three that matched to both G. lydekkeri and T. grangeri and two that matched to only G. lydekkeri among the fossil species. The three bovid synapomorphies found in G. lydekkeri and T. grangeri were a P2 of equivalent size to P3 (char. 119, state 2), a gap between the nasal and maxilla bones (char. 150, state 0), and horn cores that are centered on a point above the posterior margin of the orbit in lateral view (char. 187, state 0) (Bärmann, 2013, 2014). G. lydekkeri also had mesodont molars (based on the hypsodonty index of m3; char. 138, state 1), and a tympanic bulla that is inflated to the level of the occipital condyles (char. 176, state 1) (Bärmann, 2013, 2014). No fossil species had character 45, state 1, intermediate-sized foramina ovalia (Gentry, 1992; Thomas, 1994; Bärmann, 2013, 2014), or character 154, state 1, presence of only one lacrimal foramen on each side of the cranium (Bärmann, 2013, 2014), the last two bovid synapomorphies identified in this analysis.

FIGURE 3.

Morphological data maximum likelihood topology. Node labels are bootstrap supports over 50. Extinct species are labeled Fossil.

FIGURE 4.

Mitochondrial DNA maximum likelihood topology. Node labels are bootstrap supports over 50. Extinct species are labeled Ancient DNA.

TABLE 3.

Synapomorphic characters for Bovidae and its tribes.

Topology

Systematic positions of fossil species were generally poorly supported on the maximum likelihood total evidence topology, likely due to the lack of clear phylogenetic signal in the morphological data matrix. This discussion will focus on the maximum likelihood total evidence topology because it was used for morphological character optimizations. Of the four species that aligned with their most likely tribes, their positions within the tribes were not well supported. Ugandax gautieri often aligns with the African and water buffaloes (Bibi, 2007), while here it was the sister taxon to the Bos + Bison clade. Plesiaddax depereti likewise was located in the wrong subtribe of Caprini. Expected to align with Ovibovina (Gentry, 1992), it instead nested in Caprina in the genus Capra. In the Reduncini, Kobus subdolus was sister taxon to Kobus leche, and had the second-highest bootstrap support of the relationships within Kobus (although support was low throughout the genus), in contrast to other analyses that place K. subdolus in Redunca (Vrba, 2006; Vrba and Haile-Selassie, 2006).

Among the fossils placed outside, but near, their likely tribes, Miotragocerus amalthea offers a tantalizing hint for revising relationships within the subfamily Bovinae. Whether fossil boselaphines actually belong in the tribe is unclear (Bibi, 2007; Bibi et al., 2009); M. amalthea as the sister taxon to the bovine Pseudoryx nghetinhensis, the saola, had higher support relative to the other fossil species' placements (albeit still low at 16%). The saola itself is morphologically disparate compared with other bovines (Gatesy and Arctander, 2000), and it is possible that the two are linked here because they are derived relative to other members of Bovinae. However, M. amalthea had more characters in common with P. nghetinhensis (12 of the 31 characters encoded for it) than with either living boselaphine or most other bovines. If the position of M. amalthea is robust to the addition of more morphological characters and more fossil boselaphines, it could help fill in the history of morphological changes that distinguish the enigmatic saola from other members of Bovini.

FIGURE 5.

Topology generated from combined maximum likelihood analysis of nuclear and mitochondrial DNA. Node labels are bootstrap supports over 50. Extinct species are labeled Ancient DNA.

Of the nonmonophyletic genera on the total evidence topology, several are well documented to lack monophyly when analyzed with certain data. Lack of reciprocal monophyly for Bison and Bos on this total evidence phylogeny and each separate analysis of the molecular data is a common pattern in phylogenetic analyses. Multiple studies have found evidence for introgression among species of these two genera (Verkaar et al., 2004; Hassanin et al., 2013; Wu et al., 2018). Given the consistent lack of monophyly of Bison and Bos across multiple studies, it is clear that a thorough reassessment of these genera drawing on the many fossil species not included here and a larger set of morphological characters is required. Such an analysis could potentially support the inclusion of Bison as a junior synonym of Bos. Sylvicapra could also be considered a junior synonym of Cephalophus, and S. grimmia could be transferred to Cephalophus, given its well-supported placement within that genus on this phylogeny and in previous research (Hassanin et al., 2012; Johnston and Anthony, 2012; Bibi, 2013).

The relationship between Oryx dammah and Addax nasomaculatus, although strongly supported in this analysis, is contradicted by other studies (Bibi, 2013; Hassanin et al., 2013). Fertile hybrids between O. dammah and A. nasomaculatus are possible (Engel, 2004; Gilbert, 2017), which could support an interpretation that their pairing here is evidence for introgression captured in the molecular data. Monophyly for various clades within Caprini (the tribe containing Capricornis, Naemorhedus, Ovis, Capra, Nilgiritragus, and Hemitragus) has also been inconsistent in past studies (Hassanin et al., 1998, 2012). However, this total evidence topology has more nonmonophyletic genera than previous studies, which mostly agree on the placement of H. jemlahicus as sister taxon to C. sibirica (Hassanin et al., 1998, 2012; Bibi, 2013). Although previous research suggests that Neotragus is a polyphyletic genus (Bärmann and Schikora, 2014), the total evidence topology presented here does not support that interpretation; all three members of the genus form a single clade with moderate bootstrap support.

The relationships of tribes in Bovinae were similar in this analysis to previous molecular phylogenies: Tragelaphini and Bovini were sister taxa, with Boselaphini as the sister taxon of that pairing (Marcot, 2007; Decker et al., 2009; Hassanin et al., 2012; Chen et al., 2019). This grouping is contradicted in another analysis of combined morphological and molecular data with more limited taxonomic sampling (Gatesy and Arctander, 2000), as well as two of the morphological analyses that provided characters analyzed in this study (Gentry, 1992; Thomas, 1994), which group Bovini with Boselaphini as sister taxa instead. The extant species of Boselaphini are morphologically disparate, and their relationship to fossil boselaphine species is unclear (Bibi, 2007; Bibi et al., 2009). Although the analysis presented here agreed with the placement of Boselaphini as the sister taxon to the rest of Bovinae, including the other fossils attributed to the tribe and subfamily and a more robust character set in future analyses may help determine which position for the tribe should be favored.

Tribal relationships within Antilopinae differ greatly in this analysis from those recovered by previous molecular and supermatrix studies. Previous analyses of molecular data supported a Hippotragini + Alcelaphini clade with Caprini as its sister taxon, similar to the one found in this total evidence topology (Marcot, 2007; Hassanin et al., 2012). The sister taxon to the Hippotragini + Alcelaphini + Caprini in a supermatrix analysis of nuclear and mitochondrial genes was Antilopini (Marcot, 2007), but an analysis of mitochondrial genomes alone placed Cephalophini + Oreotragini as the sister taxon to this clade (Hassanin et al., 2012). The total evidence analysis in this paper proposes yet another grouping of these tribes. The Hippotragini + Alcelaphini + Caprini clade is linked to Reduncini, as in Bibi's (2013) mitochondrial DNA analysis. Oreotragini + Cephalophini is the sister clade to Antilopini on the total evidence topology, which in turn is sister taxon to the clade containing Hippotragini, Alcelaphini, Caprini, and Reduncini. The position of Aepycerotini and Neotragini as a sister clade to the rest of the tribes in Antilopinae, however, is found in most recent large-scale molecular analyses and may simply be a case of two highly derived families clustering with each other (Hassanin et al., 2012; Bibi, 2013; Chen et al., 2019). Including fossil species from either tribe in a future analysis would be a good test of whether fossil species represented only by morphology in total evidence analyses can accomplish the increased taxon sampling required to break up long branches (Hillis, 1996) and support a different position for these families.

The similarities of the total evidence phylogeny to other maximum likelihood analyses of mitochondrial genomes or combined mitochondrial and nuclear data (e.g., Hassanin et al., 2012; Bibi, 2013; Bärmann et al., 2013) suggest that the mitochondrial genomes have a strong effect on topology, which may be unsurprising, given that they formed the largest partition in the total evidence matrix. The absence of mitochondrial data for several species, however, did not seem to preclude placing those species with the correct tribe or genus. Instead, the differences in the availability of nuclear data may have been more determinative. More nuclear genes were available for Bison bonasus than for Bison bison; half the species of Bos that grouped with B. bonasus had sequences for AMELX, while none of those that grouped with B. bison did. For the Cephalophini, C. zebra is represented in the matrix only by sequence data for TSHB, while S. grimmia and nearly all other species of Cephalophus also had PRKCI sequences and mitochondrial genomes. For the other nonmonophyletic genera, however, there were no clear patterns between species represented by different nuclear genes and the lack of monophyly.

The lack of morphological data does not appear responsible for the misplacement of species, as multiple extant species that were present only in the morphological matrix (e.g. Neotragus pygmaeus, Procapra picticaudata, Gazella marica) nonetheless joined their congenerics in monophyletic groups, albeit with lower bootstrap support values. Low bootstrap support was likely related more to conflicting phylogenetic signal in the nuclear and mitochondrial data or the nearly complete lack of signal in the morphological data than to the presence of species with missing data. Missing data alone do not prevent accurate placement of species in a phylogeny when there are many characters in the matrix (Wiens, 2003; Wiens and Moen, 2008; Wiens and Morrill, 2011). Within total evidence analyses, missing morphological data can influence recovery of the “correct” topology (Sansom and Wills, 2013, 2017; Pattinson et al., 2015; Guillerme and Cooper, 2016); however, the inclusion of morphological data for at least half the living species in the total evidence matrix can remedy this issue (Guillerme and Cooper, 2016). The total evidence matrix here includes morphological data for 75 extant species (52% of all living species included in the matrix). Different types of morphological characters (i.e., dental, osteological, soft tissue) display varying levels of homoplasy within mammals, and thus may provide conflicting evidence of phylogenetic relationships when some types of data are missing for fossil species (Sansom and Wills, 2013, 2017; Sansom et al., 2017). Coding more morphological characters for more extant taxa could help alleviate this issue, especially by filling in gaps where the characters available for fossil species may not have been coded for living species to which they are related. The extant species correctly placed on the basis of morphological characters alone suggests that the failure of some fossil species to align in this study with their most likely tribes may be attributed to the lack of close relatives and evolutionary distance between them and extant representatives of their tribes.

FIGURE 6.

Maximum clade credibility Bayesian total evidence topology (opposite page, above, and next page). Extinct species are labeled Ancient DNA or Fossil to denote how they appear in the matrix. Posterior probabilities for each node over 0.5 are represented as scaled circles. Divergence dates estimated from fossil calibrations are represented in millions of years. Note the estimated divergence time of the clade Oreotragini + Cephalophini + Antilopini was older than the estimated divergence of the rest of Antilopinae in which it belongs.

Continued

Continued

Synapomorphies

The ultimate goal of this study was to find hard-tissue synapomorphies for the family Bovidae by optimizing morphological characters on a total evidence phylogeny, testing whether a total evidence phylogeny could provide synapomorphies for the clades that are of use in determining fossil relationships to Bovidae. Historically, the defining synapomorphy for Bovidae has been the presence of permanent, unbranched horn cores covered by permanent, unbranched keratin sheaths (Gray, 1821; Pilgrim, 1946; Janis and Scott, 1987; Bibi et al., 2009). The lack of preservation of horn sheaths in the fossil record, the great amount of homoplasy in the family, and the variable presence of horns in females have made it difficult not only to place fossil species in the family, but also to expand the list of diagnostic characters beyond the presence of permanent, unbranched horn cores (Bibi et al., 2009). An analysis by Pilgrim (1946) aligned Eotragus with Bovidae based on the presence of horns cores, selenodont molars with four distinct crescents, upper premolars with a single enamel crescent, and complex enamel folds on the lower premolars, which he proposed distinguish Bovidae from Gelocidae (Pilgrim, 1946).

Janis and Scott (1987) found several derived characters for Bovidae in their analysis, but rejected them as synapomorphies because of homoplasy with other ruminant families. These characters included hypsodont molars and elongated limbs without lateral digits, which are common among other nonbovid ruminants in open habitats; “fused metatarsals with an open gully,” shared with members of Giraffoidea; and entostyles on upper molars, which bovids share with brachydont members of Cervoidea (Janis and Scott, 1987). Although the shape and position of the auditory bulla may be diagnostic for some bovids, they noted that many species do not have an expanded bulla at all; based on their interpretation of these characters and the scope of their morphological data, they determined that the presence of horns was the only diagnostic character of Bovidae (Janis and Scott, 1987).

Gentry's (1992) analysis of morphological characters produced several synapomorphies for tribes and other subclades of Bovidae, but the topology proposed by his analysis diverges strongly from the topology produced by molecular data; as a result, comparing his tribal level synapomorphies to the results of this study presents several challenges. On his topology, the clade most resembling Bovidae for which he reports synapomorphies excludes S. grimmia, T. scriptus, and T. quadricornis. Seven of his characters map as synapomorphies for this clade: an inclined cranial roof, a braincase with parallel sides or sides that widen anteriorly, large postglenoid foramen, foramina ovalia opening laterally, lack of basal pillars on all upper and some lower molars, a large indent on the femur between the articular head and the greater trochanter, and distal flanges on the radius (Gentry, 1992).

My analysis produced two sets of characters that map as bovid synapomorphies: one set based on the full topology, and one on the topology with incorrectly placed fossil taxa removed. I found the full-topology synapomorphies, although fewer than those from the reduced topology, conformed better to the character states found in the data matrix. These characters were useful for aligning only two of the fossil species, G. lydekkeri and T. grangeri, with Bovidae because they were not scored for the other fossil taxa in this matrix; however, their broad distribution among the extant species of Bovidae supports their value for future studies.

Two of the synapomorphies identified here share some similarities with the hypsodont molar and anteriorly positioned tympanic bulla that Janis and Scott (1987) proposed as potential diagnostic traits for Bovidae. The mesodont molar synapomorphy of my analysis may simply be a different way of defining and coding one of these characters; their paper distinguishes between only brachydont and hypsodont crown heights, defining any tooth with crown height greater than the tooth width as hypsodont (Janis and Scott, 1987), whereas in the data analyzed here, mesodonty is judged by a range of hypsodonty index values from 1.5 to 3 (Bärmann, 2013, 2014). Their broader definition of hypsodonty thus encompasses the mesodonty of this paper, and potentially obscures a character that could serve as a synapomorphy. Hypsodonty is rightly recognized a highly plesiomorphic character, thus using this synapomorphy for distinguishing between bovids and other ruminants requires great care (Janis and Scott, 1987; Janis and Theodor, 2014). I recommend considering hypsodonty only in combination with other synapomorphies for the family. The auditory bulla character state that mapped as a synapomorphy in this analysis also differentiated between more levels of inflation than Janis and Scott (1987). Although it is consistent across many of the species in this matrix, the character is not coded for several tribes, and so should be treated as a tentative synapomorphy from this analysis.

Synapomorphies from this study could link fossil species with Aepycerotini, Antilopini, Bovini, Caprini, and Reduncini, out of the nine tribes with fossil representatives. Long, uncompressed, lyre-shaped horn cores and a large, deep postcornual fossa were among the characteristics that diagnose Aepyceros, the only genus in Aepycerotini (Gentry, 2011). These characters both mapped as synapomorphies for the tribe and supported including A. shungurae in the clade. A transverse metaconid ridge on p3 (char. 55, state 1), which mapped as a synapomorphy in this study, could expand the list of features that define Aepycerotini. Gentry (1992) found few monophyletic tribes that agree with large-scale molecular or combined-evidence analyses, but his analysis did reconstruct monophyletic Antilopini. One character, fusion of the metaconid and entoconid of p4, was a synapomorphy on both his morphological phylogeny and my total evidence phylogeny, which could link G. lydekkeri to Antilopini. This character was not exclusive to Antilopini on this phylogeny, thus it requires corroboration from other characters when assessing fossil relationships. Geraads (1992) proposed an inclined horn core synapomorphy for Bovini. Although the character differs in describing the inclination in my analysis, it linked U. gautieri to Bovini. This synapomorphy, along with a suite of other horn core characters, also supports the inclusion of the genus Brabovus, not analyzed here, in Bovini (Geraads, 1992). Short metacarpals are a traditional character defining Caprini (Bibi et al., 2012). The analysis in this study found several new craniodental synapomorphies that could be used for fossil placements, and could especially benefit attempts to reassess the large number of fossil caprines potentially incorrectly referred to the family (Gentry, 2000). Previous research identifies at least four cranial and dental characters that support Reduncini: the high anterior angle of the horn cores, the low angle of the face relative to the braincase in lateral view, prominent “goat fold” anterior cingulids (sensu Bärmann and Rössner, 2011) on the molars, and prominent ectostyles/-stylids (Vrba et al., 1994). The character identified in this study that could link K. subdolus and S. palaeindicus to Reduncini, a braincase the same width along the anteroposterior axis, thus represents a new potential synapomorphy for the clade.

Although character optimizations here produced synapomorphies throughout Bovidae, they must be tested against other fossil species in future analyses. For example, it is unclear whether some species of Eotragus, among other purported early bovids, are stem or crown bovids (Solounias et al., 1995; Bibi, 2007; Bibi et al., 2009), and assignment of some fossil species to certain tribes are contradicted by different analyses. There are many fossil taxa that could benefit from inclusion in a large-scale analysis with a goal toward refining their position in Bovidae, such as Brabovus nanincisus, which could belong in Bovini, Tragelaphini, or Cephalophini (Gentry, 2011). Likewise, using these characters in efforts to revise the positions of taxa currently considered stem species of Moschidae or other ruminant families could test the value of these characters as bovid synapomorphies. Recent studies of these fossils provide rich sets of morphological characters, unfortunately few of which overlap with the morphological synapomorphies in this study. The presence of P2 and P3 of equivalent size (char. 119, state 2) may be found in Propalaeoryx, considered part of the Giraffoid family Climacoceratidae, and Namacerus, a likely stem bovid (Morales et al., 2003). Presence of a single lacrimal foramen (char. 154, state 1) is shared by several fossil species (e.g., Dremotherium, Amphimoschus, Ampelomeryx) from across the ruminant tree, but not cervids and antilocaprids (Mennecart et al., 2021). In a future analysis of these characters in the broader context of ruminants, this character may not optimize as a synapomorphy of Bovidae, but could instead serve as morphological support for the relationship between Bovidae and Moschidae. Two synapomorphies for Bovidae, inflation of the tympanic bulla to the level of the occipital condyle (char. 176, state 1) and mesodont molars (char. 138, state 1), are similar to descriptions of fossils in recent studies, but require clarification before they could be used to support inferences of relationships. Sperrgebietomeryx, likely part of the Climacoceratidae, may have an inflated tympanic bulla (Morales et al., 2008). Namacerus is described as having “slightly hypsodont” molars (Morales et al., 2003), and Namibiomeryx exhibits “incipient hypsodonty” (Morales et al., 2008). Reexamining these species with quantitative versions of these characters (e.g., defining hypsodonty based on a range of hypsodonty index values) could better elucidate the anatomy's variation, evolution, and value as a diagnostic synapomorphy. By examining these species for synapomorphies identified in this study and expanding the dataset to include these species in further phylogenetic analyses, we could achieve a more complete picture of bovid, and even ruminant, evolutionary history.