Genomic diversity dynamics in conserved chicken populations are revealed by genome-wide SNPs

Background Maintaining maximum genetic diversity and preserving breed viability in conserved populations necessitates the rigorous evaluation of conservation schemes. Three chicken breeds (Baier Yellow Chicken (BEC), Beijing You Chicken (BYC) and Langshan Chicken (LSC)) are currently in conservation programs in China. Changes in genetic diversity were measured by heterozygosity, genomic inbreeding coefficients, and autozygosity, using estimates derived from runs of homozygosity (ROH) that were identified using SNPs. Results Ninety DNA samples were collected from three generations for each breed. In the most recent generation, the highest genetic diversity was observed in LSC, followed by BEC and BYC. Inbreeding coefficients based on ROH for the three breeds declined slightly between the first and middle generations, and then rapidly increased. No inbreeding coefficients exceeded 0.1. Population structure assessments using neighbor-joining tree analysis, principal components analysis, and STRUCTURE analysis indicated that no genetic differentiation existed within breeds. LD decay and ROH at different cut-off lengths were used to identify traces left by recent or ancient inbreeding. Few long ROH were identified, indicating that inbreeding has been largely avoided with the current conservation strategy. The observed losses in genetic diversity and occurrences of inbreeding might be consequences of sub-optimal effective population sizes. Conclusions The conserved Chinese chicken populations have high genomic diversity under the current conservation program (R: F). Furthermore, this study highlights the need to monitor dynamic changes in genetic diversity in conserved breeds over successive generations. Our research provides new insights into genetic diversity dynamics in conserved populations, and lays a solid foundation for improving conservation schemes. Electronic supplementary material The online version of this article (10.1186/s12864-018-4973-6) contains supplementary material, which is available to authorized users.


Background
As one of the earliest centers of domestication for chickens, China has the most abundant chicken genetic resources in the world, with 107 indigenous chicken breeds. These breeds play an essential role in the Chinese poultry industry due to the popularity of traditional cuisine. Because these chickens typically exhibit high adaptability to variable environments, strong disease-resistance, and produce high-quality meat and eggs, they are an important breeding resource to meet future market demands [1]. However, a sizeable fraction of indigenous chicken breeds (21.3% in the world) are under threat of extinction because breeding populations are too small for genetic sustainability [2]. Twenty-three Chinese indigenous chicken breeds have been listed in the national conservation catalogue and are currently managed under what is thought to be an optimal conservation scheme. However, the effectiveness of the scheme has never been evaluated over the long term.
An efficient in situ conservation scheme relies on an effective population size, as well as an effective selection and mating strategy [1]. The recommended effective population size is 50, which is not only sufficient to maintain population fitness, but is also small enough to be monitored and managed easily [1,3]. The principle of a selection and mating strategy is to minimize the average kinship between selected parents [1]. The mating systems that are in use are: (i) random mating and random selection (R: R), (ii) random mating within families, with one son kept per sire family and one daughter kept per dam family (R: F), and (iii) family rotational mating (F: R) [4].
Results from a study that compared the effectiveness of these three mating systems suggest that F: R can sustain 90% of genetic diversity in a livestock population for more than 100 years [4]. Other studies demonstrate that F: R can reduce inbreeding in populations [5,6], however, implementation of the F: R scheme requires substantial effort. The contemporary conservation scheme in China relies mainly on the R: F mating system. Although simulation experiments suggest that R: F and F: R perform similarly in maintaining genetic diversity in a conserved population [4], few empirical studies have evaluated their effectiveness.
The genetic diversity dynamics in successive generations within a conserved population directly reflect the effectiveness of a conservation scheme. DNA markers can be used efficiently to estimate the genetic diversity within and between conserved populations [7]. Various DNA marker systems have been used to assess chicken genetic diversity, such as RAPDs [8], AFLPs [9], and microsatellites [10][11][12]. SNPs, which are densely distributed across genomes, have also been used to estimate genetic diversity and population structure with high accuracy [13,14].
In order to investigate the effectiveness of the schemes that are currently used to conserve chicken genetic resources in China, we evaluated the genetic diversity of three conserved indigenous chicken breeds. SNPs, identified using high-throughput DNA sequencing in samples obtained from three generations per breed, were used to estimate diversity and track changes across time. These data can be used to evaluate and improve ongoing conservation efforts.

