Sequence analysis of four caprine mitochondria DNA lineages

Abstract

The complete mitochondrial DNA (mtDNA) (16640bp in length) was sequenced from four Chinese goat lineages representing the four major mtDNA haplogroups in goats. A total of 124 single nucleotide polymorphisms (SNPs) were found in encoding regions, and the overall ratio of transitions:transversions was 40:1 revealing a heavy transition/transversion rate in domestic goats. Eighteen non-synonymous sites were found for the total number of SNPs; the sites did not affect the predicted functions of protein for these four goat mtDNA lineages. In the region for coding tRNA and rRNA, SNPs occurred in loops, unstructured single strand and stems that were conformed with the principle of G-U pairing. We came to the conclusion that these substitutions could not change secondary structure of RNAs, and there was no positive selection on goat mitochondrial coding region according to the result of dN/dS (0.0399–0.1529) by comparing the goat with other reported mitochondrial genomes.

Key Words:

• Goat
• Mitochondrial genome
• Phylogenetics.

Introduction

Genetic analysis of mitochondrial (mt) DNA has been an important tool in understanding goat evolution. Attributes such as high copy number per cell, lack of recombination (Olivio et al., 1983), high substitution rate (Brown et al., 1979) and uniparental inheritance (Giles et al., 1980; Gyllensten et al., 1985) have made mitochondrial DNA (mtDNA) a source of information for phylogenetic studies. The circular mitochondrial genome of goat encodes 13 polypeptides, 22 transfer RNA (tRNA) genes and 2 ribosomal RNA (rRNA) genes necessary for transcription and translation of the mitochondrial genome. Among these genes, only NADH dehydrogenase subunit 6 and eight tRNAs are encoded in the L-strand while the others are encoded by the H-strand. With the exception of the ND2, ND3 and ND5 genes, all peptide-coding genes of the goat mtDNA have a methionine start codon (Parma et al., 2003).

Comprehensive studies have demonstrated that goat mtDNA is geographically structured and may be classified into four mtDNA haplogroups observed in China (Chen et al., 2005). Analyses of the previous studies focused on the d-loop that occupies less than 7% of the mtDNA genome (Chen et al., 2005; Joshi et al., 2004; Luikart et al., 2001; Mannen et al., 2001; Naderi et al., 2007) and is a major non-coding region regulating replication and transcription of mtDNAs (Clayton, 1992). Therefore, analysis of the whole mtDNA sequences from various goat lineages would more realistically characterise the evolution of diverse goat lineages. Mitochondrial genome sequences have been widely used in phylogenetic relationship analyses. Currently, the complete mitochondrial genome DNA of human, cattle, pig, horse, donkey, chicken, sheep, goat and yak have been sequenced (Anderson et al., 1981, 1982; Gu et al., 2006; Hiendleder et al., 1998; Lin et al., 1999; Nishibori et al., 2003; Parma et al., 2003; Xu and Arnason, 1994; Xu et al., 1996). Non-synonymous-synonymous substitution rate ratio (dN/dS) is important to evaluate selective pressure based on the protein-coding sequences. When dN/dS is over 1, positive Darwinian selection is inferred; when dN/dS is 1.0, strict neutrality is inferred, and when dN/dS is less than 1, purifying selection is inferred. Regions where dN/dS is over 1 are potential sites of positive selection.

In this study, the complete mitochondrial DNA (mtDNA) sequences of four goat lineages were sequenced. Characteristics of the goat mitochondrial genome were analysed and the structure of mitochondrial RNA and protein were predicted. Comparisons of the mitochondrial sequence of goat with that of other reported species were also made to analyse the molecular evolution of the goat mitochondrial DNA.

Materials and methods

Sampling and DNA extraction

Goat ear tissues were collected from small remote villages in five provinces or regions of China. Samples of the goat ear tissues were stored at −70°C before DNA extraction. Total genomic DNA was extracted from ear tissue by a standard phenol-chloroform extraction method (Sambrook and Russell, 2000). Details of the breeds, regions and sample sizes are given in Table 1.

Table 1 Information of the samples.

