Molecular evidence of introgressive hybridization between related species Jankowski's Bunting (Emberiza jankowskii) and Meadow Bunting (Emberiza cioides) (Aves: Passeriformes)

1.   Introduction

Natural hybridization often occurs between closely related species exhibiting sympatric or parapatric distributions (Moore, 1998; Abbott et al., 2013). Introgressive hybridization in particular plays an important role in animal evolution; at least 10% of animal species in nature are involved in potential introgressive hybridization with other species (Mallet, 2005). Introgressive hybridization often leads to negative results, such as the loss of species-specific alleles, which tends to make different species more genetically homogeneous (Jacquemyn et al., 2012; Brower, 2013; Charles and Stehlik, 2021). However, hybridization is also an important source of genetic variation within populations (Barton, 2001) and can produce new gene combinations adapted to different niches (Mallet, 2007; Araya-Anchetta et al., 2013). When the F1 generation resulting from crosses between incompletely differentiated species are backcrossed to the parents, gene introgression may occur, resulting in the integration of alleles from one species into the gene pool of another (Arnold, 1997). Interspecific hybridization can thus supplement the impoverished gene pools of endangered species, which is beneficial to their continued survival.

Jankowski's Bunting (Emberiza jankowskii) and Meadow Bunting (E. cioides) are two closely related Old World buntings (Päckert et al., 2015). E. cioides is widely distributed in eastern and central Asia; two subspecies, E. cioides cioides and E. cioides weigoldi, are found in Northeast China (Zheng, 2005). E. jankowskii has been listed as 'Endangered' on the IUCN Red List of Threatened Species since 2010 because of its narrow distribution and rapid population decline over the past five decades (Jiang et al., 2008; Wang et al., 2010). E. jankowskii is currently distributed mainly in eastern Inner Mongolia (Han et al., 2017). The two species are sympatric in this area, similar in morphology, and occupy overlapping niches. In phylogenetic studies, the two species always cluster in the same branch with high confidence (Liang et al., 2008; Päckert et al., 2015; Cai et al., 2021). A previous report revealed that the existence of identical or similar alleles that verified trans-species polymorphisms (TSPs) between the two species (Li et al., 2017). These findings suggest the possibility of hybridization between the two species. We speculate that gene introgression from E. cioides may have contributed to the maintenance of unexpectedly high genetic diversity in endangered E. jankowskii (Li et al., 2021). Therefore, it is necessary to investigate possible introgressive hybridization between these two species.

The core issue of hybridization research is how to accurately identify hybrids from a mixed natural population in the wild. Initial generations of hybrids are usually mixed with their parents to form a hybrid zone (Nolte and Tautz, 2010). It has been found in some studies that hybridization can form new morphological phenotypes or extremely high hybrid variation (Ait Belkacem et al., 2016; Wolfgramm et al., 2021). Nevertheless, it is difficult to accurately identify hybrid individuals by morphological traits, as phenotypes are easily affected by the environment (Rieseberg et al., 1999). Molecular methods, which can sensitively identify genome-wide micro-variations and distinguish hybrids, have been widely used in interspecies hybridization studies (Anderson, 2008; Wang et al., 2014). However, one of the major challenges for molecular methods is to distinguish hybridization from incomplete lineage sorting (ILS) because they leave similar imprints on the genome (Holder et al., 2001; Qu et al., 2012). ILS leads to the retention of common ancestral polymorphisms in differentiated related species, which can seriously interfere with phylogenetic inference and hybrid identification (Kutschera et al., 2014; Wang et al., 2014; Scornavacca and Galtier, 2017).

The application of highly variable co-dominant markers (such as microsatellites) and powerful Bayesian clustering methods facilitates the detection of hybridization. Different types of molecular markers have different degrees of variability, which can reflect introgressive hybridization at different spatial and temporal scales. Highly variable nuclear microsatellites combined with mitochondrial DNA markers can be used to describe past and ongoing gene introgression (Lexer et al., 2005; Wang et al., 2014). Bayesian methods are usually implemented as Markov Chain Monte Carlo (MCMC) algorithms using multi-locus data. This method does not rely on prior taxonomic information, but calculates the posterior probability (PP) of population clusters through various models and hypotheses to distinguish different hybridization statuses (Anderson, 2008). Robust algorithms offer the possibility for coalescent simulation and deep mining and utilization of genetic information.