Ethics statement
Sample collection procedures strictly followed protocols approved by the Animal Welfare Committee of China Agricultural University (Approval Number: XK257).

Sampling
Three Chinese indigenous chicken breeds that have been enrolled in conservation programs were used in this study: Baier Yellow Chicken (BEC), Beijing You Chicken (BYC), and Langshan Chicken (LSC). Different geographical and environmental factors have contributed to the unique characteristics of these breeds. The Baier Yellow chicken, mainly produced in Jiangxi Province, has a distinct appearance with white ears and yellow feathers, beak, and shanks. It is a rare prematurity and egg-type breed. The Langshan chicken is a classic dual-purpose breed that originated in Jiangsu Province. This is one of the oldest breeds and is unusually tall, with long legs and a tail carried at a high angle. The Beijing You chicken is an ancient breed that originated during the Qing Dynasty in Beijing. Known as "royal chickens", they are valued for their high-quality meat and eggs, and are uniquely marked by a crest on the head, a beard under the lower jaw, and feathers on both shanks [15][16][17].
These breeds have been managed for optimal conservation as part of the National Chicken Genetics Resources program (Jiangsu, NCGR). Briefly, the conservation goals are that population sizes should be kept constant across generations (30 males and 300 females), and random mating should be enforced within families, with one son kept per sire family and one daughter kept per dam family (R: F). Samples from 270 individuals were collected from NCGR, with three generations per breed and 30 individuals per generation. Samples are identified by breed abbreviation and the last two digits of the year in which samples were obtained (e.g., samples collected from BEC in 2007 are designated BEC07). Breed and sampling information is summarized in Table 1. Blood samples were collected from the wing vein and stored at − 20°C. Genomic DNA was extracted following the protocol accompanying the DNeasy Blood & Tissue Kit (Qiagen Inc., Valencia, California, USA). 3 μg high quality DNA was used to construct sequencing libraries for each sample.

Inbreeding coefficient
Two measures of inbreeding coefficient were calculated for each chicken population.
Inbreeding coefficient based on the mating plan (F ES ): The estimation of effective population size (Ne) was based on number of sires and dams, following Wright's model [23]. Computation of Ne requires the numbers of males (N m ) and females (N f ) in each population that participated in the R: F program, and is calculated using the equation: Ne ¼ 3Nf þNm 16NmNf . The increment of hypothetical inbreeding (ΔF) was calculated using the equation: Inbreeding coefficient based on runs of homozygosity (F ROH ): A run of homozygosity is defined as a region > 100 Kb containing > 50 SNPs. F ROH was calculated using PLINK v1.90 (with parameters --file BEC07_qc --ibc --allow-extra-chr --chr-set 28 --out BEC07) and is the fraction of the genome spanned by runs of homozygosity [24].

Analysis of population structure
To reduce noise due to linkage disequilibrium, SNPs with a pair-wise genotype r 2 value ≥0.2 were removed from the data set. A principal component analysis (PCA) [25] was conducted using PLINK and visualized with the SNPRelate R package [26]. A neighbor-joining (NJ) tree was constructed with Nei's genetic distances [27] using the phylogeny program MEGA v7.0 [28] and displayed with FigTree v1.4.3 [29]. The genetic structures of the 9 sub-populations described above were analyzed with STRUCTURE v2.3.4 [30], using admixture and a correlated allele model [30,31]. Ten independent runs were performed with K ranging from 1 to 10, with a burn-in period length of 10,000, followed by 100,000 Markov chain Monte Carlo (MCMC) repetitions, and 20 replications for each K value. STRUCTURE HARVESTER [32] was utilized to determine the optimal K value by comparing the likelihood of the data (LnK) for different values of K [lnP(X|K)] and by examining the second-order rate change of lnP(X|K),ΔK [33,34]. Results for K = 2 to K = 9 are included in this report.