Lineage	Breed	Sample size	Resource
A	Langkazi	1	Tibet
B	Yunling Black	1	Yunnan
C	Wuzhumuqin White Cashmere	1	Neimeng
D	Tsaidam Cashmere	1	Qinghai

DNA amplification and sequencing

Based on the first hypervariable region of mitochondrial DNA study, 4 samples (one per lineage) were collected and sequenced according to the MJ network of 145 domestic goats (Wu et al., 2009). Fifteen pairs of primers (Appendix Table 1) overlapping in mitochondrial DNA were designed with Premier primer 5.0 and the known goat mtDNA sequence (Parma et al., 2003). Polymerase chain reaction (PCR) amplifications were carried out in 25 µL reaction mixtures including each primer (1 µL of a 10 µM solution), dNTPs (1 µL of a 2.5 mM solution), 2.5 µL of 10 × buffer and 0.25 µL of 5 U/µL Taq DNA polymerase (Tiangen). The PCR programme consisted of an initial denaturing step at 95°C for 5 min, 33 amplification cycles (94°C for 30 s, 58.5°C for 30 s and 72°C for 30 s) and a final extension at 72°C for 10 min in a PTC200® Programmable Thermal Controller (MJ Research Inc., Waltham, MA, USA). The PCR products were sequenced by Shanghai Bioasia Co. (Shanghai, China). The complete sequences of goat mtDNA lineages were deposited in the GenBank database (Accession n. GU229278-GU229281). Hassanin et al. reveal the presence of numts and multiple sequencing errors, and suggested that the sequences were probably contaminated by a numt fragment of 1881 bp based on a pair of primers located around 12260–13280 bp of the complete goat mtDNA sequences. We carried out further repeated sequencing using a number of short segment primers and the problem was avoided.

Analysis of sequence data

Overlapping contigs were compiled using DNAstar Seqman software (DNASTAR Inc., Madison, WI, USA) to assemble a continuous sequence. Sequence alignments with mitochondrial DNA sequence by Bioedit package (Hall, 1999) were used to annotate goat mitochondrial genes. The influence of amino acid sequence mutation on the structure and function of protein was analysed using online alignment of amino acid sequence. The structure of mitochondrial RNA was predicted using RNAstructure version 4.5 and ViennaRNA version 1.5 software. The neighbour-joining (NJ) tree was constructed using the Mega 3.0 software programme (Kumar et al., 2004), with a Kimura 2-parameter model and a bootstrap (1000 replications) test. dN/dS(ω) values were calculated using PAML according to Yang (1997).

Results

Characteristics of four goat mtDNA lineages

We constructed the complete nucleotide sequence of four goat mtDNA lineages. The goat mtDNA was 16640–16642bp in length, the average contents of the base A,T,C and G were 29.5%, 25.8%, 31.9% and 12.% in the goat mtDNA encoding region, respectively. Average AT contents in the first, second and third codon positions of the protein genes were 55%, 62.1% and 66.1%, respectively, and the third codon position showed obvious AT bias. In total, 124 SNPs were found among the encoding regions of four lineages, 104 of which were distributed in protein-coding sequences. The overall ratio of transitions:transversions (40:1) revealed a heavy transition bias in domestic goats. The d-loop region variations among four lineages have been reported in previous studies (Wu et al., 2009).

Amino acid sequence analysis

Thirteen protein-coding genes with initial codon of methionine (10 ATG and 3 ATA) were found in the goat mtDNA, among which three genes appeared to end in an incomplete stop codon of TNN, where NN was the 5 terminus of the adjacent tRNA gene, which presumably formed a stop codon by post-transcriptional polyadenylation (Anderson et al., 1981). This observation was consistent with characteristics of the goat mitochondrial genome (Parma et al., 2003).

Eighteen non-synonymous sites were found by aligning the amino acid sequence of four goat mtDNA lineages. The influence of amino acid sequence mutation on the structure and function of protein were analysed using online alignment of amino acid sequence (Table 2). The results of polyphen were benign PSIC score: 0.055–1.282, indicating that amino acid mutation of four goat mtDNA lineages did not affect function of protein.

Table 2 Amino acid substitutions among goat lineages.

