Evidence for Genetic Hybridization between Released and Wild Game Birds: Phylogeography and Genetic Structure of Chukar Partridge, Alectoris chukar, in Turkey

: The Chukar Partridge ( Alectoris chukar , Galliformes) is one of the most important game birds in its native range, spanning from the Balkans to eastern Asia, and the regions of Europe, North America and New Zealand where it was introduced. Previous studies found two main genetic lineages of the species forming an eastern and a western clade. Chukar Partridges are raised in game farms and released to supplement natural populations for shooting in the USA, Canada, Greece, and Turkey. To explore intraspeciﬁc genetic structure, phylogeography, and possible genetic admixture events of A. chukar in Turkey, we genotyped individuals from fourteen wild and ﬁve captive populations at two mitochondrial and ten microsatellite DNA loci in. Wild and farmed Chukar Partridge samples were analyzed together to investigate possible inﬂuences of intraspeciﬁc hybridizations. We found that the farmed chukars, which mainly (85%) cluster into the eastern clade, and wild ones were genetically distinct. The latter could be separated into six management units (MUs), with partridges from Gökçeada Island in the Aegean Sea forming the most divergent population. Intraspeciﬁc hybridization was detected between wild and captive populations. This phenomenon causes rampant introgression and homogenization. The phylogeographic analysis revealed admixture among wild populations; nevertheless, this did not impair pointing to Anatolia as likely having a “refugia-within-refugia” structure. We recommend that the genetic structure of Chukar Partridge and its MUs be taken into account when developing the policy of hunting, production, and release to preserve the genetic integrity of this species.


Introduction
Genetic diversity plays a crucial role in the adaptation and survival of species [1], while their phylogeographic structure reflects the complex relationships between historical and ongoing evolutionary processes in a spatial framework. Phylogeographic studies helped elucidate gene flow patterns, hybridization, range expansion, and speciation among many bird species [2].
Climatic fluctuations have been occurring in the last three million years, alternating warm and cold periods in the Northern Hemisphere, especially in mountainous regions, Muscle tissue samples were collected from wild and captive individuals (n = 362) sampled during the 2018-2019 hunting seasons in fourteen localities throughout Turkey and six breeding stations (Table 1). Captive adult individuals were randomly selected in each breeding station. The MAKU-HADYEK-169 protocol controlled all the experiments on Chukar Partridges by MAKU, Local Ethical Committee on Animal Experiments regulations. All samples were preserved at room temperature in absolute ethanol. Total DNA was extracted using GeneJET Genomic DNA Purification Kit (Thermo Fisher Scientific, Waltham, MA, USA) or Dneasy Blood & Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions.  19 16 This study 3 Burdur (BUR) 19 18 This study 4 Eskişehir (ESK) 21 20 This study 5 Çankırı (CKR) 18 18 This study 6 Nigde (NIG) 22 21 This study 7 Sivas (SIV) 20 13 This study 8 Kayramanmaraş (KAH) 29 28 This study 9 Bayburt (BAY) 14 14 This study 10 Erzurum (ERZ) 19 19 This study 11 Kars (KAR) 18 18 This study 12 Bitlis BIT) 13 13 This study 13 Van

DNA Amplification and Sequencing
The partial cytochrome-b (Cyt-b, 1092 bp) and the entire Control region (CR, about 1155 bp) of the mitochondrial DNA (mtDNA) were amplified for all samples following Barbanera et al. [14]. PCR products were purified and sequenced on both strands at Macrogen (Seoul, Korea). Sequences were aligned in GENEIOUS PRIME 2021.2.2 with the MUSCLE plugin and further proofed manually [25]. The sequences were deposited in GenBank with accession numbers MZ706294 to MZ706461. We added 86 partial Cyt-b and CR GenBank sequences of Chukar Partridge from Europe and Asia. The accession numbers of the outgroup Rock Partridge (Alectoris greaca) and GenBank sequences were given in the tree of Supplementary Material S1 ( Table 1).