In the present study, we first used mitochondrial NADH dehydrogenase subunit 2 (ND2) gene and 15 nuclear microsatellite markers to compare the genetic diversity of E. jankowskii and E. cioides, and then detected TSPs between the two species by phylogenetic reconstruction and Bayesian cluster analysis. Having determined the hybridization status, we then examined the direction and intensity of gene introgression, as well as the impact on genetic variation within the population. This study will deepen our understanding of the relationship between introgressive hybridization, genetic diversity, and conservation of endangered species.

2.   Materials and methods

2.1   Sampling and DNA extraction

Fieldwork was carried out mainly during the breeding season from May to July. Mist nets were used to capture adults during the nestling stage (after hatching). All nestlings in the same nest and captured adults were recorded and collected. From 2012 to 2019, blood samples of 390 field-caught birds (328 E. jankowskii individuals and 62 E. cioides individuals, including nestlings and adults) were collected by wing vein puncture, a minimally invasive sampling method, from five locations representing potential hybrid zones in eastern Inner Mongolia (Fig. 1; Table 1). Details of the individuals sampled and DNA samples used for ND2 and microsatellite analyses are shown in Table 1. All sampling sites were located in the breeding areas of E. jankowskii, which is a subset of the distribution area of E. cioides (Fig. 1). E. cioides was sampled only at site Gahaitu. All adults were identified as one of the two species based on phenotype traits, and no intermediate phenotype were observed. For nestlings, the species was identified after confirming that parents were active around the nest or the female was hatching in the nest. All samples collected in the field were temporarily stored in liquid nitrogen and transferred to the laboratory at −80 ​℃. Genomic DNA was extracted using an AxyPrep Blood Genomic DNA Miniprep Kit (Aisijin, Hangzhou, China) following the manufacturer's instructions.

Figure  1.  Locations of origin of the samples analyzed in this study. The yellow and blue color mean the current worldwide distribution of Emberiza jankowskii and Emberiza cioides, respectively. The sampling sites are shown as pentagrams. The size of the circle is proportional to the sample size (number) at each site (see color legend). E. cioides was sampled only at GHT. ALH: Aolunhua, CG: Chagantaoligao, GHT: Gahaitu, KD: Kundu, LB: Lubei. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

Table  1.  Details of the individuals sampled and DNA samples used for ND2 and microsatellite analyses.

Blood samples		Number of sampled individuals	ND2 samples	Microsatellite samples
E.jankowskii
	Adult	176	176	137
	Nestling	152	55	152
	Total	328	231	289
E.cioidesi
	Adult	43	36	43
	Nestling	19	18	19
	Total	62	54	62

 | Show Table

DownLoad: CSV

2.2   PCR, ND2 sequencing, and microsatellite genotyping

Because the mitochondrial gene ND2 follows maternal inheritance and individuals from the same nest share the same haplotype, we only selected blood samples from unrelated individuals in each nest for subsequent analysis to avoid underestimation of genetic diversity. A total of 212 samples were used for ND2 sequencing, including 176 E. jankowskii individuals and 36 E. cioides individuals (Table 1). We used the primers L5215 (5′-TATCGGGCCCATACCCCGAAAAT-3′) and 1064 (5′-CTTTGAAGGCCTTCGGTTTA-3′) (Irwin et al., 2009) to amplify approximately 1047 base pairs (bp) of the ND2 gene. PCR amplifications were performed in 25-μL volumes containing 1 ​μL template DNA, 1 ​μL of each primer, 12.5 ​μL of 2 ​× ​Taq PCR Master Mix (Sangon Biotech, Shanghai, China) and 9.5 ​μL H₂O. The PCR cycling program was as follows: an initial denaturation step at 94 ​℃ for 3 ​min, followed by 35 cycles of denaturation at 94 ​℃ for 30 ​s, annealing at 55 ​℃ for 30 ​s, extension at 72 ​℃ for 1 ​min, and a final elongation step at 72 ​℃ for 10 ​min, before being held at 4 ​℃. The PCR products were subsequently purified and sequenced by Sangon Biotech (Shanghai, China). We then conducted a BLAST search in GenBank to determine the identity of raw sequences.

Fifteen microsatellite loci and 351 samples (289 E. jankowskii individuals and 62 E. cioides individuals) were used in this study (Table 1). The specific information of each microsatellite primer is listed in Appendix Table S1. Amplification of microsatellite loci followed the procedure previously described for mtDNA amplification, and we optimized the reaction volumes and PCR protocols for each different primer pair. The PCR products were electrophoresed in 1% agarose gels. Final fragment analysis and scoring were performed with GeneMarker v1.91 (SoftGenetics LLC, State College, PA, USA) by Sangon Biotech (Shanghai, China).

2.3   Descriptive statistics and admixture analysis