Gene	Position	SNP of 4 lineage				Amino acid substitutions	Results of polyphen	PSIC score
A	B	C	D
ND1	2943	T	T	C	T	Ile - Thr	Benign	0.089
COX1	6573	A	A	A	G	Ile - Val	Benign	0.233
COX2	7051	A	G	G	A	Thr - Ala	Benign	0.055
COX2	7462	A	A	A	G	Ile - Val	Benign	0.356
ATP6	8054	G	A	A	A	Val - Ile	Benign	0.361
ATP6	8103	G	G	A	G	Ser - Asn	Benign	0.180
COX3	8729	C	T	C	C	Ala - Val	Benign	0.400
COX3	9060	T	T	C	T	Met - Thr	Benign	0.807
ND3	9546	T	T	C	T	Val - Ala	Benign	0.276
ND4L	10000	A	A	A	G	Ile - Val	Benign	0.619
ND5	11807	T	T	T	C	Met - Thr	Benign	0.991
ND5	11870	G	A	G	G	Cys - Tyr	Benign	0.112
ND5	12004	G	G	G	A	Val - Met	Benign	0.851
ND5	12190	A	A	G	A	Ile - Val	Benign	1.282
ND5	12838	G	G	G	A	Ala - Thr	Benign	0.443
ND5	13476	A	C	A	A	Met - Ile	Benign	0.083
ND6	13581	T	G	G	G	Val - Gly	Benign	0.334
CYTB	14161	T	T	C	T	Ile - Thr	Benign	0.694

RNA secondary structure analysis

Twenty-two tRNA genes and 2 rRNA (12S and 16S rRNA) were found in the goat mtDNA (Parma et al., 2003). The tRNA genes were dispersed along the whole mtDNA genome, RNA secondary structure was more conserved than the primary structure and it contained information concerning the phylogenetic relationship. In this study, 22 tRNA and 2 rRNA sequences of four goat mtDNA lineages were aligned by the mega 3 software package. The result revealed variations of sequences in the 6 tRNA and 2 rRNA. The inferred secondary structures of the 6 tRNA genes showed that the variation was A/G transition. Among 12S rRNA and 16S rRNA of four goat mtDNA lineages, there were 4 and 9 differences, respectively. In the region for coding tRNA and rRNA, SNP occurred in loops, unstructured single strand, and stems that conformed with the principle of G-U pairing (Table 3). It can, therefore, be concluded that these substitutions could not change secondary structure of RNAs.

Table 3 Single nucleotide polymorphisms in goat mitochondrial RNA.

Gene	Position lineage	Variation substitution	Base	Domain
tRNA	Ile	3738	D	A-G	Connection of anticodon arm and T arm
	Ala	5083
L strand	D	A-G	Stems (principle of G-U pairing)
Trp	4962	B	A-G	D loop
Gly	9415	D	A-G	Connection of anticodon arm and D arm
His	11567	C	A-G	Connection of anticodon arm and D arm
Leu	11689	B	G-A	D loop
Leu	11700	B	A-G	Connection of anticodon arm and D arm
12SrRNA	165	C	T-C	Loops
203	B	C-T	Loops
211	C	C-T	Loops
292	D	G-A	Loops
16SrRNA	1190	C	T-C	Loops
1367	C	C-T	Loops
1848	BC	T-C	Loops
2090	BC	T-C	Loops
2117	D	T-C	Unstructured single strand
2206	D	T-C	Loops
2415	B	G-A	Loops
2515	B	C-T	Loops
2637	D	A-G	Stems (principle of G-U pairing)

Detection of positive selection

The dN/dS(ω) values were calculated according to Yang (1997). The ratios (dN/dS) of the 4 goat lineages and 9 mammals (sheep, cow, yak, pig, horse, donkey, dog, cat, rat; Accession nos. AF010406, J01394 AY684273 AJ002189 NC_001640 X97337 NC_002008 U20753 DQ673907) in all 13 protein coding genes were calculated (Appendix Table 2). The result indicated the base substitutions among 4 goat lineages and 9 mammals were not affected by Darwinian positive selection (ω=0.0399–0.1529). The same conclusion has also been drawn using the Nei and Gojobori method (Nei and Gojobori, 1986). Phylogenetic reconstruction of 10 mammals was performed using the mtDNA coding region sequences collected in this study and the published sequences of another 9 mammals (Figure 1).