Phylogeographic Analysis
The Cyt-b and CR sequences were concatenated and aligned. The phylogenetic relationships were reconstructed in MEGA X [28] and BEAST 2 [29] using the Maximum Likelihood (ML) method. The TN93 + G + I algorithm was selected using MEGA X and following the Akaike Information Criterion (AIC). Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Tamura-Nei model and then selecting the topology with a superior log-likelihood value. There were a total of 2247 positions in the final dataset. Genetic differentiation among populations (F ST ) was evaluated by analyzing molecular variance (AMOVA) in ARLEQUIN 3.5 with significances assessed by 10,100 permutations [30] and a spatial analysis of molecular variance (SAMOVA) was performed using SAMOVA 1.0 [31]. This method is based on a simulated annealing procedure aimed at identifying geographically homogeneous populations and maximally differentiated in terms of among-group components (F CT ) of the overall genetic variance without the prior assumption of group composition AMOVA relies on. The program was run for 10,000 iterations from each of 100 random initial conditions and tested all the grouping options (predefined number of groups K ranging from 2 to 18). The optimal number of groups (K): F CT values (proportion of genetic variation among groups) reached a maximum, or a plateau was selected. A median-joining network was created to visualize haplotype relationships using Network 10 [32]. Haplotypes and mismatch distributions of demographic/population expansion were defined by DnaSP 6 [33] and the number of polymorphic sites (S), haplotype diversity (Hd), nucleotide diversity (π), average number of nucleotide differences (K), and number of haplotypes (h) were calculated using DnaSP or ARLEQUIN for the mtDNA dataset.
Microsatellites: Departure from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium were calculated for each microsatellite locus and population with an exact test using GENEPOP 4.7.5 [34]. The mean number of alleles (A), observed (H O ), expected (H E ) heterozygosity, and F ST distances were calculated using ARLEQUIN [30].
We used the Bayesian clustering method implemented in STRUCTURE 2.3.4 [35] to infer the population structure. Ten independent runs with K = 1-10, where K is the different number of subpopulations, were used with an admixture model taking sampling locations as priors and correlated allele frequencies between populations. Throughout the analysis, the burn-in period was fixed at 50,000, and the number of MCMC runs at 20,000. Besides, SAMOVA was performed to identify groups using SAMOVA [31]. The most likely number of groups was determined by 100 repeatedly running with two to 10 groups and choosing those partitions with a maximum F CT value and STRUCTURE HARVESTER [36] according to the method of Evanno et al. [37].