Genetic diversity was evaluated for each population separately. Mitochondrial ND2 gene sequences were edited using BioEdit Sequence Alignment Editor v. 7.2.5 (Hall, 1999), and then aligned using ClustalX v. 1.83 (Jeanmougin et al., 1998). Estimations of genetic diversity indices including number of haplotypes (H), average number of nucleotide differences (k), haplotype diversity (H_d), and nucleotide diversity (π) were calculated applying DnaSP v. 6 (Rozas et al., 2017). We used PopART v. 1.7 (Bandelt et al., 1999) to construct a median joining (MJ) network based on haplotypes identified previously. MJ, which combines the minimum spanning network and quasi-median network algorithms to reconstruct phylogenies from recombination-free intraspecific data such as mitochondrial DNA variation, is increasingly employed in phylogeographical analysis (Zhai et al., 2019; Frankowski et al., 2020).

For microsatellite analysis, 180 unrelated samples were used for the calculation of genetic diversity, including 137 E. jankowskii samples and 43 E. cioides samples (Table 1). Subsequent population differentiation and admixture analyses used all 351 microsatelite samples to find hybrids as much as possible. Deviations from Hardy–Weinberg equilibrium (HWE) at each microsatellite locus were calculated using Genepop v. 4.7 (Raymond and Rousset, 1995). Micro-Checker v. 2.2.1 (Oosterhout et al., 2010) was used to evaluate possible genotyping errors, allele dropout, or null alleles. Fisher's exact test for linkage disequilibrium between all pairs of microsatellite loci was performed using Genepop v. 4.7. Associated probability values were corrected by sequential Bonferroni correction with a significance level of 0.05. A set of statistical tests on microsatellite loci, including number of alleles (N_a), polymorphism information content (PIC), and observed and expected heterozygosity (H_o and H_e) were performed in FSTAT v. 2.9.4 (Goudet, 1995). To measure the extent of differentiation and distribution of genetic variance within and between species, Arlequin v. 3.5 (Excoffier and Lischer, 2010) was used to perform analysis of molecular variance (AMOVA) and calculate F-statistics based on ND2 gene and microsatellite data. One thousand nonparametric permutations were used to determine statistical significance.

To identify hybrid individuals and assess the admixture of genotypes caused by interspecific hybridization, we carried out three cluster analyses of microsatellite data (351 samples): Principal Coordinate Analysis (PCoA) using GenAlEx v. 6.5 (Peakall and Smouse, 2012) and two different Bayesian methods, as performed in STRUCTURE v. 2.3.4 (Pritchard et al., 2000) and NewHybrids v. 1.1 beta (Anderson and Thompson, 2002). PCoA was conducted to visualize similarities or dissimilarities of individual data based on a distance matrix without prior information and identify individuals separated from pure populations, thus indicating admixed individuals. We used STRUCTURE v. 2.3.4 to assign individuals to K (unknown) populations in the absence of prior information on individual classification or origin. The admixture model with correlated allele frequencies was run in this analysis, and the number of a priori clusters (K) was assumed to be 2 (both species contribute to the individual gene pool), which was verified using the ΔK statistical method (Evanno et al., 2005). Each run was performed for 10 replicates of 1, 000, 000 MCMC iterations with a burn-in period of 100, 000. We evaluated the mean membership coefficient (Q) of each sampled population for the inferred clusters and assigned each genotype to an inferred cluster according to the threshold values of the individual proportion of membership (Q_i). We adopted a threshold to identify true-breeding individuals (Q_i ​≤ ​0.10 or Q_i ​≥ ​0.90) and mixed individuals (0.10 ​ < ​Q_i ​ < ​0.90) (Vähä and Primmer, 2006; Hoban et al., 2009).

NewHybrids v. 1.1 beta assumes two clusters with recent admixture and calculates the PP of categorical membership for each individual being assigned to six possible categories (Pure species A, Pure species B, F1, F2, Backcross to A, and Backcross to B). We calculated the posterior distribution using the default settings of the program, also without any prior information on individual or allele frequencies. To identify potential hybrids, we used the same threshold as for STRUCTURE analysis to confidently designate true-bred individuals. For hybrids, an individual with 0.50 ​≤ ​PP ​ < ​0.90 was classified into the corresponding cluster; otherwise, the individual was not assigned to any identified hybrid category. We ran a 100, 000-iteration burn-in period, followed by 900, 000 MCMC repetitions.

3.   Results

3.1   Genetic diversity and species differentiation