Figure 1 Radio-protective mechanism of melatonin.

Discussion

Estimates of divergence time

Molecular clock analyses are highly sensitive to evolutionary distance between main clades. Furthermore, the choice of an outgroup is important (Rogaev et al., 2006). The Stationary Markov Model (Saccone et al., 1990) was used to measure the sequence divergence, K [substitutions (subs)/site], of the different functional regions for each pair of species. The absolute nucleotide substitution rate, V (subs/site year⁻¹), was calculated as V=K/2T, where T is the divergence time (Pesole et al., 1999). The various functional regions of mtDNA evolved at different rates.

The non-synonymous sites, D-loop central domain, and tRNA and rRNA genes changed much more slowly than the synonymous sites. The nucleotide substitutions among 7 species (goat, sheep, cow, yak, pig, horse and donkey) in 22 tRNA coding sequences and the 12S and 16S rRNA genes were calculated. The tRNA component had a slightly lower average sequence divergence than either the 12S or the 16S rRNA genes (Table 4), which was also observed in other mammalian lineages (Gu et al., 2006; Pesole et al., 1999).

Table 4 Estimates of time of divergence between mtDNA clades present in different populations.

RNA	Mammalian mtDNA specific mutation rate°	Estimated time of divergence^#
		Goat vs sheep	Goat/sheep vs cow/yak	Goat/sheep vs pig	Goat/sheep vs horse/donkey
12rRNA	3.4±1.2	9.53	17.63	27.69	33.57
16rRNA	4.9±1.4	9.86	13.25	20.55	21.52
tRNA	3.4±1.1	8.68	12.83	26.85	29.08

° Estimated rate of nucleotide substitution (substitutions/site/year×10⁻⁹) specific for each functional component of the mammalian mtDNA genome (Pesole et al., 1999).

# The time since divergence (years×10⁶) was estimated on the basis of data from each functional component using the pairwise sequence divergence values appearing in Tables 3–11 in Wu, 2008.

The time since divergence between the goat and other mammals was based on previously observed substitution rates specific for each mammalian mtDNA functional component (Pesole et al., 1999). The time since divergence was found to be 9.53–9.86 million years between goat and sheep, 12.83–17.63 million years between goat/sheep and cow/yak, 20.55–27.69 million years between goat/sheep and pig, and 21.52–33.57 million years between goat/sheep and horse/donkey (Table 4). Gu et al. (2006) estimated that divergence between goat/sheep and cow/yak occurred at 13.14–27.99 million years before the present day. This is consistent with the data of 12.83–17.63 million years in the study.

Complete mitochondrial DNA and molecular evolution

The study of molecular evolution is made more difficult by the fact that only a few full nuclear genome sequences have been completed in eukaryotes and prokaryotes. Mitochondrial genomes offer a broad range of characters to study phylogenetic relationships of animal taxa. Besides nucleotide and amino acid sequences, tRNA secondary structures (Macey et al., 2000), deviations from the universal genetic code (Castresana et al., 1998; Telford et al., 2000), as well as changes in the mitochondrial gene order (Boore et al., 1995; Boore et al., 1998) have been successfully used as characters in phylogenetic inference. In particular, the changes in gene order prove to be extremely reliable phylogenetic characteristics because the probability that homoplastic translocations occur in closely related taxa is very low (Pesole et al., 1999).

There is increasing archaeological evidence suffesting that goats, in the form of their wild progenitor the bezoar (Capra aegagrus), were the first wild herbivores to be domesticated. The studies in literature suggest that this happened 10,000 years ago at the dawn of the Neolithic period, in the Middle East region known as the Fertile Crescent (Legge, 1996; Zeder and Hesse, 2000). Parma et al. (2003) reported that the whole mitochondrial genome sequence of goat was completed by using the available goat sequences aligned with the cattle mtDNA. The resulting gaps were then amplified by PCR using goat specific primers; the complete sequence did not belong to any lineage. In the study, the complete nucleotide sequence of four goat mtDNA lineages was sequenced. The 4 sequences represented four lineages. The goat mtDNA sequences submitted to GenBank had different lengths (16,640; 16,641 and 16,642) as a result of base insert/delete. The complete sequence of four goat mtDNA lineages will be useful for further study of goat evolution and genetics.

