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

Objective The analysis of runs of homozygosity (ROH) has been applied to assess the level of inbreeding and identify selection signatures in various livestock species. The objectives of this study were to characterize the ROH pattern, estimate the rate of inbreeding, and identify signatures of selection in the red-brown Korean native chickens. Methods The Illumina 60K single nucleotide polymorphism chip data of 651 chickens was used in the analysis. Runs of homozygosity were analysed using the PLINK v1.9 software. Inbreeding coefficients were estimated using the GCTA software and their correlations were examined. Genomic regions with high levels of ROH were explored to identify selection signatures. Results A total of 32,176 ROH segments were detected in this study. The majority of the ROH segments were shorter than 4 Mb. The average ROH inbreeding coefficients (FROH) varied with the length of ROH segments. The means of inbreeding coefficients calculated from different methods were also variable. The correlations between different inbreeding coefficients were positive and highly variable (r = 0.18–1). Five ROH islands harbouring important quantitative trait loci were identified. Conclusion This study assessed the level of inbreeding and patterns of homozygosity in Red-brown native Korean chickens. The results of this study suggest that the level of recent inbreeding is low which indicates substantial progress in the conservation of red-brown Korean native chickens. Additionally, Candidate genomic regions associated with important production traits were detected in homozygous regions.


INTRODUCTION
Indigenous chickens are important genetic resources due to their adaptive capacity and disease resilience.Thus, their conservation is essential for the improvement of poultry breeds.However, the existence of these indigenous breeds is threatened by the rapid introduction of high-producing commercial chicken breeds [1].The Korean native chickens have inhabited the Korean Peninsula for many years, becoming a remarkable genetic resource in the region [2].The chickens underwent a severe bottleneck during the Korean War, threatening their existence.Since then, restoration and conservation programs have been implemented [3].
The red-brown chicken breed is among the five lines of Korean native chickens which are distinguishable by their feather colour [4].Despite their slow growth rate, Korean native chickens have recently gained popularity due to their superior meat quality traits.Thus, they are potential resources in meeting the expanding demand for poultry products [2].
The red-brown line is reported to have favourable carcass quality traits and has been crossbred with exotic broiler chickens to improve meat production [4,5].
Genetic conservation programs should preserve alleles and allelic combinations and maximise heterozygosity in populations [4,6].The progress in genetic conservation is assessed through the estimation of genetic diversity parameters [7].The inbreeding coefficient is widely used in the estimation of genetic variability in livestock populations [8].It is defined as the probability that two alleles at a locus are inherited from a common ancestor [9].
Conventionally, inbreeding coefficients in livestock populations have been calculated using pedigree data.However, these estimates are often unreliable due to the scarcity of pedigree information and inaccurate records [10].The impreciseness of pedigree-based inbreeding estimation and the availability of large amounts of genomic data have inspired the application of genomic-based estimation of inbreeding coefficients [10,11].Genomic inbreeding coefficients are more reliable than pedigree-based coefficients because they capture ancient inbreeding and account for the random inheritance of alleles [11].Additionally, the level of genetic variability in specific genomic regions can be revealed.
Several approaches have been applied in the estimation of genomic inbreeding coefficients.These include methods based on runs of homozygosity (ROH), excess homozygosity, diagonal elements of genomic relationship matrices and correlation between alleles in uniting gametes [10].ROHbased inbreeding coefficients are widely used because they tend to have a strong correlation with pedigree-based inbreeding coefficients and reflect the homozygosity level due to inbreeding [12].
ROHs are long homozygous DNA segments that form due to the inheritance of identical haplotypes from the same ancestor [13].Their characterization is important in revealing the breeding history of a population.The analysis of ROH has been conducted on various livestock species, including chickens [14][15][16][17][18]. Indeed, the impact of inbreeding on important livestock traits has been demonstrated through studies on ROH.For instance, the negative effects of inbreeding on production and fertility traits in dairy cattle have been investigated through studies on ROH [15][16][17][18].
ROH analyses have also been employed to explore the impact of selection on the genetic diversity of animals.In chickens, investigations of the ROH have been conducted to assess the patterns of inbreeding and identify selection signatures in various native and commercial lines [16,17].Using high-density single nucleotide polymorphisms (SNPs) markers, the patterns of ROH in Chinese indigenous chickens were characterized and candidate genes and quantitative trait locus (QTL) related to body weight, feed intake and muscle development on the ROH islands were identified [19].
The genetic diversity of Korean indigenous chickens has been examined using microsatellite and SNP markers [7,20].Previous studies have mainly focused on the characterization of population structure, linkage disequilibrium decay and phylogenetics.These parameters reflect the Euclidean distance between individuals in a population [21].Estimation of relatedness and inbreeding coefficients in Korean native chickens using genomic data has not been comprehensively conducted.In this regard, this study aimed to characterise the ROH, estimate inbreeding coefficients, and identify candidate genes and QTLs in homozygous genomic regions in the red-brown Korean native chicken.

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 -6 ) 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.

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.measurements of inbreeding were estimated using the (--ibc) function in 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: hreshold of overlapping homozygous windows was set to 0.05.The the sliding window was estimated based on the formula proposed by easurements of inbreeding were estimated using the (--ibc) function in e include: 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: Inbreeding coefficient calculated from the diagonals of the genomic kinship matrix (FGRM) [29] as 155 follows: Where Z is a covariance matrix constructed by subtracting twice the minor allele frequency from 160 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: Inbreeding coefficient calculated from the diagonals of the genomic kinship matrix (FGRM) [29] as Wright's inbreeding coefficients calculated as the correlation between alleles in uniting gametes 168 (FUNI) [28]: Where x is the number of versions of the reference allele, p and q are the reference and alternative 173 allele frequencies of m th SNP respectively.The correlation between FROH coefficients and other 174 inbreeding coefficients was estimated using Pearson's correlation and plotted using the Corrplot 175 package in R software [26].The trend in average inbreeding coefficients across the four generations 176 was also evaluated.Wright's inbreeding coefficients calculated as the correlation between alleles in uniting gametes (F UNI ) [28]: Where Z is a covariance matrix constructed by subtracting twice the minor allele frequenc 160 the raw marker score.The covariance matrix consists of values: (0-2p) for homozygous, ( Wright's inbreeding coefficients calculated as the correlation between alleles in uniting g 168 (FUNI) [28]: Where x is the number of versions of the reference allele, p and q are the reference and alte 173 allele frequencies of m th SNP respectively.The correlation between FROH coefficients and 174 inbreeding coefficients was estimated using Pearson's correlation and plotted using the C 175 package in R software [26].The trend in average inbreeding coefficients across the four gene 176 was also evaluated.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).
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).

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 be-

