Mitochondrial DNA diversity and demographic history of Black-boned chickens in China

Abstract Black-boned chickens (Gallus domesticus, herein abbreviated BBCs) are well known for their unique appearance and medicinal properties and have a long breeding history in China. However, the genetic diversity and demographic history of BBCs remain unclear. In this study, we analyzed 844 mitochondrial DNA D-loop sequences, including 346 de novo sequences and 498 previously published sequences from 20 BBC breeds. We detected a generally high level of genetic diversity among the BBCs, with average haplotype and nucleotide diversities of 0.917 ± 0.0049 and 0.01422, respectively. Nucleotide diversity was highest in populations from Southwest China (0.01549 ± 0.00026), particularly in Yunnan Province (0.01624 ± 0.00025). Significant genetic divergence was detected between most breeds, particularly between Yunnan chickens and those from all other provinces. Haplogroups F and G had the highest levels of genetic diversity and were restricted to Southwest China, particularly Yunnan Province. Based on neutrality tests and mismatch distribution analyses, we did not obtain evidence for rapid population expansions and observed similar demographic histories in BBCs and local non-BBCs. Our results suggest that Chinese BBCs have complex breeding histories and may be selected in situ from local domestic chickens. These results improve our understanding of the genetic heritage and breeding histories of these desirable chickens.


Introduction
Black-boned chickens (Gallus gallus domesticus; herein abbreviated BBCs), renowned for their characteristic black skin, bone, and muscle, have a long history, appearing in the ancient Chinese herbology volume "Bencao Gangmu" written around 1578 C.E. (Li 2005). In China, over 18 BBC breeds have been recorded across 11 provinces (China National Commission of Animal Genetic Resources 2011). For example, Silkie, described in the thirteenth century travelogue "The Travels of Marco Polo" (Benedetto 2014), is a famous breed defined by a set of ten traits: walnut-shaped comb, dark wattles, turquoise-blue earlobes, bearded, silky feathers, five toes, black skin, bones, and meat, and booted feet (China National Commission of Animal Genetic Resources 2011). Studies of BBCs have focused on the chemical properties of the meat (Jaturasitha et al. 2008;Tian et al. 2011), distribution of melanin pigmentation (Nganvongpanit et al. 2020), molecular mechanism underlying melanin deposition (Dorshorst et al. 2011;Shinomiya et al. 2012;Yu et al. 2018;Li et al. 2019), as well as the origin and evolution of fibromelanosis Sohn et al. 2018). Additionally, genetic analyses based on genetic markers, such as microsatellites (Tang et al. 2005;Qu et al. 2006;Yu et al. 2006) and mitochondrial DNA (mtDNA) (Zhu et al. 2014;Guo et al. 2017;Jia et al. 2017;Zhang et al. 2018;Liu et al. 2018;Weng et al. 2019) have attempted to unravel the population genetic history of BBCs. However, the scarcity of sampled breeds is a major limitation of these studies, and this issue is further compounded by the complexity of chicken demographics and domestication histories (Tixier-Boichard et al. 2011;Miao et al. 2013;Lan et al. 2017;Huang et al. 2018a). Consequently, the patterns of genetic diversity and population history of BBCs, including the degree of similarity in demographic trajectories among breeds endemic to different regions, remain unclear.
Since the 1990s, mtDNA D-loop has been widely used to trace the history of chicken domestication owing to its high mutation rate, lack of recombination, and maternal inheritance (Di Lorenzo et al. 2015;Lan et al. 2017). In the current study, we evaluated an extensive mtDNA D-loop dataset for various BBC breeds spanning a wide geographical distribution to reassess genetic diversity and divergence across China. Additionally, a D-loop dataset for local non-BBC populations was used for comparative analyses to infer the demographic histories of BBCs.