Mitochondrial Nucleotide Sequences
The alignment of concatenated mt DNA loci for 354 individuals had a length of 2247 nucleotides, indels included. A total of 169 haplotypes were found, 146 belonging to the wild populations and 30 to the captive populations. Unique haplotypes (n = 148) were mostly found in the wild populations (n = 139), whereas captive populations yielded only seven haplotypes ( Table 2). The most frequent haplotypes were Hap20 (n = 48) and Hap74 (n = 25) found only in captive partridges except for one wild individual (NIG) at Hap74 and Hap 26 (n = 14) in only wild, respectively. The ML tree with the highest log likelihood (−6412.72) and posterior probability is shown in Figure 1. All haplotypes fell into one of the two main clades, i.e., either the western or the eastern one ( Figure 1 and Supplementary Material). Besides, the median-joining haplotype network showed two main groups where captive and wild populations clearly clustered apart, even if evidences of admixture were also flagged ( Figure 2. When star contraction of 169 haplotypes was applied, 96 haplotypes remained, with partridges from breeding stations and wild populations well separated from each other ( Figure 2). The largest haplogroup included H19, H20, H74, and H88 haplotypes. These were mostly held by partridges from breeding stations (and in 60.5% of all captive individuals as opposed to only 0.9% of wild ones).
Basic summary statistics, including sample sizes, haplotype and nucleotide diversities, are provided in Table 2. We found the highest Hd in KAH and KAR, followed by BAY and central populations, while the lowest was recorded in BSG, BSY, and BSA. Eastern haplotypes were found in the wild CAN, KAH, NIG populations, and all breeding stations. Nearly all samples from breeding stations (85.3%) and some from wild populations (2.1%) belonged to the eastern clade (Table 3).  Figure 1. ML tree built in MEGA X for the aligned haplotypes using the TN93 + G + I model. Th posterior probability of trees in which the associated taxa clustered together is shown next to th branches. Haplotype details are given as Supplementary Documents S1 and S2. Admixed individuals were excluded from this analysis.
Basic summary statistics, including sample sizes, haplotype and nucleotide diversities, are provided in Table 2. We found the highest Hd in KAH and KAR, followed Haplogroup Haplogroup Figure 1. ML tree built in MEGA X for the aligned haplotypes using the TN93 + G + I model. The posterior probability of trees in which the associated taxa clustered together is shown next to the branches. Haplotype details are given as Supplementary Documents S1 and S2. Admixed individuals were excluded from this analysis. other except for BSK (FST = 0.29-0.40, p < 0.001). Some wild populations were not distinguished from the other wild populations (p < 0.001). However, CAN, an island population, KAH, BAY, and ESK were differentiated from the other wild populations (p < 0.001; Figure 3).  Mismatch distribution was unimodal, and population expansion was accepted for wild populations, while for captive populations were multimodal and demographic expansion was not supported (Figure 4).  Table 3. Spatial analysis of molecular variance (SAMOVA) of Alectoris chukar for the mtDNA. K, number of groups. F CT , the proportion of total genetic variance due to the differences between groups. Captive and wild populations' K was given in the first column. Some population differentiation was observed in the wild and captive Chukar Partridges in the SAMOVA for mt DNA. The highest F CT value was found at K = 5 (F CT = 0.56; p < 0.0001). One group included all wild populations, and the breeding stations were divided into four groups (Tables 2 and 3). When only wild populations were analyzed, we did not observe an F CT plateau, rather it increased with the number of K. We selected K = 6 due to the highest F CT differentiation between K = 6 and K = 7 (F CT = 0.11; p < 0.001; Figure 3, Table 3). The first population to split apart was CAN (K = 2), followed by KAH (K = 3).  Mismatch distribution was unimodal, and population expansion was accepted for wild populations, while for captive populations were multimodal and demographic expansion was not supported (Figure 4). The AMOVA revealed that differences among populations accounted for 48.37% of the overall genetic variance observed and differences within populations for 51.63%. The differentiation between wild and captive populations was moderate yet statistically significant (F ST = 0.48, p < 0.01). The differentiation among individuals from only wild or captive populations was low statistically significant (F ST = 0.14 and F ST = 0.24, p < 0.01). The differences among and within wild populations accounted for 13.66% and 86.34%, while in the case of captive populations these figures were 24.71% and 75.29%.
Pairwise F ST estimates ( Figure 3) revealed several well-distinct groups. Breeding station differing by high levels of divergence from wild populations but not from each other except for BSK (F ST = 0.29-0.40, p < 0.001). Some wild populations were not distinguished from the other wild populations (p < 0.001). However, CAN, an island population, KAH, BAY, and ESK were differentiated from the other wild populations (p < 0.001; Figure 3).
Mismatch distribution was unimodal, and population expansion was accepted for wild populations, while for captive populations were multimodal and demographic expansion was not supported (Figure 4).

Microsatellite Analysis
The mean number of alleles per locus varied from 6.6 to 13.2 across wild populations and 6.9 to 8.3 in captive ones (Table 4). A total of 203 alleles were found of which 43 at Aru1E97, followed by 35 at AruF25 and 24 at Aru1B3 (Table 5). While private alleles were found in 11 populations, KAH was the one yielding the highest number (5 alleles; Table  5). The mean expected heterozygosity ranged from 0.69 to 0.87 in wild populations and from 0.71 to 0.79 in captive ones; observed heterozygosity ranged from 0.62 to 0.79 in the former and from 0.64 to 0.73 in the latter, respectively (Table 4).

Microsatellite Analysis
The mean number of alleles per locus varied from 6.6 to 13.2 across wild populations and 6.9 to 8.3 in captive ones (Table 4). A total of 203 alleles were found of which 43 at Aru1E97, followed by 35 at AruF25 and 24 at Aru1B3 (Table 5). While private alleles were found in 11 populations, KAH was the one yielding the highest number (5 alleles; Table 5). The mean expected heterozygosity ranged from 0.69 to 0.87 in wild populations and from 0.71 to 0.79 in captive ones; observed heterozygosity ranged from 0.62 to 0.79 in the former and from 0.64 to 0.73 in the latter, respectively (Table 4).   The linkage equilibrium was rejected for only 22 out of 900 pairs of alleles. However, after sequential Bonferroni correction, exact tests for genotypic linkage disequilibrium was non-significant. These results indicated that the loci used segregate independently. Hardy-Weinberg Equilibrium (HWE) was not accepted for 54 out of 200 the loci in all localities (Table 5). Deviation from HWE was found in all populations except BSA and BSG (Table 4), which might be indicative of inbreeding, assortative mating or null alleles. Heterozygote deficiency appeared in one to 14 loci (Table 5). Heterozygote excess (HE) occurred in NIG at Aru2D020. STRUCTURE HARVESTER indicated that the most likely number of clusters was K = 2 using the log-likelihood (L(K)) concept ( Figure 5). Wild and captive populations separated at K = 2; inferences of K = 3 to K = 6 were similar, revealing three main groups with CAN always well differentiated. Captive samples showed evidence of admixture ( Figure 6). Birds from the CAN island population are separated from the other wild populations at K = 3 to K = 6. By arbitrarily defining an individual as belonging to a specific cluster when assignment probability (q) was above 0.8, 346 of 362 individuals clustered together at K = 3 in Structure ( Table 6). The highest percentages of admixed (less than 80% of the individuals assigned to the cluster) individuals between wild and captive clusters were observed in SIV, VAN, KAR, ESK, ERZ, and KAH (Table 4).  (Table 5). Deviation from HWE was found in all populations except BSA and BSG (Table 4), which might be indicative of inbreeding, assortative mating or null alleles. Heterozygote deficiency appeared in one to 14 loci (Table 5). Heterozygote excess (HE) occurred in NIG at Aru2D020. STRUCTURE HARVESTER indicated that the most likely number of clusters was K = 2 using the log-likelihood (L(K)) concept ( Figure 5). Wild and captive populations separated at K = 2; inferences of K = 3 to K = 6 were similar, revealing three main groups with CAN always well differentiated. Captive samples showed evidence of admixture ( Figure 6). Birds from the CAN island population are separated from the other wild populations at K = 3 to K = 6. By arbitrarily defining an individual as belonging to a specific cluster when assignment probability (q) was above 0.8, 346 of 362 individuals clustered together at K = 3 in Structure ( Table 6). The highest percentages of admixed (less than 80% of the individuals assigned to the cluster) individuals between wild and captive clusters were observed in SIV, VAN, KAR, ESK, ERZ, and KAH (Table 4).   (Table 5). Deviation from HWE was found in all populations except BSA and BSG (Table 4), which might be indicative of inbreeding, assortative mating or null alleles. Heterozygote deficiency appeared in one to 14 loci (Table 5). Heterozygote excess (HE) occurred in NIG at Aru2D020. STRUCTURE HARVESTER indicated that the most likely number of clusters was K = 2 using the log-likelihood (L(K)) concept ( Figure 5). Wild and captive populations separated at K = 2; inferences of K = 3 to K = 6 were similar, revealing three main groups with CAN always well differentiated. Captive samples showed evidence of admixture ( Figure 6). Birds from the CAN island population are separated from the other wild populations at K = 3 to K = 6. By arbitrarily defining an individual as belonging to a specific cluster when assignment probability (q) was above 0.8, 346 of 362 individuals clustered together at K = 3 in Structure ( Table 6). The highest percentages of admixed (less than 80% of the individuals assigned to the cluster) individuals between wild and captive clusters were observed in SIV, VAN, KAR, ESK, ERZ, and KAH (Table 4).   Comparable population differentiation was observed in wild and captive Chukar Partridges in the SAMOVA at ten microsatellite loci, and K = 4 distinguished among CAN, the other wild populations, BSU, and the other breeding stations (F CT = 0.12; p < 0.01; Table 3). When only wild populations were analyzed, we did not observe an F CT plateau; rather it decreased with the number of K ( Table 3). The first population to emerge as distinct was CAN (K = 2), followed by HAK (K = 3), MUG (K = 4), and BIT (K = 5), respectively. Pairwise F ST values based on microsatellites showed significant genetic differentiation among most localities ( Figure 2). Non-significant values were obtained between neighboring east Anatolian localities.

Discussion
The goal of this study was threefold. First, we aimed to determine population structure and phylogeography of Chukar Partridges in Turkey; second, to investigate whether there is any difference between wild and farmed individuals; and third, to search for possible signatures of admixture between them.

Population Genetic Structure
Our mtDNA analyses of Chukar Partridges from Anatolia showed that farmed Chukars are genetically different from wild ones as well as that the two clusters they belong to fall within the western and eastern clade, respectively, emerged in previous studies [23,38] (Figure 1). Concordantly, microsatellites structure showed wild and captive birds to group in two distinct clusters (K = 2). CAN, an island population (Gökçeada Island), emerged as the most genetically differentiated one on the basis of mtDNA and microsatellites (K = 3-6) among the wild populations, and the other follow on the same order they are listed in Table 3. Noteworthy, this genetic picture emerged from F ST and SAMOVA of both genetic systems used (Figures 2, 3 and 5). Even if the wild populations clustered together, it is still possible to detect some internal differentiation among eastern Anatolian, central Anatolian and ESK, KAH, plus BAY populations as well as that evidences of admixture between them occur with the exception of CAN. The differentiation between captive and wild populations was in line with the dissimilar shape and size of bill of their individuals [39]. When wild and captive populations were analyzed separately, it was found that the captive populations were highly differentiated from each other (captive F ST = 0.24, wild F ST = 0.14). This may be due to the fact that bloodlines used at breeding stations are sometimes reinforced with confiscated Chukar Partridges from illegal hunters.
We found some wild individuals falling in the eastern clade, CAN, KAH, and NIG (Table 2). Also, we have determined some genetic admixture between farmed (of eastern origin) and wild individuals in six wild populations and four farms. While 2% of hybrid individuals were found in the wild population, a higher hybridization rate (5%) occurred in farms (Table 6). If this process continues, these admixtures might significantly alter the gene pool of wild populations, possibly impair their fitness and affect female reproduction due to low carotenoid levels in blood plasma (as observed in the Red-Legged Partridges (Alectoris rufa) [48]). A similar genetic homogenization was found in the Mallard, another popular game bird in Europe, with captive-bred individuals changing the gene pool of wild populations [23]. Casas et al. [24] showed that extinction risk of wild and genetically preserved Red-Legged Partridge populations through releases of farmed hybrids is a possibility. Our results show that high haplotype and nucleotide diversity exists in wild Chukar Partridge populations in Turkey as opposed to farm populations (Table 2). Nevertheless, introgressive hybridization might reduce the distinctiveness and diversity of wild populations, impacting their fitness in the near future.
The unimodal mismatch distribution results indicated that all wild populations together experienced recent demographic expansion. However, when taken separately their multimodal mismatch distribution suggests that admixture of haplotypes from three previously isolated lineages (one of them was captive) might have occurred. Also, multimodal mismatch distribution of population separated by SAMOVA at K = 6 indicated previously isolated lineages. Although Anatolia is not covering all the range of Chukar Partridge, these linages may be considered as potential refugia within Anatolia, one of the most important unglaciated areas in Western Palearctic during the Pleistocene. The phylogeographic analysis showed that possibly Anatolia might have been a refugium with "refugia-within-refugia" structure. This model was supported by previous studies indicating range shifts within this region, as in case of Kurper's Nuthatch, Sitta krueperi [49]. A dense forest cover existed in the northern Anatolia and its coastal belts [50], which, according to Albayrak et al. [4], might have hosted three refugia for Kruper's Nuthatch in Last Glacial Maximum (LGM), with the late Quaternary glacial-interglacial cycles shaping subsequent demographic expansion. Overall, it is assumed that many Anatolian species underwent population expansion before the Last Interglacial (LIG) [51,52], or between LGM and LIG [49].

Heterozygosity and Inbreeding
Estimates of observed heterozygosity are significantly lower than expected, except in captive populations, BSG and BSA (Table 4). Widespread heterozygote deficiency (Table 5) might be indicative of a genetic diversity loss in wild and farmed Chukar Partridges. Inbreeding is confirmed, especially in MUG and BUR (indicated by F IS value higher than 0.2). The positive F IS is an indication of decreasing heterozygosity due to null alleles. Similarly, positive F IS was observed in the historical wild group of Mallards in Europe [23].

Taxonomic and Conservation Implications
Fourteen Chukar Partridge subspecies are recognized worldwide, and two of them, A. c. cypriotes and A. c. kurdestanica, occur in southwestern/south-central Turkey, respectively [13]. Our finding supports the two described subspecies (see Figure 3; depicted in blue and green, respectively). Moreover, (SAMOVA and F ST results) might possibly support a new subspecies distributed in Gökçeada Island (CAN).
To preserve the genetic diversity of the Chukar Partridge in Turkey, a country where the release of captive individuals is a common practice, six management units (MUs) should be taken into account: CAN, KAH, BAY, ESK, south-eastern Turkey, and central Anatolia separately (Figure 3). Specific conservation efforts should be made for the population of Gökçeada Island (CAN), where partridge releases should be banned due to the occurrence of a genetically distinct and well-preserved A.chukar population. We advocate for the genetic identity of Chukar Partridges and their six MUs to be considered when developing hunting, production, and releasing policies to preserve the integrity and internal diversity in perpetuity.

Data Availability Statement:
The data presented in this study are available in GenBank with accession numbers MZ706294 to MZ706461.