Linkage disequilibrium decay
LD was evaluated as the correlation coefficient (r 2 ) between alleles at two separate SNP loci [37]. Within each population, all pairs of autosomal SNPs with MAF > 0.05 and Hardy-Weinberg equilibrium P-value >10E-6 were used to calculate r 2 with Haploview [38]. Inter-SNP distances from 0 kb to 500 kb were consolidated into 5 bins.

Effective population size
Effective population size (Ne) was estimated according to the random mating model of linkage disequilibrium, using default parameters in N E ESTIMATOR v2.01 [39].
Ne estimates for each breed were calculated as the average of the estimates for macrochromosomes (gga1-gga5) [40] (Axelsson et al., 2005).

Runs of homozygosity
Runs of homozygosity (ROH) were identified for each of the 9 sub-populations using PLINK v1.90. The ROH program slides a moving window of 1 Mb across the genome to estimate homozygosity. One heterozygous and five missing calls per window were allowed to avoid false negatives caused by occasional genotyping errors or missing genotypes. The minimum length and SNP counts required for each ROH were 100 kb and 50 SNPs, respectively. Additional statistical significant tests were conducted to detect the differences in genome-wide homozygosity levels among populations with three measures (NSEG, KB, KBAVG).

Descriptive statistics
To assess genetic diversity, DNA samples from three indigenous chicken breeds (BEC, BYC, and LSC) were subjected to high-throughput DNA sequencing.  Table S1). The distribution of SNPs across all chromosomes is shown in Additional file 2: Figure S1.
Genetic diversity within the BEC, BYC and LSC breeds All three breeds maintained relatively high genetic diversity in the R: F conservation scheme (see scheme definitions in Materials and Methods). LSC exhibited the highest genetic diversity as measured by Ho (0.2348), He (0.2250), A R (1.235), and P N (0.8130) ( Table 2). An analysis across three generations within the same breed showed that genetic diversity was highest for LSC15 with Ho (0.2379), He (0.2281) and A R (1.238). The highest proportion of polymorphic markers (P N ) was observed in BEC07 (82.41%), while BYC15 exhibited the lowest genetic diversity. As expected, BEC and BYC showed decreasing levels of diversity from the initial generation (BEC07/BYC07) to the current generation (BEC15/BYC15). Conversely, LSC displayed increasing diversity with the implementation of a conservation program (Fig. 1). However, dynamic changes in genetic diversity within breeds were less than 10% throughout the sampled generations (Additional file 3: Figure S2).

Estimation of inbreeding coefficients
Estimated inbreeding coefficients varied between breeds and conservation methods. Average F ES ranged from 0.0789 in BEC to 0.2010 in BYC. As expected, F ES values across generations increased as conservation procedures were maintained. In contrast, average F ROH tended to be lower than F ES , ranging from 0.0511 in BEC to 0.0745 in BYC, and F ROH did not exhibit the steady increase that was observed for F ES . Maximum inbreeding was observed in BYC15 (F ES = 0.2010 and F ROH = 0.0925). Correlation between F ES and F ROH was strongly positive (r 2 = 0.76).

Population structure analysis
Population structures of the three native chicken breeds, comprising 9 conservation sub-populations, were analyzed using PCA, NJ tree, and STRUCTURE. PCA showed that the first two principal components account for 17.77% (PC1) and 15.01% (PC2) of the total variability. Individuals from the 9 sub-populations clearly group into their respective breeds (Fig. 2a). The results of the , inbreeding coefficient based on runs of homozygosity NJ tree analysis were consistent with those obtained by PCA (Fig. 2b). Figure 2c shows an admixture plot representing the 270 sampled chickens, generated using a model-based clustering approach. At a low value of K (K = 2), two distinct ancestors are apparent (BYC and LSC). BEC appears to include both LSC (as the majority component) and BYC. At K = 3, individuals cluster strongly into the three corresponding breeds, consistent with the PCA and NJ tree results (as shown in Fig. 2a and b). All generations within each breed show the same pattern. The optimum population structure inferred using the admixture model in STRUCTURE was subdivided into three sub-populations based on both LnP(D) and Evanno's ΔK method (K = 3; Fig. 3). At K = 4, BEC splits into its two main ancestors. For K = 5 to 8, LSC appears to include two or more distinct ancestors, but at K = 9, it groups again into one common ancestor. Finally, the, BYC breed always exhibits homogeneity, except for K = 9.