(A) (B)
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.
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).

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.

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, popu-lation 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][40][41].Another important gene located on the ROH island of GGA5 was the leucinerich 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.

7
a covariance matrix constructed by subtracting twice the minor allele frequency from 160 the raw marker score.The covariance matrix consists of values: (0-2p) for homozygous, (1-2p) for 161 heterozygous and (2-2p) for alternate homozygous loci, where p is the minor allele frequency.162 The Wright's inbreeding coefficient (FHOM) which compares the expected and the observed 163 homozygous genotypes for each sample [28] as follows: 164 165  ��� �    ℎ  �    ℎ      �    ℎ  166 167 annotation of genes and QTL in ROH islands 179 annotation of genes and QTL in ROH islands 179

Figure 1 .
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 3 .Figure 3 .
Figure 3. Percentage of ROH segments in each ROH category per generation.ROH, runs of homozygosity.

Figure 2 .
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.

Figure 2 .
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.

Figure 4 .
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.

Figure 5 .
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).

Figure 5 .
Figure 5. Pairwise correlations between estimated inbreeding coefficients in red-brown

Figure 6 .
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.

28 612Figure 6 .
Figure 6.Manhattan plot for the occurrence of SNPs in ROH islands in the red-brown

Figure 7 . 29 625Figure 7 .
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.

Table 1 .
Number of samples per generation 1-2   161heterozygous and (2-2p) for alternate homozygous loci, where p is the minor allele frequency.

Table 2 .
Summary of runs of homozygosity segments categorized into different classes based on length Macharia et al (2024) Anim Biosci 37:1355-1366 tween [-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

Table 3 .
Summary of inbreeding coefficients calculated using different methods and runs of homozygosity classes ROH, runs of homozygosity; F GRM , genomic relationship matrix; F HOM , excess homozygosity; F UNI , the correlation of alleles between uniting gametes; F ROH , inbreeding coefficient from ROH.

Table 4 .
Identified genes in the runs of homozygosity islands of red-brown Korean