The level of genetic diversity observed in E. jankowskii was unexpectedly higher than that in E. cioides (Table 2). The H_d and π of E. jankowskii based on ND2 were 0.910 and 0.0057, respectively, which were higher than those of E. cioides (H_d ​= ​0.865, π ​= ​0.0019; Table 2). Microsatellite genotyping results illustrated that all 15 loci were successfully amplified and showed polymorphism (Table 3). Each locus was inherited independently, and no linkage disequilibrium was found between all pairs of loci. Significant deviations from HWE were confirmed at loci Mme8, Escμ03, and LS2 in E. jankowskii after Bonferroni correction. We failed to detect a large allele dropout or scoring errors caused by stuttering; however, Mme8 and LS2 in E. jankowskii and Sosp13, Escμ03, and Escμ04 in E. cioides all showed the possible presence of null alleles. Analysis indicated heterozygote deficiency at these loci. All 15 loci were considered suitable for subsequent analysis. The N_a ranged from 4 to 28 in E. jankowskii and from 3 to 19 in E. cioides. Observed heterozygosity (H_o) in E. jankowskii ranged from 0.365 to 0.898 (mean ​= ​0.696) and expected heterozygosity (H_e) from 0.443 to 0.956 (mean ​= ​0.769), while H_o in E. cioides varied from 0.140 to 0.930 (mean ​= ​0.609) and H_e from 0.133 to 0.934 (mean ​= ​0.652). In comparison, we mainly consider π and H_e, which are not significantly correlated with sample size. E. jankowskii showed higher genetic diversity than E. cioides regardless of whether the whole or subdivided population was considered (Table 2). Among the five subdivided populations of E. jankowskii, GHT and KD populations had higher nucleotide diversity, which were 0.0096 and 0.0119 respectively.

Table  2.  Summary of ND2 and microsatellite diversity in local populations.

Species		Population		ND2						Microsatellite
				N	S	H	Hd	π		N	Na	PIC	Ho	He
E.jankowskii
		ALH		37	14	10	0.874	0.0030		23	7.7	0.653	0.697	0.705
		CG		71	38	33	0.885	0.0032		67	12.0	0.738	0.684	0.773
		GHT		63	75	13	0.886	0.0096		42	8.0	0.707	0.742	0.761
		KD		3	18	3	1.000	0.0119		3	3.4	0.515	0.689	0.680
		LB		2	4	2	1.000	0.0039		2	4.3	0.614	0.707	0.730
		Total		176	97	38	0.910	0.0057		137	13.4	0.737	0.696	0.769
E.cioides
	Total (GHT)	36	23	20	0.865	0.0019	43	8.7	0.604	0.609	0.652
N: Number of samples; S: Number of polymorphic sites; H: Number of haplotypes; Hd: Haplotype diversity; π: Nucleotide diversity; Na: Number of alleles; PIC: Polymorphism information content; Ho: Observed heterozygosity; He: Expected heterozygosity.

 | Show Table

DownLoad: CSV

Table  3.  Genetic diversity of Emberiza jankowskii and Emberiza cioides based on 15 microsatellite loci.

Locus	E.jankowskii							E.cioides
	Na	PIC	Ho	He	HWE	F(null)		Na	PIC	Ho	He	HWE	F(null)
Sosp01	10	0.738	0.672	0.777	NS	0.074		10	0.732	0.651	0.769	NS	0.081
Sosp04	9	0.601	0.693	0.658	NS	−0.032		4	0.439	0.512	0.549	NS	0.032
Sosp14	14	0.860	0.818	0.876	NS	0.032		11	0.644	0.651	0.682	NS	0.016
Sosp13	12	0.731	0.715	0.763	NS	0.033		4	0.449	0.349	0.514	NS	0.200
Mme2	27	0.950	0.898	0.956	NS	0.029		15	0.834	0.930	0.859	NS	−0.050
Mme12	5	0.402	0.474	0.443	NS	−0.041		3	0.313	0.465	0.379	NS	−0.097
Cuμ4	28	0.913	0.898	0.922	NS	0.011		11	0.821	0.837	0.849	NS	−0.003
Mme8	4	0.533	0.365	0.589	***	0.236		5	0.468	0.628	0.568	NS	−0.077
Sosp02	17	0.896	0.869	0.907	NS	0.021		6	0.663	0.744	0.715	NS	−0.040
Dpμ16	9	0.732	0.832	0.765	NS	−0.047		8	0.612	0.605	0.658	NS	0.049
Escμ01	6	0.610	0.635	0.669	NS	0.018		5	0.507	0.512	0.563	NS	0.057
Escμ03	8	0.749	0.686	0.786	*	0.066		7	0.633	0.558	0.692	NS	0.113
Escμ04	20	0.874	0.832	0.888	NS	0.030		19	0.918	0.674	0.934	NS	0.158
Ecit2	9	0.601	0.577	0.650	NS	0.068		3	0.125	0.140	0.133	NS	−0.027
LS2	23	0.868	0.482	0.880	**	0.286		19	0.897	0.884	0.915	NS	0.012
Na: Number of alleles; PIC: Polymorphism information content; Ho: Observed heterozygosity; He: Expected heterozygosity; HWE: Hardy-Weinberg equilibrium; F(null): Null allele frequency; NS: Not significance; *: Statistical significance, p < 0.05; **: Statistical significance, p < 0.01; ***: Statistical significance, p < 0.001.

 | Show Table

