A Comparison of Genetic Diversity of COX-III Gene in Lowland Chickens and Tibetan Chickens

To obtain a full understanding of the genetic diversity of the cytochrome oxidase III gene (COX-III) and its association with high altitude adaptation in Tibetan chickens, we sequenced COX-III in 12 chicken populations (155 Tibetan chickens and 145 other domestic chickens). We identified a total of 11 single nucleotide polymorphisms (SNPs) and 12 haplotypes (Ha1–Ha12). Low genetic diversity (haplotype diversity = 0.531 ± 0.087, nucleotide diversity = 0.00125) was detected for COX-III, and haplotype diversity of Tibetan chicken populations (0.750 ± 0.018) was markedly higher than lowland chicken populations (0.570 ± 0.028). Obvious genetic differentiation (nucleotide divergence = 0.092~0.339) and conspicuous gene communication (gene flow = 0.33~32.22) among 12 populations suggested that Tianfu black-bone fowl (white feather) was possibly introduced from Tibetan chicken. SNP m.10587 T>C affects the specific functions of the COX enzyme. Haplotype Ha3 was found in Tibetan chickens, and SNP m.10115G>A caused an amino acid substitution (Val62Ile) associated with phospholipid binding, while mutations m.10017C>A and m.10555G>A and the previously reported SNP m.10065T>C reduced the hydropathy index to some extent. Together, this indicates that the mitochondrial membrane is more hydrophobic in Tibetan chickens.


Introduction
Domestic chickens fulfill various roles ranging from food and entertainment to religion and ornamentation [1]. Tibetan chicken is a widely distributed aboriginal chicken breed found at altitudes ranging from 2200 to 4100 m; it has adapted well to high altitudes after over the years of living on the plateau [2]. Generally, long-term exposure to hypoxia in animals reduces metabolic activity, retards development, and increases embryo mortality [3]. However, Tibetan chicken has developed an adaptive mechanism to hypoxia, demonstrated by its increased hatchability and survival rate compared with lowland chicken breeds in high altitude areas of Tibet. With its low weight, small size, and strong chest and legs, the appearance and behaviour of the Tibetan chicken resemble those of the Cochin-Chinese red jungle fowl (Gallus gallus gallus), making it particularly good at flying and foraging on the plateau alpine region [2][3][4][5]. The breed is also important to the resources used to expand the industry in cold areas of high altitude in China.
Avian species living at high altitudes are characterized by the high oxygen affinity of their haemoglobin [6,7]. Tibetan birds can improve their physiological performance by enhancing their oxygen transport capacity, which has yielded important insights into the genetic basis of adaptation involving haemoglobin as an oxygen carrier [8,9]. Taking into account the importance of using oxygen more efficiently under hypoxic conditions, information mining for the cellular respiratory chain is a valuable way of understanding the hypoxia response and adaptation mechanisms.
Mitochondrial DNA (mtDNA) sequences are widely used in molecular evolutionary studies. These sequences are useful for estimating times of species and population divergences, comparisons of relative rates of evolution, and phylogenetic inferences within and between vertebrate species [10]. The complete sequence of the chicken mtDNA is 16,775 base pairs 2 BioMed Research International and contains 13 protein coding genes, two rRNA genes, and 22 tRNA genes [11]. Cytochrome c oxidase (COX), the terminal enzyme of the mitochondrial respiratory chain, contains 14 protein subunits in mammals, of which three (COX-I, COX-II, and COX-III) are synthesized in the mitochondria [12]. Mitochondrial cytochrome c oxidase subunit I (COX-I) is the central catalytic subunit of cytochrome c oxidase (complex IV), and COX-I gene is used as a standard marker for DNA bar coding to enable species identification in animals [3,13]. The COX-III protein is an important element in regulating the efficiency of proton translocation in cytochrome oxidase over several turnovers [14]. The COX-III gene of bar-headed geese contained a nonsynonymous substitution (Trp-116→Arg) that resulted in a major functional change of amino acid class. This mutation was predicted by structural modeling to alter the interaction between COX-III and COX-I, which contributed to adaptation in mitochondrial enzyme kinetics and O2 transport capacity and may finally contribute to the exceptional ability of bar-headed geese to fly at extreme heights [15].
COX has been investigated in several biological studies [3,6,16,17], but studies of COX-III have rarely been reported in chickens. However, COX-III is an important component of the respiratory chain and is very conserved among species. Adaptive changes in COX activity can alter the ATP supply derived from oxidative phosphorylation during hypoxia [9]. The quaternary structure formed by different protein subunits is stabilized mainly through hydrophobic interactions in spite of hydrogen bonding and the van der Waals force is also important. We hypothesized that the stability of the COX holoenzyme three-dimensional structure would increase in line with increases in COX-III protein hydrophobicity. Therefore, in the present study, based on the assumption that cytochrome c oxidase activity is more stable in Tibetan chickens, we analyzed COX-III SNPs in 12 chicken populations (five lowland and seven highland populations) to better understand the COX-III genetic diversity and to determine the contribution of specific SNPs to high altitude adaptations in Tibetan chickens.

