The experimental protocol for this study was approved by the Animal Experiment Committee of Kasetsart University (approval numbers: ACKU65-SCI-021 and ACKU65-SCI-023) and was conducted in accordance with the Regulations on Animal Experiments at Kasetsart University and the ARRIVE guidelines (https://arriveguidelines.org).

Whole blood samples were collected from five populations of PDH chickens and two populations of SD chickens. The data of one population of PDH chickens (Phitsanulok 1: PDH1) was obtained from our previous study [8] (Supplementary Table S1). The PDH2 population was derived from Phitsanulok 2, PDH3 was from Chiang Mai, PDH4 was from Nakhon Pathom, and PDH5 was from Nonthaburi. SD1 population of SD chickens was derived from the Department of Livestock, Uthai Thani population, and SD2 was derived from Sanhawat Farm, Uthai Thani. Specimen collection, DNA extraction, and DNA qualification and quantification were performed as described previously [10].

Twenty-eight microsatellite primer sets were selected based on the recommendations of the Food and Agriculture Organization for assessing the biodiversity of chicken populations (Supplementary Table S2). The 5′-end of the forward primer of each primer set was labeled with fluorescein amidite (6-FAM) or hexachloro-fluorescein (HEX) (Macrogen Inc., Seoul, Korea). Microsatellite PCR amplification was performed using a 15 μL reaction mixture containing approximately 1× Apsalagen PCR buffer (Apsalagen Co., Ltd., Bangkok, Thailand), 25 ng genomic DNA, 1.5 mM MgCl₂, 0.2 mM dNTPs, 0.5 μM of each primer, and 0.5 U Taq polymerase (Apsalagen, Thailand). The PCR conditions were as follows: initial denaturation at 95°C for 10 min, followed by 35 cycles at 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s, and a final extension for 5 min at 72°C, as described previously [10]. To avoid false allele amplification, each individual sample was run in triplicate.

Genetic diversity parameters, such as allelic richness (AR) were calculated in FSTAT version 2.9.4 [14], polymorphic information content (PIC) was calculated in Excell Microsatellite Toolkit [15], number of alleles per population (N_a), observed and expected heterozygosity (H_o and H_e), F-statistics (F_IS and F_ST), relatedness (r), and pairwise Nei’s genetic distance were calculated in GenAlEx version 6.5 [16]. Population structure analyses, such as analysis of molecular variance (AMOVA) and principal coordinate analysis (PCoA) were estimated in GenAlEx version 6.5 [16], discriminant analysis of principal components (DAPC) was performed using package ADEGENET 2.0 [17] in R version 4.1.2 [18], and STRUCTURE analysis was performed in STRUCTURE version 2.3 [19] following Budi et al [10]. The Wilcoxon signed-rank test for detecting recent population bottlenecks was performed using a two-phased mutation model (TPM) and stepwise mutation model (SMM) to assess the probability of excess heterozygosity due to small sample sizes. Tests for relative long-term bottleneck events based on M ratio were performed following Budi et al [10]. To investigate the genetic selective sweeps in PDH and SD chickens, a selective sweep analysis was performed as described previously [10]. H_e and F_IS values of populations were plotted for each microsatellite locus (28 loci in total). High F_IS and low H_e values reflect sweeping or purifying selection signatures, whereas low F_IS and high H_e values indicates neutral or balanced selection.

To examine the occurrence of genetic exchange between PDH and SD chicken populations based on the microsatellite genotyping dataset, Bayesian interference analysis implemented in BayesAss version 3.0.5 [20] which is commonly used to determine recent migration rates between populations was performed. The Markov Chain Monte Carlo (MCMC) analysis was conducted for 10 million generations after a burn in a period of 1 million generations and sampled every 100 generations. The mixing parameters associated with migration rates (m), allele frequencies (a), and inbreeding coefficients (f) were optimized to satisfy 20% to 60% posterior distribution acceptance rates according to the recommended guideline [20]. Similarly, to assess the historical genetic exchange between PDH and SD chicken populations, migration rates between the populations and their effective population sizes were estimated by Bayesian analysis using MIGRATE-N version 4.4.3 [21,22]. In these analyses, uniform prior distributions were used for the basic microsatellite model, and 5,000 steps were recorded every 100 generations using the MCMC procedure. The first 100,000 generations were discarded as burn-ins. Estimates were calculated for the mutation-scaled immigration rate (M) and mutation-scaled population size (Θ). The number of individuals entering populations (N_m) was calculated, and the presence of gene flow between populations in the past was determined using the formula N_mi–>j = Θ_j*M_i–>j/4, where N_mi–>j represents the effective number of immigrants or gene flow rate from population i to population j per generation. Circos (version 0.69–8) was used to visualize the genetic connectivity among populations [23]. The genotypic data generated in this study have been stored in the Dryad digital repository dataset (https://doi.org/10.5061/dryad.hhmgqnkm0, accessed on 15 January 2024).

The genetic origins of PDH and SD chickens were investigated using the microsatellite genotyping data of chickens available under the SCBP (https://www.sci.ku.ac.th/scbp/; https://datadryad.org/stash/share/x2qlPmboMgCROXO8, accessed on January 15, 2024), including red junglefowl and indigenous and local chicken breeds in Thailand. Two populations of SD chickens and five populations of PDH chickens were treated as separate data sets and compared the data with those of other indigenous and local chicken breeds, as well as with red junglefowl, which was considered a different data set. Pairwise genetic distances between populations and clustering analyses using PCoA, DAPC, and STRUCTURE were performed as described previously [10].

A total of 370 alleles were observed in the five populations of PDH chickens, with the mean number of alleles per locus being 4.086±0.312, whereas 179 alleles were found in the two populations of SD chickens, whose mean number of alleles per locus was 3.607±0.349 (Supplementary Table S3). All allelic frequencies showed a significant departure from the Hardy–Weinberg equilibrium of the population in both PDH and SD chicken populations, with evidence for linkage disequilibrium (Supplementary Tables S4–S10). Null alleles were frequently found for the MCW0034, MCW0014, MCW 0067, ADL0112, and MCW0069 loci; nevertheless, all the markers were used with no distinction for estimating the genetic diversity. Three populations of PDH chickens (PDH3, PDH4, and PDH5) and two populations of SD chickens each exhibited negative F values, whereas the other two populations of PDH (PDH1 and PDH2), exhibited positive F values. The PIC of all PDH and SD chicken populations ranged from 0.137 to 0.652, whereas the Shannon’s information index (I) was from 0.061 to 1.442. The mean H_o and H_e values were 0.591±0.052 and 0.594±0.030, respectively, for PDH chickens, and 0.668±0.050 and 0.503±0.036, respectively, for SD chickens. The mean AR value of PDH chickens and SD chickens were 3.771±1.471 and 2.526±0.873, respectively. The standard genetic diversity indices are summarized in Table 1. All pairwise differences in H_o and H_e values between PDH and SD chicken populations were not significant, except for the PDH3 and SD1 populations (Supplementary Table S11).

The mean pairwise r values in PDH and SD chickens were −0.013±0.070 and −0.022±0.106, respectively, whereas the F_IS values were −0.001±0.029 and −0.374±0.060, respectively (Supplementary Tables S12 to S20). AMOVA revealed that genetic variation accounted for 0% of the total variance within the populations of PDH chickens, whereas it accounted for 8% within the populations of SD chickens, and 24% among the populations of PDH chickens, and 29% among the populations of SD chickens, (Supplementary Table S21). Nei’s genetic distances values were ranged from 0.301 to 2.549 among the five populations of PDH chickens, and 1.236 between the two populations of SD chickens (Supplementary Table S22). Pairwise differences (F_ST) based on microsatellite data ranged from 0.043 to 0.440 among populations in PDH and SD chickens (Supplementary Table S23). The Wilcoxon signed-rank tests generated the values of 0.582 and 0.973 under SMM and a TPM for PDH1 (normal L-shaped distribution), 0.493 and 0.537 for PDH3 (normal L-shaped distribution), and 0.695 and 0.327 for PDH4 (normal L-shaped distribution), 0.711 and 0.779 for SD1 (normal L-shaped distribution), and 0.425 and 0.248 for SD2 (shifted mode), respectively. The Wilcoxon signed-rank test could not be calculated for PDH2 and PDH5 due to their small sample sizes. The M ratios of all populations were below 0.68, reflecting a historical decline of populations [24] (Table 1).

PCoA and DAPC classified five PDH chicken populations into four clusters, whereas two SD chicken populations were classified separately into two clusters (Figure 1; Supplementary Figure S1). A Bayesian clustering analysis with STRUCTURE showed that PDH and SD chickens exhibited different population structure patterns at various K-values (K = 2 to 25) (Figure 2; Supplementary Figure S2). Based on Evano’s ΔK values, the optimized structure pattern was assigned to six clusters at K = 6, whereas based on ln P(K), the optimized structure pattern found at K = 7. At K = 6, all populations exhibited different gene pool patterns, except for PDH4 and PDH5. Similarly, at K = 7 or more than 7, all populations, except for PDH4 and PDH5 showed different gene pool patterns. Genetic selective sweep analysis revealed neutral or balanced selection occurred in all populations of PDH and SD chickens, which was reflected by a relatively low F_IS coupled with a high H_e (Supplementary Figure S3). Recent genetic gene exchange estimated using BayesAss revealed that migration rates range were 0.699 to 0.923 within populations and 0.0123 to 0.134 among populations in PDH and SD chickens. The asymmetric mutation rate (M) derived from MIGRATE-N ranged from 7.000 to 777.670, with the highest value observed from PDH3 to PDH2. The median scaled mutation rates (Θ) ranged from 0 (PDH4) to 0.099 (PDH1) (Supplementary Table S25). A diverse range of N_m values (0 to 0.501) was observed in PDH and SD populations (Supplementary Table S26).

Multiple population clusters were observed based on the PCoA and DAPC results (Supplementary Figure S4–S5). The PDH and SD chickens tended to cluster together with indigenous and local chicken breeds and red junglefowl. The gene pool patterns of PDH and SD chickens were compared with those of the indigenous and local chicken breeds and red junglefowl in Thailand using the dataset of microsatellite genotypes in our previous studies (SCBP, https://www.sci.ku.ac.th/scbp/). STRUCTURE analysis revealed the highest posterior probability with one peak at K = 2 based on Evano’s ΔK, whereas the mean ln P(K) showed a different peak at K = 25 (Supplementary Figure S6 to S7). At K = 2, the two populations of SD chickens shared the similar gene pool with all PDH populations, except for PDH1 population. The peak at a high K value (K = 25) indicated the presence of multiple clusters in SD and PDH populations. The SD2 population shared a partial gene pool with red junglefowl (G. gallus gallus) derived from Roi Et and Sisaket, whereas the SD1 population partially shared the gene pool with those of the PDH3, PDH1, and Lueng Hang Khao, Chee, Khaew Paree, and Decoy population (Supplementary Figure S6). The PDH4 and PDH5 populations shared similar gene pools with those of the Lao Pa Koi, Prama, Wein Chang, Chee Fah (Mae Hong Son), Fah Luang (Mae Hong Son), Mae Hong Son, and Dong Tao (Udon Thani) population. The PDH2 population showed a distinctive gene pool compared to other chicken populations. Additionally, a small portion of the gene pool of PDH2 population was observed in SD2 and PDH3 populations. The gene pools of SD and PDH populations were closely related to those of other indigenous and local chicken breeds and red junglefowl in Thailand, and no detectable genetic selective sweeps were observed in any chicken breeds and red junglefowl.

Population structure analysis using STRUCTURE revealed that PDH and SD chicken populations are separated into several different clusters at K = 6 and K = 7, which was consistent with the results by PCoA and DAPC analyses. The SD1 chicken population exhibited a distinct gene pool that was different from those of the SD2 chicken population and PDH chickens through all K-values, which suggests that the two SD chicken populations have distinct genetic origins [10]. The SD1 population exhibited a strikingly gene pool pattern similar to those of PDH1 and PDH3 populations. Genetic exchange between the SD1 population and PDH chickens was confirmed using gene flow analysis based on Bayesian interference analysis of recent migration rates between populations (Supplementary Figure S8) [29]. This suggests that PDH and SD chickens may have been originated from the same breed known as the PDH chicken breed. PDH and SD chickens have similar body sizes and weight; however, they can be distinguished by the coloration of their eyes, face, combs, and hackles. PDH chickens have red eyes, faces, and combs, and red-brown colored hackles, whereas SD chickens have typically black eyes, faces, and combs, and brownish-oak colored hackles [4,5,7]. Therefore, SD chickens have been considered to be a variety of the PDH chicken breed. Long-term selective breeding. focusing on the production of high quality of poultry products and ease of breeding, has led to the improvement of the PDH chicken breed. The sub-group that differs in either external morphologies or plumage color or both, within a breed is designated as a “variety.” These variations are found in PDH chickens, such as the Pradu Hang Dam Chiang Mai and Pradu Hang Dam KKU55 populations [30,31]. Recent selection biases have resulted in the homogenization of morphological features of the PDH chicken breed. This led to the rapid reduction in genetic diversity and acceleration of genetic similarity between varieties as observed in SD chickens.

Notably, the SD1 and SD2 chickens were probably derived from different origins and evolved into distinct subgroups as different sub-varieties. The sub-group with specific desired traits improved by selective breeding, is referred to as “strain” [32]. The SD1 chicken population, which was raised for meat production, was genetically related to PDH chickens studied here, which were improved for meat production and cockfighting. Their gene pools were partially shared with that of the SD2 chicken population, which was improved for ornamental competition with red junglefowl derived from Roi Et and Si Sa Ket, which is why there is a genetic footprint from red junglefowl in SD chickens. Genetic variability in the variety of SD chickens may be influenced by the breeders’ purposes, shaping the selection and domestication processes. However, the microsatellite DNA markers used in this study could not identify the genetic traits linked to them in each strain or population. High-throughput genotyping, such as double digest restriction-site associated (ddRAD) or whole-genome sequencing, would be alternative options for the in-depth identification of genes or genetic markers associated with the chicken breed- or variety-specific biological features.

Many chicken breeds have been developed for entertainment because of their specific morphological, external and behavioral features, such as body shapes and sizes, combs, wattles, tail feathers, plumage colors, crowing sounds, and rooster’s fighting behavior. However, according to the annual report of the Food and Agricultural Organization (FAO), 21% of 1,729 chicken breeds are on the verge of extinction [33]. These alarming statistics present the critical need to identify the factors that drive genetic erosion and extinction risks in chicken breeds, including the indigenous and local breeds and to take appropriate actions. Loss of genetic diversity within and between breeds, that is, genetic homogenization of chicken breeds, can occur due to various factors, such as economic pressure to make profit, ecological changes, social shifts, and government policies. To conserve the genetic diversity of chicken breeds effectively, determining these driving factors by population genetics analyses of genetic diversity, population structures, and population movement monitoring and taking appropriate actions are required. In the SCBP at Kasetsart University, the efforts to protect endangered breeds are enforced by collecting and preserving DNA resources of indigenous and local chicken breeds and red junglefowl in collaboration with local communities in Thailand. An alternative approach would be to emphasize the cultural and economic value of indigenous and local chicken breeds in local communities and encourage their active involvement in conservation. The SD chicken varieties, which have unique phenotypes and cultural significance are on the brink of extinction. Therefore, preserving them would not only ensure a stable supply to local communities but also potentially contribute to maintaining sustainable agriculture through the production of indigenous and local chicken breeds. Genetic resources of indigenous and local chicken breeds have the potential for meeting customers’ demands of breeding in the future, which would generate profits for both farmers and consumers. Employment generation and income distribution in local communities can be facilitated by implementing the strategies of effective conservation of useful chicken resources to support sustainable agriculture.