DownLoad: CSV

The AMOVA based on ND2 haplotype frequencies revealed that intraspecific variation accounted for 28.69% of the total genetic variation, while interspecific variation accounted for 71.31%. The fixation index (F_st) value was 0.713, which was statistically significant after Bonferroni correction (p ​ < ​0.05). Similarly, for microsatellite analysis, according to AMOVA, 7.48% of the variation was attributed to "among individuals within species", 77.00% to "within individuals", and 15.52% to "among species". The global F_st value was 0.155 (p ​ < ​0.05). In general, genetic differentiation among populations is considered low (0 ​≤ ​F_st < 0.05), medium (0.05 ​≤ ​F_st < 0.15), high (0.15 ​≤ ​F_st ≤ 0.25), or very high (F_st > 0.25) (Rousset, 1997). The aforementioned results, as shown in Table 4, revealed high genetic differentiation between the two species.

Table  4.  Analysis of molecular variance (AMOVA) and the F_st for Emberiza jankowskii and Emberiza cioides.

Maker	Source of variation	df	Sum of squares	Variance component	Percentage of variation	Fst
ND2 gene	Among species	1	1992.28	33.11	71.31
	Within species	210	2796.64	13.32	28.69
	Total	211	4788.92	46.43		0.713*
Microsatelite	Among species	1	139.69	1.02	15.52
	Among individuals within species	178	1077.11	0.49	7.48
	Within individuals	180	912.00	5.07	77.00
	Total	359	2128.80	6.58		0.155*
df: Degrees of freedom; *: Statistical significance, p < 0.05.

 | Show Table

DownLoad: CSV

3.2   Haplotype relationships

We identified 53 haplotypes from 212 ND2 sequences. All haplotypes were roughly divided into two distinct clusters; haplotypes 1–34, 36, and 38 were uniquely clustered together for E. jankowskii, while haplotypes 39–53 belonged uniquely to E. cioides (Fig. 2). The respective haplotypes in both species showed star-like clusters. Although there were at least 54 mutations between the two haplogroups, we still found two shared haplotypes (Haplotypes 35 and 37), indicating mitochondrial introgression from E. cioides to E. jankowskii. Haplotype 35 was present in four E. jankowskii individuals and 15 E. cioides individuals, while Haplotype 37 was present in one E. jankowskii individual and four E. cioides individuals.

Figure  2.  Median joining network based on 53 haplotypes, inferred from ND2 gene sequences. A bar on line between haplotypes represents one mutation. Area of circle is proportional to the number of samples with the given haplotype and the smallest size corresponds to a single sample.

DownLoad: Full-Size Img PowerPoint

3.3   Hybrid detection and admixture analyses

The PCoA results for the microsatellite data are shown in Fig. 3. Coordinate 1 explained 10.94% of the total variation, while coordinate 2 explained 5.30%. Except for one uppermost outlier (corresponding to individual 387), most individuals were divided into two clusters corresponding to the two previously assigned morphological classifications (E. jankowskii and E. cioides, respectively), indicating that interspecific differentiation was stronger than intraspecific differentiation. Although we can only observe one individual intermediate between the two clusters, we cannot rule out the possibility of later-generation hybrids.

Figure  3.  Scatterplot of PCoA performed on 15 microsatellite loci for Emberiza jankowskii and Emberiza cioides. Individuals are represented by yellow triangles Emberiza jankowskii and blue squares (Emberiza cioides) respectively. Coord.1 and Coord.2 account for 10.94 and 5.3% of the total variation, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