Materials and Methods
In all experimental populations, 5 populations (Muchuan, Emei, Jiuyuan, Black-Tianfu, and White-Tianfu) belonged to lowland chickens; the other 7 populations (Haiyan, Doilungdêqên, Ganzi, Nyingchi, Diqing, Shannan, and Shigatse) which belonged to Tibetan chickens were collected (Table 1). Blood samples were collected from the wing vein. No bird was slaughtered or unexpectedly injured during sampling. The protocol was approved by the Committee on the Care and Use of Laboratory Animals of the State-Level Animal Experimental Teaching Demonstration Center of Sichuan Agricultural University (Approval ID: Decree number S20160906).
2.1. DNA Extraction, Amplification, and Sequencing. We extracted mtDNA by salt extraction method [18]. PCR used the known primer pairs F9797: 5 -ACCAATAATACCATC-AATCTCC-3 and R10830: 5 CGCTTAGTAGAAAGG-ATAGTGAG-3 [19,20]. PCR amplification was performed in a 50 l volume with 100-150 ng of genomic DNA, 25 mM MgCl2, 2.5 mM of dNTP mixture, 2 mM each primer, 5 l of 10x buffer, and 1.25 U LA Taq polymerase (Takara, Dalian, China) under the following conditions: denaturation at 94 ∘ C for 5 min, then 35 cycles of 94 ∘ C for 30 s, 55 ∘ C for 30 s, and 72 ∘ C for 60 s, followed by a final extension at 72 ∘ C for 7 min [19]. PCR products were verified on 1.5% agarose gels, and specific bands were purified using the TIANgel Midi Purification Kit (Tiangen Biotech, Beijing, China). Purified PCR products were sequenced in both directions using the Big Dye Terminator v. 3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA) on the ABI Prism 3100 DNA sequencer (Applied Biosystems) according to the manufacturer's instructions.

