Inferring the phylogeny of bovidae using mitochondrial DNA sequences: resolving power of individual genes relative to complete genomes.

Molecular techniques that assess biodiversity through the analysis of a small segment of mitochondrial genome have been getting wide attention for inferring the mammalian diversity. Due to their highly conserved nature, specific mitochondrial genes offer a promising tool for phylogenetic analysis. However, there is no established criteria for selecting the typical mitochondrial DNA (mtDNA) segments to achieve a greater resolving power. We therefore chose the family Bovidae as a model and compared the tree-topologies resulting from the commonly used and phylogenetically-informative genes including 16S rRNA, 12S rRNA, COI, Cyt b and D-loop with respect to complete mitochondrial genome. The tree topologies from the whole mitochondrial genome of 12 species were not identical albeit similar with those resulting from the five individual genes mentioned above. High bootstrap values were observed for mtDNA compared with that of any single gene. The average pair-wise sequence divergence using different genetic modes was found to be: D-loop (0.229) > Cyt b (0.159) > COI or complete mtDNA (0.143) > 12S rRNA (0.094) > 16S rRNA (0.091). The tree resulting from complete mtDNA clearly separated the 12 taxa of Bovidae into 3 major clusters, one cluster each for subfamily Cervinae and Bovinae and the third cluster comprised the distinctive clades of Caprinae and Antilopinae. However, jumping clades of Antilopinae were observed while using the individual genes. This study showed that Bison bison and Bos Taurus have very close phylogenetic relationship compared to Bubalus bubalis (Bovinae), irrespective of the method used. Our findings suggest that complete mtDNA genome provides most reliable understanding of complex phylogenetic relationships while the reliability of individual gene trees should be verified with high bootstrap support.

Introduction

Members of the family bovidae (order: Artiodactyla) include bison, buffalo, antelopes, gazelles, sheep, goats, muskox and domestic cattle, which are distinguished by the presence of permanent hollow horns. The family bovidae has a great variety of morphologies with 137 living and more than 300 fossil species have been described.1 The phylogenetic relationships and taxonomy of this family have been controversial for a long time. Undisputed divisions of bovids include (i) Bovinae (for example cattle, nilgai and eland), (ii) Cephalophinae (duikers), (iii) Caprinae (sheep, goats and related animals), (iv) Hippotraginae (roan antelope) and (v) Antilopinae (gazelles, chiru and blackbuck).2 Specifically, there is only one morphological character that unambiguously defines the bovids: their non-deciduous horn cores and horn sheaths.3 More than half a century ago, the systematics of Bovids was extremely difficult and Bovidae was considered as one of the most troublesome groups of mammals to classify.4 However, the new advents in sequencing analysis and bioinformatics have simplified the molecular systematics of Bovidae to some extent. The application of mitochondrial genes such as 12S rRNA, Cytochrome b (Cyt-b) and displacement loop (D-loop) has been getting wide interest in phylogenetic analysis of diverse taxa.5–8

Mitochondrial DNA (mtDNA) has a relatively fast mutation rate, which results in significant variation in mtDNA sequences between species and in principle, a comparatively small variance within species.9 Mitochondrial 16S rRNA gene sequence has been used for the identification of 182 vertebrates and 103 invertebrates, while a single locus appeared to be sufficient for the identification of most of the species.10 Recently, the reliability of mitochondrial gene barcodes has been determined for diverse clades of birds.11–15 Fernandez & Vrba16 have provided a complete estimate of phylogenetic relationship among the members of Ruminantia. However, the efficiencies of single mitochondrial genes versus complete mitochondrial genome for inferring complex phylogenies have not been thoroughly investigated. In this study, we addressed the question whether the individual mitochondrial genes (16S rRNA, 12S rRNA, COI, Cyt-b and D-loop) provide the same phylogenetic information as compared to complete mitochondrial DNA (mt-DNA). We used 12 representative species from the mammalian family Bovidae to test our hypothesis.

Materials and Methods