According to the peak value of ΔK, we determined that K ​= ​2 was the most likely number of clusters, representing the two respective species (Appendix Fig. S1). Most individuals of E. jankowskii and E. cioides were classified into their respective pure populations with high probability (Q_i ​≤ ​0.10 or Q_i ​≥ ​0.90; Fig. 4). Only a few individuals had a probability (0.10 ​ < ​Q_i ​ < ​0.90) of being assigned to the admixed cluster. Among them, individual 36, assigned with Q_i of 0.58 and 0.42 to E. jankowskii and E. cioides, respectively, was considered to have multiple genetic admixtures, and the remaining individuals (individuals 4, 7, 13, 20, 31, 281, and 387) were considered admixed individuals with 0.75 ​ < ​Q_i ​ < ​0.90.

Figure  4.  Comparison of admixture analysis based on ND2 gene and 15 microsatellites. For each sample, its identity according to ND2 genes is shown above the STRUCTURE and NewHybrids diagrams. Bar plots represent STRUCTURE inferences of individual assignments (K ​= ​2). Each column (representing the individual) is divided into K segments of different colors, whose length is proportional to its membership in the corresponding genetic cluster. In NewHybrids diagram, each individual corresponds to a vertical bar, which is divided into segments of different colors and the length of each segment is proportional to the posterior probability of belonging to each category (see color legend). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

The NewHybrids analysis classified all individuals into six categories (Fig. 4). The results showed that most individuals morphologically identified as E. jankowskii or E. cioides were indeed classified as corresponding true-bred individuals with PP ​≥ ​0.90. Seven of 351 individuals (excluding missing data) could not be identified as pure parental species according to the previously defined threshold. Individuals 36 and 387 were assigned by NewHybrids to the F2 generation with PP ​ > ​0.85. However, F1-hybrids appear to be quite rare and were not detected in our analyses. Because the probability of two F1-hybrids interbreeding is thus extremely low, we have reclassified these F2-hybrids as backcrosses. The PPs of individuals 384 and 386 being classified as hybrids were 0.38 and 0.46, respectively, and the remaining three individuals (individuals 4, 13, and 20) ranged from 0.15 to 0.20, but there was no evident definition of the six genotype classes. No individual was assigned to the F1 hybrid class. In terms of phenotype, four of the seven possible hybrids are E. jankowskii and the other three are E. cioides. Combined with the admixture of mitochondrial ND2 gene, only one E. cioides individual among the five individuals in Haplotype 37 (Fig. 2) was classified as a hybrid in the NewHybrids analysis (corresponding to individual 384), and the remaining two E. cioides individuals in Haplotype 37, and 15 E. cioides individuals in Haplotype 35 were classified as true-bred. However, four E. jankowskii individuals of the two shared haplotypes (three in Haplotype 35 and one in Haplotype 37) lacked corresponding microsatellite data (red bar, Fig. 4). The percentage of admixed individuals is indicated in Fig. 4. Both STRUCTURE and NewHybrids results suggested weak gene introgression between the two species, and seven individuals (4, 13, 20, 36, 384, 386, and 387) were considered possible hybrids by NewHybrids analysis.

Almost all possible hybrids came from ALH and GHT populations, and only one hybrid was found in LB. This may because we have fewer samples at LB (Table 2).

4.   Discussion

4.1   High genetic diversity

Genetic diversity is related to a population's environmental adaptability and evolutionary potential, which further affects the long-term survival of the species (Frankel and Soule, 1981). In general, endangered species have depleted genetic resources and often exhibit low genetic diversity (Keller and Waller, 2002). However, this study reconfirmed the high genetic diversity of endangered species E. jankowskii, which is higher than that of the more common and sympatric E. cioides (Li et al., 2021). Microsatellite analysis has revealed that populations of endangered species often do not conform to HWE because of subpopulation structure, inbreeding, and the presence of null alleles (Lade et al., 1996; Paxton et al., 1996; Ardren et al., 1999). Given the small population size and narrow distribution of E. jankowskii, we considered that deviations from HWE and high null allele frequencies could result from inbreeding and recent demographic fluctuations. We retained problematic loci (Sosp13, Mme8, Escμ03, Escμ04, and LS2) in our analyses because these deviations are unlikely to significantly affect subsequent genetic analysis (Chapuis and Estoup, 2007; Murphy et al., 2019). According to current data, genetic diversity estimates based on 15 microsatellite loci showed that PIC, H_o, and H_e were higher in E. jankowskii than in sympatric E. cioides. Of course, we recognize that even if the sampling strategy of unrelated individuals is adopted, the genetic diversity of E. cioides may still be underestimated to some extent under this biased sampling mode (sample size difference). In addition, the effect of genetic drift at a species' range margins can not be ignored. According to the abundant center hypothesis, marginal populations tend to have lower diversity indices than central populations (Beatty et al., 2010; Langin et al., 2017; Rönkä et al., 2019). Our sampling site is located at the edge of the published range map of E. cioides. Therefore, the E. cioides population we sampled may have a poor gene pool. Nevertheless, E. jankowskii still exhibits higher genetic diversity than some other endangered species or even least concern species, such as Taeniopygia guttata (H_d ​= ​0.692, π ​= ​0.0060), Parus montanus (π ​= ​0.0021), and Neothraupis fasciata (π ​= ​0.0045, H_d ​= ​0.868) (Zink, 2005; Newhouse and Balakrishnan, 2015; Lima-Rezende et al., 2019). On the one hand, the species' large historical population size and relatively stable demographic history may have contributed to current genetic variation (Li et al., 2021), and historical population declines may not have been severe enough to cause appreciable effects on genetic variation. Previous studies have revealed the absence of population structure in this metapopulation, suggesting good connectivity among local populations (Li et al., 2021). Hence, gene flow within this metapopulation brings several subpopulations together to form a large gene pool that increases genetic stability (Lamy et al., 2012; Jangjoo et al., 2016; Wesselmann et al., 2018; Peng et al., 2019), buffering the impacts of population decline on genetic diversity. However, maintenance of genetic diversity relies on the balance of influx versus efflux of genetic variation; here, exogenous gene flow compensates for the loss of genetic variation in the E. jankowskii population. Although gene flow may have arisen from unknown populations, comprehensive field records and the special habitat preference of this species almost deny this possibility. Additionally, the high diversity may be due to gene introgression from close relatives.

