Implications
This work adds a substantial amount of data to the existing records of bovine mitochondrial diversity, which is important for developing DNA-based traceability systems. Based on these data, we also address the issue of mitochondrial recombination in cattle. We present the evidence that support the possibility that this diversity arose from recombinations: haplotypes exceed mutations in number; and the genomic encoding for enzymes such as mitofusins and DNA ligase 3 that mediate mtDNA recombination. Whether mtDNA recombination takes place at all remains an open and fundamental question for basic research of the mammalian mitochondrion.
Introduction
The first reported complete sequence of bovine mitochondrial (mt) DNA (GenBank Accession no. V00654; Anderson et al., Reference Anderson, Debruijn, Coulson, Eperon, Sanger and Young1982) was 16 338 bp long and is often used, as in this work, as the reference genome for numbering the positions of elements in the bovine mitochondrial genome. MtDNA is maternally inherited and can be employed in identity and maternity control using a stretch of 910 bp in the control region of the mitochondrial genome. This highly polymorphic region is also called the D-loop as it contains the origin of replication, which forms a displacement-loop structure when a third strand of DNA is taken up by the double-stranded DNA during replication of the circular mitochondrial chromosome. Genetic markers, including variation in the D-loop region, provide the most reliable way of DNA-based traceability linking the source of the carcass to meat products (reviewed by Dalvit et al., Reference Dalvit, De Marchi and Cassandro2007). Specific advantages of using D-loop mtDNA for cattle identification are: (a) its large allelic variation, which has not been affected by the loss of autosomal genetic variability in cattle since the introduction of artificial insemination (AI) (Taberlet et al., Reference Taberlet, Valentini, Rezaei, Naderi, Pompanon, Negrini and Ajmone-Marsan2008); (b) the copy number of bovine mtDNA molecules, ranging from 220 to 1720 per cell, and its applicability in severely decomposed samples in which nuclear DNA has already degraded (Xu et al., Reference Xu, Reuter, Xu, Alexander, Gilroyed, Jin, Stanford, Larney and McAllister2009); (c) typing is efficiently performed by the widely used method of PCR-direct sequencing; (d) similar to microsatellites, it is a multiallelic marker, but its direct sequencing is less prone to error than microsatellite genotyping and unlike single nucleotide polymorphisms (SNPs), it can quickly detect the samples containing DNA from multiple individuals by revealing more than two alleles per marker (Karniol et al., Reference Karniol, Shirak, Baruch, Singrun, Tal, Cahana, Kam, Skalski, Brem, Weller, Ron and Seroussi2009). The high copy number of the mitochondrial genome may become a disadvantage in direct sequencing since this method produces ambiguous results following the indel position in cases of length heteroplasmy, which has been reported to be widespread in cows (Wu et al., Reference Wu, Smith, Freeman, Beitz, McDaniel and Lindberg2000). Moreover, variation of heteroplasmy levels within organs and tissues of the same individual (Wu et al., Reference Wu, Smith, Freeman, Beitz, McDaniel and Lindberg2000) can further complicate the use of the mtDNA D-loop code for individual identification. The distribution of mtDNA haplotypes and the possible sub-population structures within the relevant populations are important for identity control (Bar et al., Reference Bar, Brinkmann, Budowle, Carracedo, Gill, Holland, Lincoln, Mayr, Morling, Olaisen, Schneider, Tully and Wilson2000). To gain insight into this distribution, we compared 1515 bovine mtDNA D-loop sequences available from GenBank to 219 sequences derived from Israeli domestic and introduced cattle. Most of the sample consisted of Holsteins, including international AI sires. The local breed originated from an initial import of Dutch and Friesian cattle in the early 1920s, which was crossbred with local Baladi and Syrian Shami cattle (http://www.icba-israel.com). Beef cattle samples included Hereford, Simmental, Charolais and Limousin breeds, which were similarly crossbred with local breeds. Aberdeen Angus and Brahman were recently introduced from Australia. Based on the D-loop sequence data obtained and genomic evidence for the coding for the molecular machinery for recombination of mtDNA in the bovine genome, we also address the issue of mitochondrial recombination, which remains an open question for basic research of the mammalian mitochondrion.
Material and methods
Animals
DNA from the 219 individuals was extracted from either semen or ear tissue as previously described (Karniol et al., Reference Karniol, Shirak, Baruch, Singrun, Tal, Cahana, Kam, Skalski, Brem, Weller, Ron and Seroussi2009). This sample included sire semen of Israeli (n = 68), United States (n = 13), Dutch (n = 12), French (n = 5), German (n = 4), Danish (n = 2), Italian (n = 1), British (n = 1) and Canadian (n = 1) origins. The actual identities of the sampled sires are detailed in their GenBank records and can also be determined through their four-digit identity number from The Israeli Cattle Breeders Association (http://www.icba.org.il/cgi-bin/bulls/en/bl_main.htm). The rest of the sample (n = 112) was obtained from ear tissue of the following breeds: Holstein (n = 38), Angus (n = 23), Simmental (n = 21), Charolais (n = 9), Limousin (n = 3), Hereford (n = 9) and Brahman (n = 9).
Identification of polymorphism
Template DNA was PCR-amplified using the forward primer: 5′-GCCCATACACAGACCACAGA-3′ and reverse primer: 5′-GGGGCCTGCGTTTATATATTG-3′; and BioMix™ Red kit (Bioline Ltd, London, UK) according to the manufacturer’s instructions and the following conditions: 30 cycles for 40 s at 92°C, 60 s at 63°C and 60 s at 72°C. PCR products were separated on agarose gels, excised, purified with AccuPrep® Gel Purification Kit (BioNeer Corp., Seol, Korea) and then bidirectionally sequenced using an ABI3730 sequencer and the above primers. Sequences were assembled and compared using the GAP4 program (Staden et al., Reference Staden, Beal and Bonfield2000).
Results
Cost effectiveness is important for cattle traceability and the entire mtDNA D-loop (910 bp) is too large to be sequence-genotyped in a single economically feasible reaction. Therefore, we limited our analysis to the effective range of a single-sequence read by designing primers that were capable of PCR-amplifying most of the variable regions (791 bp), which have been proven useful for traceability (Figure 1b). To ensure efficient amplification of the target DNAs, these designed primers were placed in conserved regions according to the sequence information that was available when this project was first initiated (2005). Most (25/41) nucleotide positions corresponding to these primers remained perfectly conserved across the 1337 sequence submissions that had become available when this report was prepared (Figure 1b and Supplementary Table 1, positions 15921-15940 and 334-355, see http://cowry.agri.huji.ac.il/SupTable1.xls).
(Colour online) Variation in Bos taurus type D-loop mtDNA. (a) Delineation of the D-loop region corresponding to positions 15921-16338 and 1-355 (Anderson et al., Reference Anderson, Debruijn, Coulson, Eperon, Sanger and Young1982) with annotation of the functionally important domains following Steinborn et al., Reference Steinborn, Muller and Brem1998, including the origin of replication (red), stem-loop forming regions (pink), sequence elements (green), conserved sequence boxes (brown), protein-binding domains (blue), transcriptional start sites of L-strand (LSP) and H-strand (HSP) promoters (light purple) and LSP-enhancing domain (grey). (b) Distribution of point mutations observed across the control region in alignment with the delineation in (a). The columns, which indicate the total number of occurrences in each site, correspond to the logarithmic scale while the distribution within each column of mutations detected in this study (red series) and previously reported (blue series) is linear. (c) Reverse (upper) and forward (lower) sequence chromatograms representing ambiguous trace derived from Holstein international AI sire McCloe-Pond Trent, USA, 17226843 (IL no. 5247). (d) Reverse sequence chromatograms of common (upper) and rare (lower) haplotypes. Arrows point to the two novel mutations that coincide in this sire sequence (GenBank accession no. FN557401).
Eight novel mtDNA D-loop mutations
Using these primers, we sequence-analyzed the 219 individuals described above. The sequence data of 209 individuals were submitted to EMBL (GenBank accession numbers FN556936-FN556985; FN557366-FN557440; FN562583-FN562665; FN563411). We excluded 10 individuals that gave pronounced ambiguous traces, detected mostly following the C-stretch in positions 216 to 221 (Figure 1b) and evident from both the forward and reverse sequence chromatograms (Figure 1c). The presence of indels in these individuals was confirmed using ShiftDetector (Seroussi et al., Reference Seroussi, Ron and Kedra2002). When assembled, the rest of the data consisted of 143 different alleles (Supplementary Table 2, http://cowry.agri.huji.ac.il/SupTable2.xls), most of them derived from a single (n = 85) or two (n = 19) individuals. The allele with the highest frequency was encountered 16 times.
To compare the bovine mtDNA D-loop sequences available from GenBank to the sequences we generated, we Blast-searched the ‘non-redundant’ database (http://blast.ncbi.nlm.nih.gov/Blast.cgi) using the 791-bp consensus sequence and the Entrez query limits of organism Bos taurus and sequence length >700 bp. The latter boundary aimed to limit our search to high-quality sequences. Records with ambiguous nucleotides were also excluded. Following downloading, the detected mtDNA D-loop sequences (n = 1515) were annotated using their breed and accession-number records and assembled into contigs, allowing allele identification and comparison in a unified GAP4 database of 1724 records (Supplementary Table 2). This analysis, using the normal assembly default with a maximum percent mismatch of five, divided the sequences into three contigs: the first corresponded to the B. taurus haplotypes and contained most of the sequences (n = 1532). The second (190 sequences) contained the Bos indicus haplotypes. These included our sequences, which were derived from two Brahman and 13 Angus individuals, whose breeds were introduced from Australia. As evident from this contig, there were no previous reports of Aberdeen Angus carrying the B. indicus mtDNA genome. The third contig consisted of sequences from only two individuals of the Qinghai yellow cattle (AB177769, AB177776), which matched Bos grunniens and suggested the hybridization of this species with B. taurus in the Qinghai province of China (Jia et al., Reference Jia, Haisheng, Eiji and Shinobu2004; Qi et al., Reference Qi, Jianlin, Wang, Rege and Hanotte2010). A second analysis using a maximum percent mismatch of 0 divided the sequences into 925 contigs that represented different haplotype alleles, 858 and 67 of the B. taurus and B. indicus types, respectively (Supplementary Table 2). We further examined the mutations associated with the B. taurus-type haplotypes (Figure 1b; Supplementary Table 1). Plotting the mutated site along the analyzed mtDNA D-loop region (Figure 1b) indicated the influence of functional elements drawn above this region’s delineation (Figure 1a) on the mutation rate within it. While the stem-loop-forming domain of the 5′-cloverleaf around the replication origin was nearly conserved, most mutations clustered in variable regions. Eight novel mutations were observed (Table 1), including a case in which two of these mutations coincided in the same individual (Figure 1d). Most mutations detected in this study (red series, Figure 1b) occurred in previously reported sites (blue series, Figure 1b). In total, just 302 variable sites formed the 858 haplotypes detected. However, 55 sites displayed a complex level of variation (Supplementary Table 1), including one position (332) that reached the maximal five levels of variation (the four possible base nucleotides and an indel status). Taking into account that each level represents an independent mutation event, a total of 399 mutations were traced, which could potentially explain independent formation of less than half (47%) of the haplotypes encountered for the B. taurus type: most of the haplotypes were derived from different combinations of these mutations.
Discussion
The extent of the effect of recombination on animal mtDNA variability remains an open question (see for review Rokas et al., Reference Rokas, Ladoukakis and Zouros2003). As in other mammals, mtDNA recombination is thought to be absent in cattle, yet mitofusins (Entrez GeneIDs: 515180; 534574), which catalyze mitochondrial fusion, are present in the bovine genome, as are the enzymes for mtDNA recombination that are part of the mtDNA repair and replication machinery, such as bovine DNA ligase III (LIG3, Entrez GeneID: 514719). The mitochondrial nature of this latter enzyme is evident from its first 87 amino acids, which are analogous to the double-methionine structure (86 residues, 65% identity and 77% similarity) at the amino terminus of human LIG3. The second methionine has been designated as the first residue for the nuclear LIG3 protein, while the first one is followed by a putative mitochondrial-targeting sequence and is exploited for the mitochondrial isoform (Lakshmipathy and Campbell, Reference Lakshmipathy and Campbell1999). Another line of evidence for a possible effect of mtDNA recombination is the demonstrated decline in linkage disequilibrium as a function of distance between polymorphic sites in some ruminant populations (Slate and Phua, Reference Slate and Phua2003). Recently, Neiman and Taylor, Reference Neiman and Taylor2009 reviewed and analyzed the accumulation rate of deleterious mutations in mtDNA concluding that a complete lack of recombination will ultimately lead to extinction when compensatory mutation is rare, and that undetectable levels of recombination may effectively counter this mutation accumulation.
The basis for the bovine mtDNA diversity can be further assessed by analyzing the 1724 sequences available from GenBank and derived from the variable D-loop region, including the 209 submissions that were deposited as part of this work. In the analyzed D-loop region (791 bp), we encountered 302 variable sites and 399 variations, which in different combinations produced 858 different haplotypes. Besides recombination, a likely explanation for the haplotype number exceeding the number of traced mutation events is that the distribution of mutations is not random along the D-loop. Thus, independent recurrences and reversals of the same mutations in some variable sites may produce an abundance of haplotypes, using a relatively limited number of mutation sites. As shown in Figure 1, this hypothesis is supported by the data as most mutations clustered in variable regions, while the existence of the functional element delineated above (Figure 1a) limited the observed variation in the corresponding regions (Figure 1b). The mutational hotspot explanation and the impact of functional elements on mtDNA D-loop patterns in cattle have been thoroughly discussed by (Slate and Phua, Reference Slate and Phua2003), who concluded that they should be considered a serious alternative to recombination.
Regardless of their cause, the diverse haplotypes of the analyzed 791-bp mtDNA D-loop region present an effective alternative for identity control, due to the theoretical number of combinations (5302) generated by the five potential levels of mutation in the 302 variable sites. Indeed, a considerable proportion (41%) of the individuals we tested displayed a unique novel haplotype. However, some haplotypes were much more frequent than others: the most frequent haplotype among the sequences generated in this study (Supplementary Table 1) was encountered 16 times (7%) in four different breeds and only once in previous GenBank submissions, underscoring the importance of investigating different allele frequencies and their distribution within the relevant population (Bar et al., Reference Bar, Brinkmann, Budowle, Carracedo, Gill, Holland, Lincoln, Mayr, Morling, Olaisen, Schneider, Tully and Wilson2000). The commonly observed haplotype may thus originate from a local Middle-Eastern breed.
The exemplified degree of polymorphism in the D-loop suggests that direct sequencing may be an efficient method of typing D-loop mtDNA. Data analysis is amenable to automation using algorithms for sequence comparisons, but might be complicated by ambiguous sequence results. We removed 10 individuals (5%) from our submissions, which had pronounced ambiguous traces (Figure 1C). Base ambiguity can be explained either by heteroplasmy or, more commonly, length heteroplasmy, which has been suggested to occur naturally due to DNA-strand slippage during replication of the homopolymer regions (Laipis et al., Reference Laipis, Vandewalle and Hauswirth1988; Wu et al., Reference Wu, Smith, Freeman, Beitz, McDaniel and Lindberg2000), or by experimental artifacts. The latter can result from similar in vitro slippage during PCR or from sample contamination. The effect of contamination would be expected to be most pronounced when the target DNA has variations within the regions to which the PCR primers anneal, giving the contaminating sequence a potential replication advantage over the original template. It should be noted that the source of such contamination could also potentially reside in the autosomal DNA, for example, a stretch with 78% identity to the 791-bp D-loop segment was present in build 4.0 on BTA2 positions 2093660-2094385.
In heteroplasmic lineages, unequal partitioning of mitochondrial genotypes among siblings and tissues (Laipis et al., Reference Laipis, Vandewalle and Hauswirth1988) may further complicate identity and maternity control. Thus, although mtDNA typing is a valuable tool for traceability, it should be supplemented with other methods, such as SNP genotyping, which are suitable for automation (Karniol et al., Reference Karniol, Shirak, Baruch, Singrun, Tal, Cahana, Kam, Skalski, Brem, Weller, Ron and Seroussi2009).
Acknowledgements
This work was supported by the Israel milk marketing board. Contribution from the Agricultural Research Organization, Institute of Animal Science, Bet Dagan, Israel, no. 554/09.