Population differentiation analysis
To investigate the extent of population differentiation between different generations within breeds and between breeds, F ST values were calculated using the filtered genotype data (Additional file 4: Table S2). The F ST values for all pair-wise population comparisons are shown in Fig. 4a. For the entire population, F ST values varied from 0.0046 to 0.1530, and F ST values between breeds ranged from 0.1127 to 0.1243 (Fig. 4b). F ST values are expected to be significantly higher between breeds than between generations within a breed. All F ST values between generations within a breed were below 0.05 (from 0.0046 to 0.0423), indicating that no obvious genetic differentiation appeared within any breed (Fig. 4a). However, these F ST values increased during conservation. For example, for BEC, the F ST value was 0.0046 between BEC07 and BEC10, and increased to 0.0329 between BEC07 and BEC15. Similar trends were observed between BEC07 vs BEC15 (0.0329) and BEC10 vs BEC15 (0.0285).
Linkage disequilibrium decay (LD decay) and effective population size LD for each population was estimated as the physical genomic distance at which the genotypic association (r 2 ) decays to less than half of its maximum value. Short-range LD was always observed in each of the three different generations within same breed (Fig. 5). As expected, LD values tended to increase as conservation continued. For example, LD values for LSC10, 12 and 15 were 13.19 kb, 16.99 kb, and 20.10 kb, respectively. The BEC pattern was similar: 11.56 kb (BEC07), 13.55 kb (BEC10), and 20.77 kb (BEC15). The highest LD value in BYC occurred in the last generation (BYC15, 24.61 kb), but BYC07 (23.99 kb) exceeded BYC10 (17.95 kb).
Effective population size (Ne) was estimated for autosomal chromosomes gga1 through gga28 based on linkage disequilibrium (Table 3). Average Ne differed amongst the 9 sub-populations (from 19.33 to 34.85). Within macrochromosomes (gga1-5), BEC07 had the highest estimated Ne (81.52) within a breed, and Ne declined in BEC as conservation continued. In contrast, Ne was lower (55.36) in LSC12 than in LSC10 (70.28) or LSC15 (73.74). Ne estimates for BYC fluctuated, with a high value (75.04) in BYC10 and lower values in BYC07 and BYC15.

Runs of homozygosity (ROH)
Runs of homozygosity (ROH) were identified in the genomes of the 9 sub-populations from all three breeds (Table 4). A genome-wide survey for autozygosity was conducted to identify regions with signatures that reflect ancient or recent inbreeding effects. We estimated F ROH , and found that the maximum values occurred in the  Table 2). BYC15 had the highest level of inbreeding (0.0925), while the BEC and LSC breeds had similar and lower inbreeding levels (~0.05).
ROH was then assessed to determine whether any populations exhibited evidence of recent inbreeding. All three generations in BYC had higher ROH values, suggesting that recent inbreeding had occurred in this breed (Table 4 and Fig. 6). Consistent with the LD decay analysis, the highest ROH was observed in BYC15 (r 2 0.1248 = 24.61 kb), followed by BYC07 and BYC10. We speculate that BYC inbreeding might reflect the small effective population size of this breed ( Table 3). The ROH values for BEC were highest in BEC15, followed by BEC07 and BEC10, which parallels the trend observed for LD decay and F ROH . Similar patterns were observed for LSC. Homozygosity was also measured between the individuals in each sub-population using three methods (NSEG, KB, and KB AVER ) (Additional file 5: Table S3). All three measures varied significantly between BYC generations (P<0.01). NSEG and KB also showed significant differences (P<0.0001) between generations in BEC and NLS.