Sample collection
A total of 346 wing-vein blood samples were collected and stored at À80 C in the sample database of Jiaying University. These samples were obtained from 13 BBC breeds distributed across nine provinces in China (Table S1). Animal handling and experimentation followed the animal experimental procedures and guidelines approved by the Ethics Committee of Jiaying University (#20151103). Genomic DNA was extracted using the standard phenol-chloroform method, and the DNA concentration and purity were assessed by gel electrophoresis and using a NanoDrop Spectrophotometer 2000 (ThermoFisher Scientific, Waltham, MA). Additionally, 498 published mtDNA D-loop sequences for 15 BBC breeds in China (Table S2) as well as mtDNA Dloop sequences of local non-BBCs (Table S3) were retrieved from GenBank (www.ncbi.nlm.nih.gov).
The PCR products were separated by electrophoresis on a 1.5% agarose gel containing GelRed (Biotium Inc., Fremont, CA), and visualized under ultraviolet light. Sequencing of the D-loop was carried out at Guangzhou IGE Biotechnology Co., Ltd. (Guangzhou, China) using an ABI 3730xl Analyzer (Applied Biosystems, Foster, CA).

Data analyses
DNA sequences were manually checked using BioEdit (Hall 1999), and then aligned using the ClustalW algorithm in MEGA 6.0 (Tamura et al. 2013). All sequences were aligned and trimmed to 518 bp, corresponding to nucleotide positions (nps) 1-518 of the red junglefowl reference sequence AP003321. Haplotypes were defined using DnaSP 6.10.01 (Rozas et al. 2017). These sequences were then assigned to specific haplogroups according to DomeTree (Peng et al. 2015a). Nucleotide (p) and haplotype diversities (Hd) were calculated using DnaSP 6.10.01 (Rozas et al. 2017). A medianjoining network of D-loop sequences was constructed using NETWORK 5.0 (Bandelt et al. 1999). To assess population genetic differentiation and gene flow among sampling locations, pairwise F ST values were computed using ARLEQUIN 3.5 with 10,000 permutations (Excoffier and Lischer 2010). These pairwise F ST values were then used to generate non-metric multidimensional scaling plots using SPSS 19.0 (IBM Corp., Armonk, NY). To further explore geographical structure, an analysis of molecular variance (AMOVA) was also calculated in ARLEQUIN 3.5 (Excoffier and Lischer 2010) using breeds (populations), provinces or regions as groups. A mismatch distribution analysis with the expectations of a sudden expansion model was performed and population expansion statistics (Fu's Fs and Tajima's D) were calculated using ARLEQUIN 3.5 (Excoffier and Lischer 2010) with 10,000 permutations.

Geographical distribution of mitochondrial haplogroups
Our analysis included 844 mtDNA D-loop sequences (346 de novo; GenBank accession numbers: MH923578-MH923923) belonging to 20 BBC breeds in China. In total, 116 haplotypes were identified (nps 1-518), including 90 haplotypes exclusive to a single breed and 26 haplotypes shared by two or more breeds. Yunnan Province had the most haplotypes (54), followed by Zhejiang (33), Guizhou (24), and Sichuan (15). BBCs from Yunnan and Zhejiang provinces had more private haplotypes (Table S4). The 116 haplotypes belonged to haplogroups A-G and Z (Table S4). The geographical distribution of mtDNA haplogroups of BBCs is shown in Table S5. Interestingly, a similar pattern (haplogroups F and G) was also observed in non-BBCs in China (Table S6). Overall, all seven haplogroups were found in Yunnan Province, whereas the other provinces lacked two or more haplogroups.
A median-joining phylogenetic network based on 844 mtDNA D-loop sequences showed a clear geographical pattern ( Figure 1). Haplogroups A, B, C, E, and G formed star-like patterns in the network, each of which had a predominant haplotype and few derived haplotypes. Haplogroups D and Z were rare, only occurring in Henan, Jiangxi, and Zhejiang Province. Haplogroup F, detected only in Yunnan Province, was found at a particularly high frequency in the Tengchong white breed, while haplogroup G was mainly observed in the Yanjin Black-boned breed, and only one haplotype was detected in each of Sichuan Silky and Sichuan Mountain Black-boned breeds. The median-joining network of mtDNA D-loop sequences for BBCs (n ¼ 844) and local non-BBCs (n ¼ 1663) did not demonstrate BBC-specific or non-BBC-specific patterns, and the two groups shared the same haplogroups or major haplotypes ( Figure S1).

Genetic diversity
The Hd and p values, and the average number of nucleotide differences (K) for all samples were 0.917 ± 0.0049, 0.01442, and 7.365, respectively (Table 1). The Wuliangshan Blackboned breed had the highest p value (0.01633 ± 0.00048), while the Silkies breed had the lowest (0.00733 ± 0.00063). With respect to the 11 provinces, BBCs from Yunnan Province had the highest p value (0.01624 ± 0.00025), while those from Jiangxi Province had the lowest (0.01024 ± 0.00062). Regionally, the highest p value was detected in Southwest China (0.01549 ± 0.00026) and the lowest was detected in Northwest China (0.01162 ± 0.00081). Regarding Hd, Yunnan, Guangxi, Henan, and Zhejiang Provinces showed the highest values (>0.9), while Fujian and Jiangxi Provinces had the lowest (<0.8) ( Table 1). Combining the results for Hd and p, diversity was the highest in chickens from Southwest China, especially those from Yunnan Province, and the lowest in chickens from Northwest China. Interestingly, it shows near levels of p in both BBCs and non-BBCs, with the exception of those in Yunnan Province, where it harbors higher p values in BBCs than in non-BBCs (Tables 1 and S7). Estimates of genetic diversity for major haplogroups are presented in Table S8. Briefly, p values were the highest in haplogroups F and G in Southwest China.

Genetic divergence
To evaluate genetic differentiation among chickens, pairwise F ST values were computed using ARLEQUIN 3.5 according to chicken breeds and provinces. A non-metric multidimensional scaling plot was subsequently generated based on the pairwise F ST values to ascertain population/province relationships ( Figure S2). The Yunxian Black-boned breed showed the greatest differentiation from other breeds, while Silkies, Hubei Silky black, and Sichuan Silky breeds were clustered together ( Figure S2A). The BBC populations from Yunnan Province were significantly divergent from those from other provinces ( Figure S2B).
The AMOVA showed that variance was higher within populations than among populations within groups or among groups (Table S9). This suggests that the molecular variance mainly exists within breeds, followed by provinces, and finally regions. These results are generally in accordance with the results of F ST analyses.

Demographic history
Neutrality tests and a mismatch distribution analysis of BBCs and non-BBCs were performed to detect population trajectories of BBCs. As determined by Fu's Fs neutrality statistics for BBCs, only chickens from Yunnan Province and Southwest China showed significant deviations from neutrality, but none of the breeds showed deviations from neutrality (Table 1). All Tajima's D tests for BBCs showed no significant departure from neutrality (Table 1). In the Fu's Fs neutrality tests of non-BBCs, however, departures from neutrality were detected in populations from Yunnan, Henan, and Hubei Provinces (Table S7). The mismatch distribution graph is in accord with the results of neutrality tests ( Figure S3).

Discussion
Based on the most comprehensive dataset to date, including 844 mtDNA D-loop sequences for 20 breeds from 11 provinces across China, we characterized genetic diversity and genetic differentiation in Chinese BBCs, which show interesting geographical patterns. Haplogroups A-G were disproportionately observed among the 844 chickens. Additionally, these haplogroups showed differential geographical endemicity. Common haplogroups, including A, B, and E, were found in all eleven provinces. However, haplogroup C was absent in Northwest China, haplogroup D occurred in Central and East China, while haplogroups F and G were restricted to Southwest China. Interestingly, non-BBCs in China showed a similar geographical distribution pattern to those of haplogroups F and G.
Specific mitochondrial haplogroups are widely used as candidate genetic markers to trace the demographic history of chickens (Lan et al. 2017), including chicken C1 (Huang et al. 2018a) and haplogroup D (Herrera et al. 2017;Zhang et al. 2017). Thus, the geographical distribution observed in the current study suggests that BBC breeds associated with different areas are likely to have distinct demographic histories.
To evaluate this possibility, we performed detailed genetic analyses of mtDNA D-loop sequence data for Chinese BBCs. Compared with that of other indigenous Chinese chickens, such as chicken breeds from Guangdong Province and its adjacent regions (Huang et al. 2018b), Jiangsu Province (Jia et al. 2017), East China (Jia et al. 2017), and other regions in China (Huang et al. 2018b), BBCs maintained a high level of genetic diversity. In particular, genetic diversity of BBCs from Southwest China was higher than diversity estimates for chickens from other regions, whereas BBCs from East China showed the lowest diversity. This observation may be explained by differences in traditional cultures as well as different intensities of selection. BBCs are fascinating breeds owing to their unique appearance, and local residents tend to limit crossbreeding to maintain characters. Silkies is a famous BBC breed with an atypical fluffy plumage; it was officially recognized in 1874 in North America as the Standard of Perfection (Ekarius 2007). Owing to its popularity, Silkies has been subjected to intense selection, resulting in a decline in genetic diversity (Jia et al. 2017). With a relatively large sample size, we detected few exclusive haplotypes and low genetic diversity in Silkies, confirming its breeding history. Neutrality tests indicated that BBCs in China have not undergone rapid population expansion. The mismatch distribution graph showed that BBCs and non-BBCs share similar demographic trajectories, except for those from Henan and Hubei provinces. This finding is in accordance with the results of Huang et al. (2018a). A neighbor-joining network did not indicate BBC-specific haplogroups or major haplotypes. Consistent with the estimates of genetic diversity and genetic divergence and with husbandry reports, numerous BBC breeds endemic to China may be selected in situ from local domestic chickens (Weng et al. 2019).
The significant differentiation between chickens from Southwest China, particularly Yunnan Province, and those from other regions suggests that BBCs of Southwest China experienced limited gene flow with each other and with chicken from other regions, with only weak genetic selection (Miao et al. 2015). Both Fu's Fs and pairwise F ST values indicated that no breed deviates from neutrality, whereas deviations were detected at the province/region level. It is noteworthy that genetic analyses by haplogroups showed that F and G, privatized to Southwest China, have higher p values than those of other haplogroups, including A and C. Although previous research has reported that the chickens from Yunnan Province display high levels of genetic diversity and multiple unique mtDNA haplogroups (Liu et al. 2006;Miao et al. 2013;Huang et al. 2018a), the complex history of chicken domestication (Xiang et al. 2014;Peng et al. 2015b;Peters et al. 2015Peters et al. , 2016Eda et al. 2016;Huang et al. 2018a) necessitates further research using more robust approaches and high-density markers, such as genome-wide SNPs, to establish whether the breeding center of Chinese BBCs is in Southwest China and to obtain a more in-depth understanding of BBC evolution.