Characterisation of runs of homozygosity and inbreeding coefficients in the red-brown Korean native chickens

Ethical statement

The care and handling of animals used in this study were approved by the Institute of Animal Care and Use Committee of the National Institute of Animal Science (approval No: NIAS 20212219) and Chungnam National University’s Animal Ethics Committee.

Study samples, genotyping, and quality control

The study animals were obtained from the red-brown Korean native chicken pure line reared at the National Institute of Animal Science, Republic of Korea. The study sample consisted of a total of 651 chickens sampled from four consecutive generations between the years 2019 and 2022. The number of samples per generation is shown in Table 1. Genomic DNA was then extracted from the blood samples using a commercial kit (GeNet Bio, Daejeon, Korea). The DNA samples were genotyped using the Illumina chicken 60K SNP chip (Illumina Inc., San Diego, CA, USA) containing 57,636 SNPs [22]. Sex chromosomes and markers that were not assigned to specific chromosomes were excluded from further analyses. Only SNPs with a call rate (>0.9) and Handy-Weinberg equilibrium (p-value >1.0×10⁻⁶) were used for the analysis. No further quality control was conducted during the analysis of ROH, as recommended by Meyermans et al [23]. On the other hand, genotypes with a minor allele frequency of less than 0.05 were excluded during the estimation of other inbreeding coefficients. Consequently, 53,872 SNPs were used for ROH analysis while 44,569 SNPs were available for the calculation of other inbreeding coefficients.

Table 1

Number of samples per generation

Analysis of runs of homozygosity

The calling of ROH was conducted using the (--homozyg) function in PLINK v1.9 software [24]. The analysis parameters were set to ensure maximum genome coverage and minimise calling errors. Consequently, ROH segments were called based on the following parameters: a minimum ROH length of 1,000 Kb, a sliding window with a minimum of 49 SNPs, one heterozygous SNP per ROH segment, five missing SNPs per ROH segment, a maximum gap of 1,500 Kb between consecutive homozygous SNPs and the threshold of overlapping homozygous windows was set to 0.05. The minimum number of SNPs in the sliding window was estimated based on the formula proposed by Lencz et al [25] as follows:

Where, n_s is the number of genotyped SNPs per sample, n_i is the number of genotyped individuals, α is the percentage of false positive ROH (0.05) and het is the average heterozygosity across all genotyped SNPs. The ROH segments were grouped into classes based on length. The ROH categories were: 1 to 2, 2 to 4, 4 to 8, 8 to 16 and >16 Mb. The patterns of ROH were reported in terms of the mean ROH length per individual, the minimum, average and maximum length of ROH segments, and the frequency of ROH segments in different ROH length classes. The ROH summary statistics were calculated using the detectRUNS package in R software [26]. The ROH pattern across the four generations was also scrutinized.

Estimation of inbreeding coefficients

The coefficient of inbreeding based on ROH was calculated using a formula employed by McQuillan et al [27] as follows:

Where, L_ROH is the sum of the length of ROH segments in the genome of each bird and L_AUTO is the total length of the autosomal genome. The genome length was estimated to be 0.94 Gb. The F_ROH was also calculated for the different ROH classes defined based on the ROH length.

Additionally, three other measurements of inbreeding were estimated using the (--ibc) function in the GCTA software [28]. These include:

Inbreeding coefficient calculated from the diagonals of the genomic kinship matrix (F_GRM) [29] as follows:

Where, Z is a covariance matrix constructed by subtracting twice the minor allele frequency from the raw marker score. The covariance matrix consists of values: (0-2p) for homozygous, (1-2p) for heterozygous and (2-2p) for alternate homozygous loci, where p is the minor allele frequency.

The Wright’s inbreeding coefficient (F_HOM) which compares the expected and the observed homozygous genotypes for each sample [28] as follows:

Wright’s inbreeding coefficients calculated as the correlation between alleles in uniting gametes (F_UNI) [28]:

Where, x is the number of versions of the reference allele, p and q are the reference and alternative allele frequencies of mth SNP respectively. The correlation between F_ROH coefficients and other inbreeding coefficients was estimated using Pearson’s correlation and plotted using the Corrplot package in R software [26]. The trend in average inbreeding coefficients across the four generations was also evaluated.

Identification and functional annotation of genes and quantitative trait locus in runs of homozygosity islands