We retrieved the nucleotide sequences of the whole mitochondrial genome as well as the individual mtDNA gene sequences of 16S rRNA, 12S rRNA, CO1, Cyt-b and D-loop of 12 members of the family Bovidae from the GenBank database. The details of these sequences are given in Table 1. The sequences were aligned by using CLUSTAL-X software (http://www.clustal.org), version 2.0.12.17 The sequences were then trimmed to get their equal lengths for all the species. As a result, a total of 15864 nucleotide positions for mtDNA, 1527 bp for 16S rRNA gene, 704 bp for 12S rRNA, 1545 bp for CO1, 1140 bp for Cyt-b and 806 bp for D-loop were used in the final dataset.

Table 1

Sequences of the animal species used in the comparative phylogeny of the family Bovidae.

Species	Subfamily	Accession no.	Size (bp)
			mtDNA	16S rRNA	12S rRNA	CO1	Cytb	D-loop
Ammotragus lervia	Caprinae	NC009510	16530	1569	958	1545	1140	1098
Capra hircus	Caprinae	GU068049	16642	1572	571	1545	1140	1212
Budorcas taxicolor	Caprinae	FJ006534	16667	1574	955	1545	1140	1235
Ovis aries	Caprinae	AF010406	16616	1574	958	1545	1140	1180
Capricornis crispus	Caprinae	AP003429	16453	1568	959	1545	1140	1022
Naemorhedus caudatus	Caprinae	FJ469673	16519	1557	957	1545	1140	1099
Pantholops hodgsonii	Antilopinae	NC007441	16498	1566	957	1545	1140	1067
Antilope cervicapra	Antilopinae	AP003422	16431	1570	957	1545	1140	998
Cervus nippon	Cervinae	AB211429	16663	1570	955	1545	1140	1229
Bison bison	Bovinae	EU177871	16319	1570	956	1545	1140	888
Bos taurus	Bovinae	DQ124418	16340	1571	956	1545	1140	910
Bubalus bubalis	Bovinae	AY702618	16359	1569	955	1545	1140	910

The evolutionary history was inferred using the maximum likelihood18 and neighbor-joining19 methods. All the phylogenetic analyses were conducted in MEGA software (http://www.megasoftware.net), version 4.20 Other software tools such as PHYLIP, PAUP, HyPhy, etc. may also be used to conduct phylogenetic analysis using gene sequence data. The percentage of replicate trees in which the associated taxa clustered together was determined by the bootstrap test (1000 replicates).21 Bootstrapping is a commonly used method for constructing reliable trees by subsampling from the sites in an alignment to create trees based on sub-samples. The process is iterated multiple times (preferably 1000 times) and the results are compiled to allow an estimate of the reliability of a particular grouping.22 Estimates of evolutionary divergence between sequences were determined using the maximum composite likelihood method in MEGA.23 A composite likelihood is defined as a sum of related log-likelihoods. Since all pair-wise distances in a distance matrix have correlations due to the phylogenetic relationships among the sequences, the sum of their log-likelihoods is a composite likelihood. This model assumes equality of substitution pattern among lineages and of substitution rates among sites. We also conducted maximum parsimony and minimum evolution models to verify the bootstrap supports and these findings have been shown as nodal dots in NJ trees to avoid repetition by showing the same trees for these two methods. Pair-wise base homology (%) was determined by using the formulae: (1 −evolutionary divergence between sequences) × 100. All the positions containing gaps or missing data were eliminated (complete deletion) from the dataset prior to analysis. However, MEGA software also provides alternatives to retain all such sites initially and excluding them as necessary in the pair-wise distance estimation (pairwise deletion option) or to use the partial deletion (site coverage) as a percentage.

Results and Discussion

The pair-wise sequence diversity was found to be lowest for 16S rRNA gene and highest for D-loop. The average sequence diversities for various gene segments were as follows: 16S rRNA (average 0.091, range 0.04–0.13), 12S rRNA (0.094, 0.03–0.13), mtDNA genome (0.143, 0.06–0.17), CO1 (0.143, 0.05–0.18), Cyt-b (0.159, 0.08–0.20) and D-loop (0.229, 0.10–0.32). Conversely, in terms of pair-wise sequence similarity, it was highest for 16S rRNA gene and the lowest for D-loop (Fig. 1). The overall pair-wise sequence similarities were 90.9% (range, 82.2%–96.4%), 90.6% (87.2%–97.4%), 85.7% (82.5%–93.8%), 85.5% (82.0%–94.8%), 84.1% (79.8%–92.4%), 77.1% (68.4%–89.6%) for 16S rRNA, 12S rRNA, mtDNA genome, CO1, Cyt-b and D-loop, respectively. The phylogenetic trees constructed using the whole mitochondrial genome appeared to be identical irrespective of the method used (ML versus NJ) (Figs. 2 and 3). The phylogenetic trees resulted from the sequences of individual genes (16S rRNA, 12S rRNA, CO1, Cyt-b, D-loop) were not identical with that of whole mitochondrial genome based tree; their topologies were only partially similar for both ML (Fig. 2) and NJ (Fig. 3) methods.

Figure 1

Comparative view of average pair-wise sequence similarity among 12 members of Bovidae family using individual mitochondrial genes and nearly complete mitochondrial genome. Vertical error bars show the standard error of mean.

Figure 2

Maximum likelihood dendrograms showing the phylogenetic relationship among members of the family Bovidae based on almost complete nucleotide sequences of (A) mtDNA (B) 16S rRNA (C) 12S rRNA, (D) CO1 (E) Cyt-b and (F) D-loop.

Notes: Bootstrap values (expressed as percentages of 1000 replications; >50%) are shown at branching points. Corresponding GenBank/EMBL/DDBJ accession numbers are written in the parentheses. Bars represent 1 substitution per 200 nucleotides (A, D, E), 100 nucleotides (B, C) and 500 nucleotides (F).

Figure 3

Neighbor-Joining dendrograms showing the phylogenetic relationship among members of the family Bovidae based on almost complete nucleotide sequences of (A) mtDNA genome, (B) 16S rRNA, (C) 12S rRNA, (D) CO1, (E) Cyt-b and (F) D-loop.

Notes: Bootstrap values (expressed as percentages of 1000 replications; >50%) are shown at branching points. Corresponding GenBank/EMBL/DDBJ accession numbers are written in the parentheses. Bars represent 1 substitution per 200 nucleotides (A, D, E, F) or 100 nucleotides (B, C). Nodes marked with solid circle are supported by >50% bootstrap values when analyzed by maximum parsimony and minimum evolution methods.

The phylogenetic analyses using complete mt-DNA indicated main split of the 12 members of the Bovidae into one bovine clade and one non-bovine clade, which grouped all other bovids (Figs. 2 and 3). The same broader cladistic was observed using 16S rRNA, 12S rRNA and CO1 genes but not with Cyt-b and D-loop. Further branching showed low bootstrap values (<50%) for the trees obtained using individual genes. An elaborated study using 197 species of the ruminants showed the same trends, splitting into bovine clade and non-bovine clade.16 The tree topologies unanimously showed very close relationship between Bison bison and Bos taurus under the subfamily Bovinae. Bootstrap values for the grouping of these species for all the trees were 100% (Figs. 2 and 3) with low sequence diversity (Supplementary Tables 1–6). Due to their close phylogenetic relationship, it has been suggested that Bos and Bison should be integrated into a single Bos genus.16 Bubalus bubalis was observed distantly related with Bos-Bison group as previously reported.7 Molecular data of Bovini suggest two lineages, buffalo (Bubalus sp.) versus cattle (Bos, Bison).2,24 Another report also indicated that Bison and Bos are more closely related to each other rather with Bubalus, as supported from morphological, paleontological, and reproductive data.5

Besides clearly differentiating the subfamily Bovinae, the complete mitochondrial genome tree clustered all the members of Caprinae with high bootstrap supports and sequential ancestry from the two members of Antilopinae. However, the placement of Antilopinae using the individual genes appeared to be unstable that can be seen as jumping clades in respective trees. Previous study has suggested the monophyly of Bovinae and Caprinae however Antilopinae appeared to be polyphyletic.16 All the tree analyses also showed close relationship between Capricornis crispus and Naemorhedus caudatus (subfamily: Caprinae) compared with the other members included in the tree analyses and supported by high bootstrap values (98%–100%). This corroborates with the previous reports.16,25,26

The systematic position within the tribe Bovini remains confused since the analyses of morphological characters have led to several conflicting hypotheses. Some authors have suggested that Bovinae could be comprised of hybrid species produced by the crossing of the banteng with gaur, zebu, or water buffalo.27 Systematic work on Bovids has been difficult and the Bovidae is one of the most troublesome groups of mammals to classify.4 The family Bovidae might have formed about 25 million years ago within a period of 5 million years.28 One of the most striking aspects of recent higher mammalian phylogeny indicated that artiodactyls possessed a common ancestor that is not shared by any other group.2 Genomic distances between the cattle species have been influenced by genetic exchange between neighboring ancestral populations.29 The intractability of this systematic problem is consistent with a rapid radiation of the major bovid groups6

The highly effective method for measuring support for phylogenetic relationships is boot-strapping.21 We examined whether a clade present in the tree constructed with single gene is present on the whole genome tree or not and the bootstrap value of the node. The most frequent within-method variations were related to the placement of Antilopinae. Between-methods comparison also showed variations in tree topologies resulting from all the individual genes except Cyt-b. Moreover, the trees generated from complete mtDNA genome showed greater resolution with high bootstrap support as compared to phylogeny inferred from individual genes. A large number of nucleotide sites are needed to exactly determine the whole-genome tree whereas a relatively small number of sites often results in a tree with closer topology.30 It has been shown that blocks of contiguous sites are less likely to lead to the whole-genome tree than samples composed of sites drawn individually from throughout the genome.30

In conclusion, the findings of this study showed that complete mitochondrial genome provides a greater resolution in phylogenetic analysis of complex taxonomic groups. The phylogeny of Bovidae using the sequences of individual genes (16S rRNA, 12S rRNA, CO1, Cytb and D-loop) of mtDNA failed to provide identical tree topology with that of complete mtDNA. The tree resulting from complete mtDNA clearly separated the 12 taxa of Bovidae into clusters with distinctive phylogeny however jumping clades of Antilopinae were observed while using the individual genes. The common phylogenetic inference using individual genes or complete mitochondrial genome was the placement of Cervinae and Bovinae. Thus, for understanding the complex phylogenetic relationships, the use of complete mitochondrial genome should be preferred over individual genes. Nevertheless the individual gene trees with conditional high bootstrap support may also provide useful phylogenetic information.

Table S1.

Estimates of evolutionary divergence between sequences of nearly complete mtDNA of the used taxa under the family Bovidae.

No.	Species	1	2	3	4	5	6	7	8	9	10	11
1	Ammotragus lervia	0
2	Capra hircus	0.10
3	Budorcas taxicolor	0.10	0.11
4	Ovis aries	0.10	0.11	0.11
5	Capricornis crispus	0.11	0.12	0.12	0.12
6	Naemorhedus caudatus	0.11	0.12	0.12	0.12	0.08
7	Pantholops hodgsonii	0.11	0.12	0.12	0.12	0.12	0.12
8	Antilope cervicapra	0.15	0.15	0.15	0.15	0.15	0.15	0.14
9	Cervus nippon	0.17	0.17	0.17	0.17	0.17	0.17	0.16	0.17
10	Bison bison	0.16	0.16	0.16	0.16	0.16	0.16	0.16	0.16	0.17
11	Bos taurus	0.16	0.16	0.17	0.16	0.16	0.17	0.16	0.17	0.17	0.06
12	Bubalus bubalis	0.16	0.16	0.16	0.16	0.16	0.16	0.16	0.16	0.17	0.13	0.14

Note: There were a total of 15,864 positions in the final dataset.

Table S2.

Estimates of evolutionary divergence between sequences of 16S rDNA of the used taxa under the family Bovidae.

No.	Species	1	2	3	4	5	6	7	8	9	10	11
1	Bison bison	0
2	Bos taurus	0.04
3	Bubalus bubalis	0.08	0.08
4	Capricornis crispus	0.11	0.11	0.11
5	Naemorhedus caudatus	0.12	0.12	0.11	0.04
6	Ammotragus lervia	0.11	0.11	0.11	0.05	0.06
7	Capra hircus	0.11	0.11	0.11	0.06	0.07	0.05
8	Budorcas taxicolor	0.09	0.10	0.10	0.06	0.07	0.05	0.05
9	Ovis aries	0.11	0.11	0.10	0.07	0.07	0.06	0.07	0.06
10	Pantholops hodgsonii	0.10	0.11	0.11	0.07	0.07	0.06	0.07	0.06	0.06
11	Antilope cervicapra	0.12	0.12	0.11	0.10	0.10	0.09	0.09	0.09	0.10	0.09
12	Cervus nippon	0.12	0.13	0.11	0.11	0.12	0.11	0.10	0.10	0.12	0.12	0.12

Note: There were a total of 1,527 positions in the final dataset.

Table S3.

Estimates of evolutionary divergence between sequences of 12S rDNA of the used taxa under the family Bovidae.

No.	Species	1	2	3	4	5	6	7	8	9	10	11
1	Bison bison	0
2	Bos taurus	0.03
3	Bubalus bubalis	0.08	0.08
4	Ammotragus lervia	0.09	0.09	0.09
5	Capra hircus 0.09	0.10	0.10	0.04
6	Budorcas taxicolor	0.11	0.11	0.10	0.05	0.07
7	Ovis aries 0.11	0.10	0.10	0.06	0.06	0.07
8	Pantholops hodgsonii	0.10	0.10	0.10	0.07	0.07	0.09	0.07
9	Capricornis crispus	0.10	0.10	0.10	0.05	0.06	0.08	0.07	0.08
10	Naemorhedus caudatus	0.12	0.12	0.11	0.08	0.08	0.09	0.07	0.08	0.04
11	Antilope cervicapra	0.12	0.13	0.11	0.11	0.11	0.11	0.10	0.12	0.12	0.12
12	Cervus nippon	0.12	0.13	0.11	0.10	0.10	0.12	0.11	0.13	0.12	0.12	0.13

Note: There were a total of 704 positions in the final dataset.

Table S4.

Estimates of evolutionary divergence between sequences of CO1 mtDNA of the used taxa under the family Bovidae.

No.	Species	1	2	3	4	5	6	7	8	9	10	11
1	Antilope cervicapra	0
2	Cervus nippon	0.16
3	Capra hircus	0.14	0.17
4	Ovis aries	0.16	0.17	0.10
5	Ammotragus lervia	0.17	0.17	0.11	0.11
6	Budorcas taxicolor	0.16	0.17	0.10	0.12	0.11
7	Capricornis crispus	0.16	0.17	0.12	0.12	0.13	0.11
8	Naemorhedus caudatus	0.15	0.15	0.11	0.11	0.11	0.10	0.07
9	Pantholops hodgsonii	0.15	0.16	0.12	0.13	0.13	0.12	0.13	0.12
10	Bison bison	0.17	0.18	0.17	0.15	0.17	0.17	0.16	0.17	0.18
11	Bos taurus	0.17	0.17	0.17	0.16	0.17	0.16	0.15	0.16	0.17	0.05
12	Bubalus bubalis	0.17	0.16	0.16	0.16	0.15	0.16	0.16	0.16	0.17	0.14	0.15

Note: There were a total of 1,545 positions in the final dataset.

Table S5.

Estimates of evolutionary divergence between sequences of Cytb mtDNA of the used taxa under the family Bovidae.

No.	Species	1	2	3	4	5	6	7	8	9	10	11
1	Ammotragus lervia	0
2	Capra hircus	0.11
3	Budorcas taxicolor	0.13	0.14
4	Ovis aries	0.12	0.12	0.15
5	Pantholops hodgsonii	0.11	0.12	0.13	0.12
6	Capricornis crispus	0.11	0.12	0.14	0.13	0.12
7	Naemorhedus caudatus	0.12	0.13	0.15	0.15	0.14	0.09
8	Antilope cervicapra	0.17	0.18	0.19	0.17	0.16	0.16	0.18
9	Cervus nippon	0.17	0.17	0.19	0.17	0.17	0.17	0.18	0.19
10	Bison bison	0.17	0.17	0.18	0.19	0.18	0.17	0.18	0.19	0.19
11	Bos taurus	0.17	0.17	0.20	0.20	0.18	0.17	0.19	0.19	0.18	0.08
12	Bubalus bubalis	0.16	0.18	0.19	0.18	0.17	0.17	0.18	0.18	0.17	0.14	0.14

Note: There were a total of 1,140 positions in the final dataset.

Table S6.

Estimates of evolutionary divergence between sequences of D-loop mtDNA of the used taxa under the family Bovidae.

No.	Species	1	2	3	4	5	6	7	8	9	10	11
1	Ovis aries	0
2	Cervus nippon	0.23
3	Capra hircus	0.18	0.21
4	Budorcas taxicolor	0.19	0.23	0.17
5	Ammotragus lervia	0.18	0.22	0.15	0.16
6	Capricornis crispus	0.17	0.22	0.16	0.17	0.15
7	Naemorhedus caudatus	0.19	0.23	0.18	0.18	0.16	0.11
8	Bison bison	0.27	0.23	0.24	0.25	0.25	0.25	0.26
9	Bos Taurus	0.26	0.26	0.26	0.27	0.27	0.27	0.28	0.10
10	Bubalus bubalis	0.26	0.24	0.26	0.25	0.25	0.26	0.27	0.19	0.19
11	Antilope cervicapra	0.32	0.30	0.30	0.30	0.29	0.29	0.29	0.30	0.31	0.30
12	Pantholops hodgsonii	0.20	0.26	0.20	0.17	0.19	0.19	0.19	0.23	0.25	0.24	0.25

Note: There were a total of 806 positions in the final dataset.