4.2   Significant differentiation and weak introgression

Incompletely differentiated species formed in geographic isolation often undergo interspecific hybridization and produce fertile offspring when secondary contact occurs (Mallet, 2005). Two species with divergence times of 5 Myr (passerines) to 17 Myr (nonpasserines) still retained the potential to produce offspring (Price and Bouvier, 2002). Considering the phenotypic and genetic similarities between E. jankowskii and E. cioides (Päckert et al., 2015), it is possible that introgressive hybridization occurred in sympatric regions. The haplotype network based on ND2 gene data showed two shared haplotypes for both species, which could result from either introgressive hybridization or the preservation of an ancestral polymorphism. Nevertheless, additional evidence suggests that introgressive hybridization rather than ILS is the main reason for haplotype sharing between the two species. The fixed difference of 54 mutations in ND2 gene haplotypes, high F_st value, and clear population clustering between the two species indicate that they have experienced a long period of independent evolution. As recently evolved lineages diverge, the transition to monophyly requires only 8–12 N_e generations, where N_e is the effective population size (Chesser and Baker, 1996; Hoelzer, 1997). Dated multi-locus phylogenies have shown a divergence time of more than three million years between the two species (Päckert et al., 2015, 2020; Cai et al., 2021), a period sufficient for them to develop reciprocal monophyly (Pamilo and Nei, 1988). Subsequent analysis of microsatellite data using STRUCTURE and NewHybrids also indicated that there was little genetic admixture. These are clearly not signals of ILS, as the retention of ancestral polymorphisms would allow "traces" of ancient genotypes to be evenly dispersed across the genome. Interestingly, we seem to have found some evidence of introgressive hybridization from the genetic statistics of subdivided populations. Compared with other populations, the nucleotide diversity of GHT and KD populations in E. jankowskii is relatively high (Table 2). However, the difference is that local admixture may contribute to the rich nucleotide variation of GHT population, while that of KD population may come from the statistical deviation of small samples, because hybrid individuals are not found in KD population. In conclusion, our evidence suggests that weak introgressive hybridization occurred between the two species.

An alternative approach to distinguish introgressive hybridization from ILS is based on hybrid zone studies. Multiple alleles inherited from a common ancestor will be randomly distributed among descendant populations, with no specific biogeographic pattern (Barbujani et al., 1994; Toews and Brelsford, 2012). The null hypothesis of ILS can be rejected if ancestral alleles are concentrated only within a known hybridization zone, by arguing that TSPs are derived from introgressive hybridization. However, because we did not perform allopatric sampling, we cannot yet clarify whether there is allele sharing among populations from other regions. Furthermore, local gene pools across the large range of E. cioides might differ substantially. Admixture in E. cioides populations might be underestimated because we obtained samples from only one known overlapping breeding ground (GHT), and "pure" allopatric E. cioides populations were not included in our data (Ravinet et al., 2018; Päckert et al., 2019). We therefore do not know whether the latter populations are free from any gene flow and could affect the signal in a data set. Future allopatric sampling is therefore warranted.