Genomic regions with common ROH across the population were extracted using the detectRUNS package in R software. The frequency of occurrence of SNPs in the ROH segments was estimated by enumerating the number of times each SNP was detected in the ROH. Genomic regions where the incidence of SNP in ROH exceeded a population-specific threshold were considered ROH islands. The threshold was calculated as described by Gorssen et al [30]. Accordingly, the distribution of ROH incidences was standardized using a standard normal Z-score and then the top 0.1% of SNPs with a p>0.999 were selected to form the ROH islands. The GALLO package in R software was used to identify annotated genes and QTLs in the ROH islands [31]. The QTLs in the ROH islands were annotated based on the animal QTL annotation database [32]. Similarly, candidate genes in the ROH islands were annotated based on the GRCg6a reference genome assembly. Finally, a literature search was conducted to elucidate the functions of the identified genes.

Runs of homozygosity

Genotype data was utilised to characterise the ROH and estimate inbreeding coefficients in the red-brown Korean native chicken. Figure 1 shows the distribution of the length and numbers of ROH segments. A total of 32,176 ROH segments were called, with an average of 49 segments per bird. The average length of the ROH segments was 2.7 Mb. The longest ROH segment was found on Gallus gallus chromosome (GGA) 2 with a length of 33.47 Mb, while the shortest segment was identified on GGA1 with a length of 1.00 Mb. GGA1 had the highest number of ROH segments (n = 6,832) while GGA2 had the greatest average ROH length (3.2 Mb) (Figures 2A–2B). The overall average ROH length was 134.12 Mb (Figure 1A).

Figure 1

Comparison of the length and number of ROH segments in red-brown Korean native chickens. (A) Total sum of ROH for individuals; (B)Scatterplot of number of ROH in individuals versus sum ROH length; (C) Average length of ROH segments for individuals; (D) Scatterplot of the number of ROH for individuals versus the average length of ROH segment. ROH, runs of homozygosity.

Figure 2

Characterization of ROH in chromosomes. (A) The bars represent the number of ROH while the line represents the percentage of ROH per chromosome; (B) Average length of ROH per chromosome. ROH, runs of homozygosity.

The ROH segments were further classified based on their length into 1 to 2, 2 to 4, 4 to 8, 8 to 16, and >16 Mb classes (Table 2). Short ROH segments (1 to 4 Mb) were predominant (83.84%). In contrast, long ROH segments (>8 Mb) were few (3.36%). ROH segments of moderate length (4 to 8 Mb) accounted for only 12.75% of all ROH found. There were no marked differences in the percentage of ROH segments in each ROH category across the four generations (Figure 3).

Table 2

Summary of runs of homozygosity segments categorized into different classes based on length

Figure 3

Percentage of ROH segments in each ROH category per generation. ROH, runs of homozygosity.

Genomic inbreeding coefficients

The inbreeding coefficients estimated from the different methods and classes of ROH are summarised in Table 3. The overall average inbreeding coefficient (F_ROH) ranged between 0.039 and 0.327. The other inbreeding coefficients varied between [−0.232 to 3.591], [−0.441 to 0.951], and [−0.337 to 2.271] for F_GRM, F_HOM, and F_UNI respectively. The F_ROH was also calculated for the different ROH length classes. The average F_ROH (1 to 2 Mb) was equal to the average F_ROH values. The F_ROH coefficients estimated from segments longer than 4 Mb were smaller than the overall F_ROH. Figure 4 shows the trends of average inbreeding coefficients across the four generations. There was no marked difference in average F_ROH across the four generations.

Table 3

Summary of inbreeding coefficients calculated using different methods and runs of homozygosity classes

Figure 4

Trend of average inbreeding coefficients across the generations. FROH, FGRM, and FUNI are inbreeding coefficients calculated using methods based on ROH, genomic relation matrix, excess homozygosity, and correlation between uniting gametes, respectively.

The correlations between the inbreeding coefficients were all positive but highly variable (Figure 5). All the F_ROH estimates were strongly correlated (r>0.74) while correlations between the F_ROH and other inbreeding coefficients were low to moderate. The highest correlation was found between F_ROH and F_UNI (r = 0.41) while F_ROH (4 to 8 Mb) and F_GRM had the lowest correlation (r = 0.18). There was a strong correlation between F_GRM and F_UNI (r = 0.81) and between F_HOM and F_UNI (r = 0.88). The correlation between F_GRM and F_HOM was relatively weak (r = 0.41).

Figure 5

Pairwise correlations between estimated inbreeding coefficients in red-brown Korean native chicken population. These included inbreeding coefficients calculated from the diagonal elements of the genomic relationship matrix (FGRM), excess homozygosity (FHOM), correlation of alleles between uniting gametes (FUNI), and runs of homozygosity (FROH).