Discussion
The chicken, one of the first animals to be domesticated, has been subjected to long-term natural selection, artificial selection, and genetic drift for diverse specific traits [41,42]. A variety of factors accelerated the generation of phenotypic differences and genetic variability [43,44]. However, this variability has been threatened due to ecosystem damage and commercial breeding. In China, gene pool (live) conservation and protected regions for both ex situ and in situ conservation have been established for the management of poultry genetic resources, as exemplified by the National Chicken Genetics Resources Program (Jiangsu). Because live conservation is typically implemented using small populations, it is necessary to monitor the status of each population to evaluate the effectiveness of the management strategy. In this study, we performed genotyping by sequencing (GBS) to assess the genomic diversity of different generations from three conserved breeds (Baier Yellow Chicken, Beijing You Chicken and Langshan Chicken).
The majority of studies have estimated genetic variability in Chinese indigenous chickens using microsatellites [45,46] or mtDNA [47,48], but genome-wide SNPs have seldomly been used. We estimated the genetic diversity in three chicken populations using SNP markers. The results showed that all three chicken breeds have maintained rich genetic diversity in terms of heterozygosity (Ho, He), proportion of polymorphic markers (P N ), and allelic richness (A R ), consistent with previous studies [45,46,[49][50][51]. In the most recent generations     Fig. 6 Homozygosity frequency distribution derived from runs of homozygosity (ROH) for each generation and breed sampled, LSC15 ranked first in genetic diversity, followed by BEC15 and BYC15 (Table 2). A study of the same populations in 2008 indicated that genetic diversity measured using microsatellites was highest in BYC, followed by LSC and BEC [45], suggesting that genetic diversity in BYC has decreased more rapidly than in the other breeds. We also observed this declining trend in the BYC breed. BYC diversity decreased from 2007 to 2015 (Table 2), while genetic variability in the BEC breed fluctuated, and LSC exhibited a slight increase. The BYC breed has been under conservation (~39 generations) for a longer period than either LSC or BEC, which have been conserved only since 1998 (~17 generations). The long-term practice of conservation in a small population size may reduce genetic diversity. Furthermore, all three breeds were subjected to ex situ live conservation in Jiangsu. The BYC breed, which originated in Beijing, might have adapted poorly to the environment, resulting in a loss of genetic diversity. In contrast, the LSC and BEC breeds might have adapted more easily.
The genetic diversity in all breeds changed no more than 10% between generations (Additional file 3: Figure  S2). The conservation goal is to maintain 90% of the genetic diversity from the initial population and an inbreeding coefficient less than 0.1 for 100 years [52,53]. According to our results, the genetic diversity of the three chicken populations meets conservation criteria under the current program (R: F). In particular, inbreeding events have been effectively avoided under the R: F mating system, based on assessments of population structure, genetic differentiation, LD decay, and ROH.
Nevertheless, the decline of genetic diversity should not be ignored ( Fig. 1 and Table 2). The significant differences in ROH that we observe between generations in all three breeds also suggest that these populations have not reached the desired level of genetic stability during conservation. Both the decline in genetic diversity and the high heterozygosity across generations are indicative of genetic drift, which can be reduced by enlarging the population size. In our study, the estimated effective population sizes (Ne), based on whole-genome SNPs for the conserved populations, were far below the required threshold of 50 individuals [1,3]. We also evaluated Ne according to chromosome size, using the classification proposed by the International Chicken Genome Sequencing Consortium [54]: large macrochromosomes (gga1-5), intermediate chromosomes (gga6-10) and micro-chromosomes (gga11-28). Because micro-chromosomes have high rates of recombination, we estimated Ne based on the macrochromosome class (gga1-5). The maximum Ne was 81.52 in BEC07 and the minimum was 41.56 in BEC15 (Table 3), suggesting these conserved populations are relatively stable but also at risk. We therefore recommend that ex situ live and in situ live conservation efforts be combined to help maintain high levels of genetic diversity in the long term.

Conclusions
In summary, we collected 270 samples from three successive generations of three conserved chicken breeds. We estimated dynamic changes in genetic diversity using genome-wide SNPs, making it possible to comprehensively evaluate the current conservation scheme (R: F). The results demonstrated that the conserved Chinese chicken populations have sustained high levels of genetic variability under current conservation practices. We also compared successive generations within each breed to characterize trends in genetic diversity, allowing us to assess the effects of conservation over time. Overall, this study demonstrates an efficient strategy for assessing the success a conservation program and for improving conservation and management practices.