The combined use of mitochondrial and nuclear markers has spawned a growing number of studies showing mitonuclear discordance patterns (Toews and Brelsford, 2012). In the present study, we seemingly observed mitonuclear discordance in the direction of gene introgression. In the haplotype network, five E. jankowskii individuals carrying the ND2 haplotypes of E. cioides clustered in the E. cioides branch, but no E. cioides individuals clustered in the E. jankowskii branch, indicating that E. cioides mtDNA introgressed into the E. jankowskii population unidirectionally; however, in microsatellite analysis, both species showed signals of gene introgression from each other. This introgression pattern may be mediated by female E. cioides individuals (Dabrowski et al., 2005; Gainsford et al., 2020). However, this inference of mitonuclear discordance in patterns of introgression is speculative, as we could not clearly compare the intensity of mitochondrial and nuclear gene introgression. Mostly, our data indicate that the two species are highly distinct and have limited gene flow between them. Bias in local abundance and reproductive pressure among numerically imbalanced populations could enhance hybridization among the two bunting species (Lepais et al., 2009; McCracken and Wilson, 2011), and the field findings showed that E. cioides was rarer in these peripheral populations as compared with E. jankowskii in its core distribution range. This would be an explanation for the asymmetrical introgression from E. cioides to E. jankowskii (Peters et al., 2017). In addition, according to Haldane's Rule, heterogametic sex (females in birds) tends to have lower hybrid fitness than the homogametic sex (males in birds), hence gene flow through females could be relatively low, which may be one of the reasons for the weak hybridization signal in this study (Delph and Demuth, 2016).

4.3   Conservation implications

The important role played by introgressive hybridization in genetic rescue of threatened species has attracted considerable attention (Grant and Grant, 2014, 2019; Harrisson et al., 2016; Pierce et al., 2017). Small populations are more susceptible to exogenous gene flow due to the shallow gene pool. Although the introgressive hybridization between E. jankowskii and E. cioides is so weak that we cannot clearly determine the genetic contribution of introgression, we suggest that introgression may nonetheless contribute to the maintenance of high genetic diversity in E. jankowskii (Tompkins et al., 2006; Yoshida et al., 2016; Pohjoismäki et al., 2021), and that natural rather than artificial hybridization can be used cautiously as an option for short-term genetic rescue of endangered species under extremely urgent conditions (Benson et al., 2011; Pierce et al., 2017; Chan et al., 2019; Grant and Grant, 2019). Nevertheless, the genetic consequences of introgressive hybridization under different conditions are controversial. In some cases, it may be associated with low reproductive success and even an increased risk of extinction (Quilodrán et al., 2018; Arantes et al., 2020). Although there is no evidence of hybridization threatening the genetic integrity of endangered E. jankowskii, current population declines and potential negative effects of hybridization cannot be ignored. Considering its small population size and narrow distribution, we believe that expanding the habitat size of E. jankowskii remains the most promising conservation strategy at present.

5.   Conclusions

Our study confirms that weak introgressive hybridization occurs between sympatric E. jankowskii and E. cioides, which may be mediated by female E. cioides individuals. Although the genetic contribution of introgressive hybridization on the maintenance of genetic diversity of endangered E. jankowskii may not be as significant as expected, we believe that even relatively weak gene flow is still an important source of genetic variation in small populations. Further research will focus on revealing introgressive hybridization with more sensitive molecular markers such as Single Nucleotide Polymorphisms and analysis including clear historical population dynamics.

Author contributions

LH conceived and designed the experiments, carried out the molecular genetic analysis, participated in the manuscript's design, and drafted the manuscript. SL and YJ developed the ideas and obtained funding support. LZ, and FY helped collect samples and classify species. WS and DL helped perform the microsatellite experiments. All authors read and approved the final manuscript.

Declaration of competing interest

The authors declare that they have no competing interests.

Ethics statement

This study conformed to the guidelines for the care and use of experimental animals established by the Ministry of Science and Technology of the People's Republic of China (Approval number: 2006–398). The research protocol was reviewed and approved by the Ethics Committee of Jilin Agriculture University. The reporting in the manuscript follows the recommendations in the revised Animals (Scientific Procedures) Act 1986 in the UK and Directive 2010/63/EU in Europe. Field studies comply with the IUCN Policy Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora.

Acknowledgements

We thank Xue Zhang and Yufei Cui for their support in the field. We thank Gabe Yedid, Ph. D, from Liwen Bianji (Edanz) (www.liwenbianji.cn) for editing the English text of a draft of this manuscript. This work was supported by the National Natural Science Foundation of China (Grant Nos. 31601856 and 31670398).

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.avrs.2022.100035.

Glossary