Runs of homozygosity islands and identification of selection signatures

Common homozygous genomic regions in the study population (ROH islands) were explored by selecting the top 0.1% of the SNPs found in the ROH segments in 44.7% of the population. Consequently, five ROH islands in GGA1, GGA5, GGA7, and GGA8 were identified (Figure 6). The longest ROH was on GGA8 with a length of 2.3 Mb and a total of 114 SNPs, while the shortest ROH was on GGA1 with a length of 0.05 Mb and only four SNPs. The ROH islands had 88 annotated protein-coding genes (Table 4). Functional enrichment of the identified genes did not reveal any significant gene ontology terms probably due to the small number of genes in the ROH islands. Consequently, a literature review was conducted to elucidate the functions of the genes identified. Notable genes identified in the ROH islands included NELL1, BDNF, COL6A1, AMYA1, PLAG4S, and PTGS2. The ROH islands were mapped to the chicken QTL database to identify overlapping QTLs, resulting in 11 enriched QTLs (adjusted p-value<0.05). Of the identified QTLs, 88.25% were associated with important production traits such as body weight, abdominal fat weight, carcass weight and feed conversion ratio (Figure 7A). The top enriched QTLs found per chromosome are depicted in Figure 7B. The ROH island in GGA1 overlapped with QTLs related to body weight and daily weight gain in chickens. GGA5 harbored the highest number of QTLs in the ROH islands which included the QTLs linked to body weight, carcass weight and abdominal fat content. Similarly, the ROH islands in GGA7 and GGA8 contained QTLs related to carcass weight and feed intake.

Figure 6

Manhattan plot for the occurrence of SNPs in ROH islands in the red-brown Korean native chicken population. The red line shows the threshold for selecting the top ROH islands (44.7%). SNPs, single nucleotide polymorphisms; ROH, runs of homozygosity.

Table 4

Identified genes in the runs of homozygosity islands of red-brown Korean native chickens

Figure 7

The summary of QTLs in ROH islands in the red-brown Korean native chicken population (A) Proportion of annotated QTLs for different traits, (B) Bubble plot depicting the top enriched QTLs (adjusted p-value<0.05). QTLs, quantitative trait locus.

Runs of homozygosity

The success of genetic conservation programs relies on the maintenance of genetic variability in a population. The rising demand for livestock products has led to increased selection intensity for higher production, resulting in inbreeding [33]. This can lead to low genetic variability and consequently reduce the population fitness and slow genetic gain in production traits [34]. In the current study, the patterns of ROH in the red-brown Korean native chickens were analysed to assess the extent of inbreeding and genetic variability. Additionally, the ROH islands were evaluated for candidate homozygous regions having genes and QTLs of interest in chicken production.

The length and distribution of ROH segments can give insights into the population’s genetics and breeding history [13]. Homozygous genomic segments inherited from the same ancestor are fragmented into smaller segments by recombination over generations [35]. Thus, short ROH segments indicate ancient inbreeding, while longer ROH segments are due to recent inbreeding.

In this study, short ROH segments less than 4 Mb were prevalent across the four generations, indicating that Korean native chickens may have experienced a bottleneck and inbreeding in the past. The short ROH segments may be due to the impact of other events such as recombination, population contraction and genetic drift. However, it is important to note that the use of medium-density SNP data tends to exaggerate the number of short ROH segments (<4 Mb) [36]. Therefore, the number of short ROH segments in this study might be overestimated. Generally, previous studies including those using whole genome sequences, have reported higher proportions of short ROH segments [19,37]. This suggests that short ROH segments are common in the chicken genome probably because of recombination. On the other hand, the low proportion of long ROH (>8 Mb) indicates that recent inbreeding is uncommon in the population which may suggest good progress in the conservation of genetic diversity. The ROH distribution pattern per chromosome revealed that the number and size of the ROH segments varied with chromosomal size. As expected, the macrochromosomes (GGA1–GGA5) had more and longer ROH segments compared to the smaller chromosomes.

Genomic inbreeding coefficients

Conventionally, the level of inbreeding is estimated using pedigree information [34]. However, inaccurate pedigree records and missing information can limit the accuracy of the estimated inbreeding coefficients. Furthermore, estimates from pedigree records do not account for Mendelian sampling effects. The recent advancement in genomic sequencing technology has enabled the generation of large amounts of genomic data. Thus, genomic inbreeding coefficients can be reliably estimated. Genomic inbreeding coefficients do not depend on pedigree records hence not affected by recording errors. Additionally, these estimates give account to Mendelian sampling effects, making them more realistic [12].