Conclusions

The goat mtDNA was 16640–16642bp in length, with the third codon position showing obvious AT bias. Eighteen non-synonymous sites were found by aligning amino acid sequence of four goat mtDNA lineages, and the non-synonymous sites did not affect the function of protein. In the region for coding tRNA and rRNA, SNP occurred in loops, unstructured single strand, and stems, which conformed with the principle of G-U pairing. It was, therefore, concluded that these substitutions could not change the secondary structure of RNAs. The ratios (dN/dS) of the four goat lineages and 9 mammals in all 13 protein coding genes were calculated and the results indicated that the base substitutions among four goat lineages and 9 mammals were not affected by Darwinian positive selection (ω=0.0399–0.1529).

Acknowledgements:

the authors are grateful to Dr. Lingjuan Wang Li for his fundamental review of this manuscript.

This study was supported in part by the National High Technology Research and Development Programme (863) of China (grant n. 2006AA10Z198), the National Key Technology Research and Development Program of the Ministry of Science and Technology of China (grant n. 2006BAD13B08) and the Special Funds of Basic Research Funds of the Ministry of Science and Technology of China (2006FY1107 00).

APPENDIX

Appendix Table 1 Primer sequence, corresponding polymerase chain reaction product size and position.

Code	Primer sequence	Product size	Position	Temperature
G1	5‘ CGCCCTCCAAATCAATAA 3’	639	144–782	53
5‘ CCCATTTCTTCCCATTCC 3’
G2	5‘ CCCATTTCTTCCCATTCC 3’	621	765–1391	52
5‘ TTGCCACATAGATGAGTTG 3’
G3	5‘ CAACTCATCTATGTGGCA 3’	991	1373–2363	53.2
5‘ TAGGGTAACTTGTTCCGT 3’
G4	5‘ AGTTTTGGTTGGGGTGAC 3’	1405	2237–3641	56
5‘ ATAGGGCCAGTGTTAGGG 3’
G5	5‘ TACCCCCGATTCCGTTAC 3’	1421	3569–4989	56
5‘ TTAGGGCTTTGAAGGCTC 3’
G6	5‘ ACTCACACCAATCCTCTCAAT 3’	1359	4919–6277	56
5‘ CCAGTTGGGATAGCGATAAT 3’
G7	5‘ GATACACGGGCTTACTTCAC 3’	1227	6225–7451	56.8
5‘ ATGGGTAGTACAACTCGGTT 3’
G8	5‘ AAACCTGGAGAACTACGA 3’	1427	7399–8825	55
5‘ GGAGTATGGTGACCTTGG 3’
G9	5‘ ATACCAATGATGACGAGATGTG3’	1155	8769–9923	56
5‘ AGAGATACTGCGAATGCTGTTA 3’
G10	5‘ AACAGCATTCGCAGTATCT 3’	1182	9903–11084	55
5‘ GGGTGTCTGGATAAGGATG 3’
G11	5‘ CCTGAAATCACTCATCGC 3’	1247	11009–12255	55
5‘ TTGTAGGGCTGCTGTGTT 3’
G12	5‘ GGTGGTATGGACGAGCAGA 3’	1148	12215–13362	58
5‘ GTAGGGCATTAAGCGGTGT 3’
G13	5‘ TTGGAAGCCTTTTCACAG 3’	1289	13121–14409	54
5‘ CTCCGTTTGCGTGTATGT 3’
CB	5‘ AATGATATGAAAAACCATC 3’	1344	14106–15449	50
5‘ CCACTATTAACCACATCTA 3’
DL	5‘ AGACTCAAGGAAGAAGCCA 3’	1438	15359–16640	55
5‘ GATTTGGAGGGCGTTACT 3’	1–156