Whole genome sequencing revealed genetic diversity 2 and selection of Guangxi indigenous chickens 3

Guangxi chickens play an important role in promoting the high-quality development of the broiler industry in China, but their value and potential are yet to be discovered. To determine the genetic diversity and population structure of Guangxi indigenous chicken, we analyzed the whole genomes of 185 chicken from 8 phenotypically and geographically representative Guangxi chicken breeds, together with 12 RJFt, 12 BRA and 12 WL genomes available from previous studies. Calculation of heterozygosity (Hp), nucleotide diversity (π), and LD level indicated that Guangxi populations were characterized by higher genetic diversity and lower differentiation than RJFt and commercial breeds except HGFC. Population structure analysis also confirmed the introgression from commercial broiler breeds. Each population clustered together while the overall differentiation was small, MA has the richest genetic diversity among all varieties. Selective sweep analysis revealed BCO2, EDN3 and other candidate genes have received strong selection in local breeds, these also provided novel breeding visual and data basis for future breeding.


30
Chickens are the most widely distributed livestock species globally, more than half 31 the total (53%) is found in Asia, one of the largest producers is China [1] . In China, 32 poultry meat consumption accounts for the second largest proportion after pork. People  only sporadic studies were reported on genomic information about them [2,3] . 46 A comprehensive and deep understanding of the genome diversity of the 47 indigenous breeds could reveal the genetic diversity and population structure of these 48 breeds. This study therefore investigated genetic diversity, population structure, linkage 49 disequilibrium (LD), and signature selection within Guangxi indigenous chickens using 50 genome-wide single nucleotide polymorphisms (SNPs) generated from the whole 51 genome sequencing.

52
Ethics statement 53 This study was carried out in accordance with the guidelines of the Animal 54 Experimental Ethical Inspection Form of Guangxi Research Institute (20190318).

56
Sampling and genotyping 57 A total of 185 blood samples from six breeds and two characteristic populations 58 were investigated from conservation center or breeding farms (S1 Table and  the total ROHs length of each chromosome was centered and scaled in breed's level.

83
Population structure analysis 84 To investigate the genetic background of the chickens, principal component 85 analysis (PCA) and structure analysis were conducted. SNPs in high linkage 86 disequilibrium were removed by PLINK version 1.9 [8] , based on the pruned SNP data, 87 the individual ancestries were estimated using a maximum likelihood method 88 implemented in ADMIXTURE version 1.23 [9] . region, and only 2.27% located in the exon region ( Fig 1B).

114
The genome is divided into isochrones with a sliding window of 100kb, and 115 divided into five categories (L1, L2, H1, H2 and H3) according to different GC level to 116 explore the potential impact between GC content and genetic variations ( Fig 1C and   117 1D) [12] . Our results shown that the L2 category has the largest number of isochrones, 118 covering 37% of the genomic region, and the SNPs and Indels counts peak in this 119 category. H1 category with higher GC level also contains a lot of genetic variations (S5 120   Table). In general, genomic regions with moderate GC content contain more variation. published [13] . The nucleotide diversity (π) and heterozygosity were calculated to 132 evaluate the genetic diversity of all the chicken breeds. We observed Guangxi 133 indigenous chickens harbored the higher genome-wide π than RJFt (π =0.00334) except 134 for DZAC (π = 0.00332), the lowest genome-wide π in WL (π = 0.00152), followed by 135 BRA (π = 0.0031) (Fig 2A). Unlike the results in nucleotide diversity, BRA harbored undergone higher selective pressure than autosome [14] . Linkage disequilibrium (LD) 140 analysis showed that WL population had the slowest LD decay rate, significantly slower 141 than the followed BRA. MA had a faster LD decay rate than other chicken breeds, 142 DZAC and WC have similar LD level with RJFt in second group (Fig 2C).

143
The level of ROH reflects the recent inbreeding history of a population [15] . As 144 shown in Fig 2D,   (rs313409504) was consisted with the previous report [16] . The strongest selective sweep by an asymmetric cleavage reaction [19] , is established as the causal gene for the yellow inhibit expression of BCDO2 in skin caused yellow skin, but not in other tissues [16] .

259
Fallahshahroudi's study showed the down-regulation of BCO2 in skin, muscle, and 260 adipose tissue was associated with the derived haplotype [20] . Also, BCO2 has variety 261 variants in different breeds. Wang found a G>A mutation in exon 6 to be associated 262 with the concentration of carotenoids in Guangxi-huang and Qingjiao-ma chicken [21] .

263
A GAG haplotype was fixed in commercial breeds of yellow skin [16] . We also found 264 the missense mutant at chr24:6155481 led to the mutation of threonine to alanine. adipogenesis in mammals [26] . Zhang et al.'s study provides strong in vivo evidence that 284 atg7, and by inference autophagy, is critical for normal adipogenesis [27] . AMACR 285 coding protein alpha-methylacyl-CoA racemase, this protein is involved in the pathway 286 bile acid biosynthesis, which is part of Lipid metabolism (gga00120). Bile acid is the 287 main component of bile and its main function is to promote the digestion and absorption 288 of fat. HSD17B4 codes a bifunctional enzyme mediating dehydrogenation and 289 anhydration during β-oxidation of long-chain fatty acids, and a non-synonymous SNP 290 has been reported to be related to meat-quality traits in pig [28] . PRKAA1 is associated 291 with skeletal muscle lipid accumulation [29] . and red jungle fowls were downloaded from NCBI at ERP112703 (S2 Table).