Four types of inbreeding coefficients were calculated and compared in this study. There was no significant change in average F_ROH across the four generations. This further shows that the rate of recent inbreeding in the red-brown Korean chickens is low. The estimated inbreeding coefficients were considerably variable. Unlike other coefficients, the F_ROH estimates ranged between 0 and 1 which is in line with the true definition of the inbreeding coefficient [9]. The average F_ROH calculated from ROH segments shorter than 4 Mb was higher than F_ROH from longer segments. This could be due to the overestimation of short ROH segments [36]. The F_GRM, F_HOM, and F_UNI estimates ranged between negative values and values greater than one. Negative inbreeding coefficients suggest a gain in genetic variability [10]. On the other hand, Inbreeding coefficients greater than one are difficult to interpret because it is not practical to lose more than 100% of the genetic variation [14]. Unlike F_ROH, the three SNP by SNP inbreeding coefficients are dependent on allele frequencies in the population. They only yield reasonable estimates when the allele frequencies in the population are equal to those of the base population [12]. Additionally, these methods do not differentiate genomic regions that are identical by descent from those identical by state [12]. F_ROH could be regarded as a more reliable measure of inbreeding since it is robust to allele frequencies and only reflects the homozygosity that is identical by descent [12].

All the inbreeding coefficients were positively correlated. The overall F_ROH value was strongly correlated with F_ROH (1 to 2 Mb) and F_ROH (2 to 4 Mb) (r>0.98). This suggests that a large proportion of inbreeding in the study population is ancient. This is consistent with previous studies in Chinese indigenous chickens [19]. The correlation between F_ROH and other coefficients was low to moderate. The correlation between F_ROH and F_UNI was higher compared to other inbreeding coefficients (r = 0.41). The low correlation between F_ROH and other inbreeding coefficients could be due to discrepancies in allele frequencies and the different approaches used in the calculation [12,37]. F_GRM was highly correlated with F_HOM. Likewise, the correlation between F_GRM and F_UNI was strong. This is consistent with previous reports on chicken and cattle [16,37].

Runs of homozygosity islands and identification of selection signatures

Selection pressure increases the frequency of homozygous genotypes in a population, thus increasing the homozygosity rate at a given locus and neighbouring sites [38]. Common ROH segments in a population can therefore reveal how selection has modified the variability of specific genomic regions. The current study identified five ROH islands that were common in 44.7% of the study population. Functional annotation of these regions found QTLs associated with production traits, including, body and carcass weight, breast muscle weight, abdominal fat weight and daily weight gain. Tian et al [19] mapped similar QTLs in the ROH islands of indigenous Chinese chickens. These results indicate the loss of genetic variability in genomic regions associated with production traits in chickens probably due to selection.

In this study, annotated genes coding for proteins and long non-coding RNAs were identified in the ROH islands. The ROH islands in GGA5 harboured notable genes including the neural EGFL like 1 (NELL1) and BBOX1, which are linked to the regulation of bone development and feed efficiency in chickens respectively [39–41]. Another important gene located on the ROH island of GGA5 was the leucine-rich repeat-containing G protein-coupled receptor (LGR4) which is involved in the regulation of carbohydrates, lipids and amino acid metabolic pathways [42]. In the current study, the LGR4 gene was mapped in the ROH island overlapping with QTLs associated with growth performance traits. This suggests that the gene may be important in the growth performance of chickens. The same ROH island harboured the brain-derived neurotrophic factor (BDNF) gene which is associated with responses to environmental stressors such as heat stress and light intensity in chickens [43].

The collagen type 11 alpha 1 chain (COL6A1) gene located in ROH island in GGA7 is highly expressed in the growing Graafian follicles of laying chickens [44]. This gene is reported to influence intramuscular fat deposition in chickens; thus, it may affect meat quality [45]. The ROH island in GGA8 contained notable genes, such as amylase alpha 1 (AMYA1) which encodes for the amylase enzyme that is important in starch metabolism [46]. Notably, the same ROH island overlapped with QTL associated with breast muscle weight which suggests that the AMYA1 gene could be important in the growth performance of chickens. This gene was also mapped in the ROH island of Chinese indigenous chickens [19]. Additionally, the ROH island in GGA8 contained the PLA2G4A and PTGS2 genes, which have been linked to reproductive functions in various livestock species [47]. Interestingly, some of the identified genes have been reported as selection signatures in other chicken breeds [16,19]. This suggests that there are common genomic regions that are subject to natural and artificial selection across various breeds of chicken.