Sequence Data Analysis and Statistical Analysis.
Raw sequences were aligned and edited by DNAstar software (DNAstar Inc. Madison, WI, USA). We exported all sequences as an aligned FASTA file. Sequence variations were identified using MEGA 6.0 software [21]. Standard population genetics statistics, including haplotypes and number of haplotypes, haplotype diversity within each group (Hd), nucleotide diversity (Pi), nucleotide divergence (Dxy), net genetic distance (Da), coefficient of differentiation (Gst), and gene flow (Nm), Tajima's value neutral test, were defined using DnaSP V5 software [22], whereas median joining network analysis was performed using program network 4.611 (http://www.fluxus-engineering.com/sharenet .htm). The complete mitochondrial genome sequence of the red jungle fowl was used as the reference sequence (GenBank accession number: NC 001323).
Analysis of molecular variance (AMOVA) was estimated using Arlequin 3.0 software [23]. Bayesian inference was performed as previously described [24]. The optimal model for each data set was estimated by the program Modeltest 3.7 [25]. The program BEAUti v1.5.3 (distributed with BEAST) was used to create the input file to run in BEAST (http://beast.bio.ed.ac.uk/). Samples from the posterior were summarized on the maximum clade credibility tree using the program TreeAnnotator v1.4.8 (distributed with BEAST) and visualized using the program FigTree v1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/). Statistical differences in COX-III haplotype frequencies between Tibetan chickens and lowland chickens were analyzed using Fisher's exact test; odds ratios (OR) and 95% confidence intervals (95% CIs) were also calculated. value < 0.05 was taken into account as statistical significance. The bioinformatics platform MitoTool (http://www.mitotool.org/) was used to analyze haplotype distribution frequencies between Tibetan chickens and lowland chickens [26]. < 0.05 were taken to be statistically significant.

Nucleotide Diversity of COX-III.
The length of 300 COX-III sequences was truncated into 784 bp (GenBank accession numbers: NC 001323; no insertion/deletions were detected). We calculated the overall base composition of COX-III from the 12 chicken populations using MEGA 6.0. Cytosine (C) was shown to be the rarest nucleotide (16%) with guanine (G) to the most common (32%). The A + T% was around half of COX-III (52%).
A total of 11 SNPs (including six singleton sites and five parsimony-informative sites (with no insertions/deletions)) accounted for 1.403% of the total 784 bp COX-III sequence from 300 individuals. Table 2 summarized the number of haplotypes, Hd, Pi, and other information within each group. The number of variable sites in each population varies from 1 in Diqing and Shigatse to 7 in Haiyan, and the sequences from Doilungdêqên, Ganzi, Nyingchi, Shannan, were 5, 5, 3, and 2, respectively. The highest haplotype diversity was found in Haiyan and Doilungdêqên. While the average number of nucleotide differences in Emei from lowland chickens was greater than others. Taken together, a total of 8 haplotypes were identified in lowland chickens, and the overall haplotype diversity, nucleotide diversity, and average nucleotide differences were 0.570 ± 0.028, 0.00155, and 1.216, respectively. In Tibetan chickens, a total of 9 haplotypes were identified, and the overall haplotype diversity, nucleotide diversity, and average nucleotide differences were 0.750 ± 0.018, 0.0016, and 1.250, respectively. The result showed that the genetic diversity of the Tibetan chicken was noticeably higher than that of lowland chickens. After Tajima's value neutral test, all values were greater than 0.1; therefore, the 12 populations belong to neutral mutations.

Nucleotide Divergence and Net Genetic Distance among
Populations. Nucleotide polymorphism among populations can be represented by nucleotide divergence (Dxy) and net genetic distance in nucleotides (Da). Within the 12 chicken populations, the average Dxy was 0.062% (range 0.092%-0.339%) and the average Da was 0.196% (range 0%-0.216%) ( Table 3). The largest Da (0.216%) was found between Diqing and Jiuyuan black chickens (between Tibetan chickens and lowland chickens), and the smallest Da (0%) was observed between Haiyan and Doilungdêqên chickens; the smallest Dxy (0.092%) was observed between Black-Tianfu and Muchuan black chickens (within lowland chickens) and the largest Dxy (0.339%) was found between Diqing and Emei chickens (between Tibetan chickens and lowland chickens).
These findings are conclusive of definite genetic differentiation between different groups of chickens, with the strongest differentiation seen between Tibetan chickens and lowland chickens. Net genetic distances in nucleotides (Da) are substantially consistent with the outcome of nucleotide differences (Dxy).  (Table 4). This indicates that differences in varieties caused genetic differentiation and that the introduction of varieties of different regions led to gene exchange.

Analysis of Molecular
Variance. AMOVA showed that the percentage of variation within populations (70.43%) was greater than that between populations (29.57%). Fst value was 0.2957 ( * * < 0.01) which implied that the genetic divergence within populations was significant. The results indicate that the twelve geographic populations do not produce largely genetic differentiation, while the genetic diversity in COX-III gene mainly comes from within the populations (Table 5).  Table 6 shows the observed allele frequencies in each polymorphic site between Tibetan chicken and lowland chicken breeds. After using Pearson chi-square test, we found that three SNPs compared with lowland chickens in COX-III gene (m.10115G>A, m.10270G>A, and 10587 T>C) were significantly different with the allele frequency of 9.7%, 91.0%, and 86.5% ( * < 0.05) in Tibetan chicken, respectively. Similarly, in lowland chicken, the SNP m.      Note: relative content in parentheses meant the number of birds in corresponding allele in the specific polymorphic site; * < 0.05 meant significant difference between lowland chicken and Tibetan chicken breeds; a TC was the abbreviations for Tibetan chicken; b LC was the abbreviations for lowland chicken.

Median Joining Network of Haplotypes and Phylogenetic
Analysis. We identified 12 haplotypes (Ha1-Ha12) in 300 chickens from the 12 different populations ( Table 7). The median joining network was constructed using the 12 haplotypes. Three clusters (A, B, and C) were clearly defined from the network with substantial mutation distances visible between the clusters (Figure 1). Ha1, Ha3, Ha8 We also used the Meleagris gallopavo as an outside group to construct a phylogenetic tree using the Bayesian method.
The correlation of the 12 haplotypes is shown in Figure 2, and three clusters can also be seen in the Bayesian tree.

Prediction and Analysis of Secondary Structure Changes in the COX-III Protein.
Protein sequencing showed that the Tibetan chicken-specific nonsynonymous COX-III variants   1  8  15  22  29  36  43  50  57  64  71  78  85  92  99  106  113  120  127  134  141  148  155  162  169  176  183  190  197  204  211  218  225  232  239  246   By analyzing population haplotype diversity and the structure of chicken breeds in southwestern China based on COX-III sequences, we found that the Tibetan chicken has a higher level of haplotype diversity (0.750 ± 0.018) than lowland chickens (0.570 ± 0.028). For the results observed, we think that artificial selection leads to reduced nucleotide diversity in lowland chickens. Although 11 SNP sites were identified in COX-III from 12 geographic populations, the average genetic distance was only 0.196%, revealing a low level of genomic polymorphisms in all populations. This low level of genetic diversity indicates that this gene is functionally important and hence has an evolutionary constraint. The 12 populations belong to neutral mutations ( > 0.1). This suggests that these populations are stable and have not undergone a population expansion in the past few years.
The analysis of the coefficient of differentiation and gene flow between populations suggests that some chicken populations have undergone genetic differentiation at places of equal altitude. This may reflect the introduction of varieties of different regions. For instance, the Tianfu white-bone fowl from a lower altitude has a close population genetic relationship with the Tibetan chicken (Gst = 0.01, Nm = 25.21; Gst = 0.05, Nm = 4.9), indicating that the Tianfu whitebone fowl was introduced from high altitude area. Moreover, Tibetan chickens from Haiyan and Doilungdêqên regions were most similar (Da = 0%, Dxy = 0.177%) suggesting that the genetic distance of the two regions is very close, and gene exchange is very rich (Gst = 0.01, Nm = 32.22).
The regulating mechanism in COX-III gene difference is still ambiguous between lowland and highland populations. A previous study detected SNPs in three mitochondrially encoded subunit genes of chicken COX, including only one in COX-III (m.10081A>G) between an expanded sample of 56 Tibetan chickens and 152 lowland birds [20]. Another study identified [19] eight SNPs, of which five (m.10081A>G, m.10115G>A, m.10270G>A, m.10336A>G, and m.10447C>T) showed significant differences between Tibetan chickens and lowland chickens. Only the synonymous mutation m.10081A>G was found to differ between haplotypes H4 and H5, and chickens with the A allele at m.10081A>G had a probability of being over 2.6 times better adapted to hypoxia than those with the G allele indicating that m.10081A>G may be a prerequisite for shaping high altitude adaptation-specific haplotypes. Our focus on COX-III SNPs to explore the different haplotypes detected a novel mutation associated with high altitude adaptations. AMOVA showed that COX-III variation mainly existed within a population. Of the 12 defined haplotypes, the existence of Ha4, Ha7, Ha8, and Ha10 only in highland chickens and Ha9, Ha11, and Ha12 only in lowland chickens indicates different degrees of genetic divergence between Tibetan chickens and lowland chickens. Ha1, Ha2, and Ha4 were found to be the earliest ancestors (Figure 1), while haplotype Ha1 was common to all populations suggesting that it is more stable and capable of adapting to new environmental selection. Moreover, Ha2 had significant relationship with high altitude adaptation ( value, 1.911 × 10 −14 * ; OR, 30.375, 95% CI, 7.332-125.837), with the C allele at m.10587 T>C was found to have a probability of being over 30.375 times better adapted to hypoxia than the T allele. We propose that m.10587 T>C affects the functions of the COX enzyme in a similar way to the effect of m.10081A>G [19] on high altitude adaptation of the Tibetan chicken. However, haplotype Ha1 simultaneously contained mutations m.10081A>G and m.10587 T>C and was negatively associated with high altitude adaptation. The function of these mutations warrants need to further research.
Three of the four nonsynonymous mutations identified in the present study (m.10017C>A, m.10115G>A, and m.10555G>A) were peculiar to highland chickens. Mutation m.10115G>A (with the allele frequency of 9.7%, * < 0.05), shared by populations Ganzhi and Diqing, caused the amino acid mutation Val62Ile (Figure 1). A previous study reported SNP m.10115G>A as an uncommon missense mutation in the Tibetan chicken mtDNA genome [30], while another study found that more than one-third of Tibetan chickens (44/125, 35.2%) harboured this mutation, suggesting that it might be associated with high altitude adaptation [19]. In the current study, we found that mutation m.10115G>A is a key site for phospholipid binding, suggesting that it impacts on the function of the mitochondrial membrane in Tibetan chicken. Whereas the other two nonsynonymous mutations (m.10017C>A m.10555G>A) and the reported SNP m.10065T>C [19] both reduced hydropathy index to a certain extent.

Conclusions
In conclusion, we found largest genetic differentiation between Tibetan and lowland breeds and identified haplotype Ha2 is associated with Tibetan chicken populations. The possible association between increased hydrophobicity/reduced hydrophilic characteristics of the mitochondrial membrane and high altitude adaptation could provide a theoretical reference for poultry genetics. Therefore, our results provide a theoretical basis for future research into fowl breeding.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.