Genetic variation of endangered Jankowski’s Bunting (Emberiza jankowskii): High connectivity and a moderate history of demographic decline

Introduction

Genetic diversity plays an important role in the environmental adaptation and the long-term evolution of species (Mable, 2019). Climate change and human activities since the Industrial Revolution are major contributors to the loss of genetic variation in wild populations worldwide (Leigh et al., 2019). Under the influence of habitat loss and fragmentation, endangered species with small population sizes and narrow distributions are more susceptible to genetic drift and inbreeding depression than common species, leading to the loss of intraspecific variation (Frankham et al., 2002). Furthermore, the isolation effect between habitat patches often blocks gene flow between populations, ultimately resulting in genetic isolation and accelerated extinction (Nei, 1973). Numerous empirical studies have shown that low genetic diversity and gene flow between populations appear to have become the “tags” of endangered species (Wang et al., 2015; Hague and Routman, 2016; Basiita et al., 2017). However, several recent studies in birds, amphibians, and mammals have shown that large historical effective population size (N_e) and long generation time can act as buffers against the loss of genetic variation, leading to an observed high genetic diversity in threatened populations (Brekke et al., 2011; Torres-Florez et al., 2014; Zhang et al., 2015). In addition, mild population dynamics (a little variation in N_e over time), ecological selection for increased heterozygotes, and exogenous gene flow may also explain the genetic stability of declined populations (Schroeder et al., 2010; Peters et al., 2014). Demographic processes, along with species characteristics such as mutation rates and life history, together determine the current amount and distribution of genetic variation (Hague and Routman, 2016). Tracing back demographic bottlenecks and expansions, as well as more complex fusion and differentiation events, can deepen our understanding of the current patterns of genetic variation.

Owing to excessive and intensive land-use (Yin et al., 2018), mainly involving agricultural production, overgrazing, and other human activities (Wang et al., 2011), wild species in Inner Mongolia are threatened by habitat loss, fragmentation, and degradation (Jiang et al., 2008). Jankowski’s Bunting, Emberiza jankowskii, is a small resident bird with historical breeding grounds only in northeastern China, the Far-East Russia, and north-eastern North Korea. Its population has been in decline since the 1960s, and its distribution has shrunk dramatically (Jiang et al., 2008; Wang et al., 2011; Han et al., 2017). According to a recent and extensive population survey report, the distribution of Jankowski’s Bunting population in China is limited to 16 towns near Zhaluteqi in eastern Inner Mongolia (Han et al., 2017). Nevertheless, field observations and stable isotopes (δ²H) showed that Jankowski’s Bunting often spreads in short distance or even far away from the breeding area, rather than all being sedentary (Han et al., 2017, 2021). Its wintering habitat can reach as far as Beijing (Chao et al., 2016). The Jankowski’s Bunting population exhibited an island-like distribution due to habitat fragmentation. Based on remote sensing and bird density, the population size was estimated to be approximately 9,800–12,500 individuals, while the suitable habitat area was only 280 km² (Han et al., 2017). Jankowski’s Bunting is currently listed as ‘endangered’ on the IUCN Red List of Threatened Species. The species prefers low altitude hills and grasslands with grass (such as Stipa baicalensis and Spodiopogon sibiricus) and small shrubs (such as Prunus sibirica; Gao et al., 2003). The main threats to its survival include habitat fragmentation and degradation, nest parasitism, grazing disturbance, drought, wildfire, and mercury exposure (Jiang et al., 2008; Feng et al., 2020; Zhang et al., 2020; Han et al., 2022). Currently, although the authorities have taken measures to protect grassland habitats, which will help the population recover in the short term, the genetic impact of ecological pressure on the population and the long-term population dynamics are still unclear. Previous limited studies have revealed an unexpectedly high genetic diversity and absence of genetic structure in the metapopulation consisted of five local populations around Zhaluteqi in eastern Inner Mongolia (Li et al., 2021). However, snapshots of the current genetic variation in Jankowski’s Bunting are insufficient to provide robust evidence of genetic health. Evaluation of population dynamics would contribute to the comprehensive assessment of population health and targeted implementation of conservation efforts.

Technological advances in genetic analysis enable us to infer the demographic history of a species from a contemporary single-sample (Mondin et al., 2018). Bayesian skyline plots (BSPs) inferred from mtDNA D-loop sequences have revealed a long-term stable demographic history of Jankowski’s Bunting (Li et al., 2021). This stability appears to derive from the ice-free geological history of Northeast Asia during the Pleistocene (Wang et al., 2013; Dong et al., 2017). However, BSP analysis has many limitations, such as an over-reliance on molecular clock calibration and sampling strategies (Grant, 2015). When the population is decimated by glaciation, the decline is not severe enough to qualify as a bottleneck but sufficient to erase previous demographic history (Grant et al., 2012). Additionally, simulation studies showed that small sample size cannot provide the power needed to detect population expansion, and the loss of information during population fluctuations often leads to a flat BSP being overinterpreted as population stability (Grant, 2015). Although large-scale geological history plays an irreplaceable role in shaping the current patterns of genetic variation in populations, the impact of recent habitat changes may be more pronounced (Zarlenga et al., 2014; Feng et al., 2018; Leigh et al., 2019). So far, the lack of fine-scale demographic maps of Jankowski’s Bunting has led to a vague understanding of recent history. The DIYABC program based on approximate Bayesian computation (ABC) accepts user-defined complex evolutionary scenarios with loose requirements for samples, and provides estimates of posterior probability and parameter confidence, which can be used to reconstruct robust evolutionary history using coalescent simulation (Cornuet et al., 2008). Routine detections of expansion and bottleneck events, such as neutrality test and nucleotide mismatches, which have strong statistical power for low frequency mutations in non-recombinant regions, could provide supplementary information for reconstruction of demographic history (Ramos-Onsins and Rozas, 2002; Ramírez-Soriano et al., 2008).

Assessing and increasing the connectivity between isolated populations is an important measure to help small populations resist to habitat fragmentation (Beier et al., 2008; Newmark, 2008). The commonly used estimators $4N_m=\left(1-F_{\mathrm{st}}\right)/F_{\mathrm{st}}$ and coalescent-based Migrate-n yielded extremely high estimates of historical migration rates over the past 4N_e generations; this time period may be over hundreds of years in the past for the Jankowski’s Bunting, which is much longer than the duration of recent human disturbance (Li et al., 2021). However, the impact of recent (<200 years) environmental changes on contemporary population connectivity is unknown. In addition, both estimators are based on constant migration rates and are subject to numerous assumptions, which may be unreasonable for populations in highly disturbed habitats (Wilson and Rannala, 2003). In contrast to the estimators of long-term gene flow, disequilibrium-based Bayesian estimators allow for arbitrary allele frequency distributions and deviations from Hardy–Weinberg equilibrium (HWE; Wilson and Rannala, 2003), which may be more appropriate for the populations we study. The comparison of historical and contemporary gene flow provides an opportunity to understand the impacts of recent habitat changes caused by human activity on population connectivity (Chiucchi and Gibbs, 2010; Epps et al., 2013a,b). Due to the high mobility of birds, as well as the narrow distribution of local populations and the lack of geographic barriers, we expect that recent habitat changes may not have a significant impact on population connectivity.

Evolutionary forces act differently on different molecular markers (Hoarau et al., 2004), thus a single marker cannot fully reflect the entire evolutionary process of a population. The combination of molecular markers with different mutation rates can provide a comprehensive insight into the evolution history of the studied population (Hindrikson et al., 2017). In this study, we used a mitochondrial control region (partial mtDNA CR) and 15 nuclear microsatellite markers to explore the recent demographic history of Jankowski’s Bunting, and we compared the historical and contemporary gene flows between populations to reveal the impact of habitat change on population connectivity. Specifically, we aimed to test the following hypotheses: (1) Jankowski’s Bunting has a large historical N_e and a moderate demographic history; and (2) recent habitat change might have no significant impact on the species’ population connectivity. Despite the lack of historical census data, the high genetic diversity, absence of genetic structure, historical frequent records, and the availability of grassland-dominated habitats may support the rationality of the above assumptions (Li H. et al., 2019; Li et al., 2021). This study aims to provide a clearer picture of the demographic history and genetic imprints of recent environmental changes in the species, as well as lead to effective conservation efforts.

Materials and methods

Sample collection and DNA extraction

Blood samples from 117 Jankowski’s Buntings were used in this study [including 59 samples described in Li et al., 2021], which were collected from 2012 to 2019 at five localities in Inner Mongolia: Aolunhua (ALH), Gahaitu (GHT), Lubei (LB), Kundu (KD), and Chagantaoligao (CG; Figure 1). These sampling sites represented the known breeding grounds of Jankowski’s Bunting. During the breeding season each year, between May and July, nests were located and recorded by visual search and bird calls at the above locations. To minimize the interference of sampling on reproduction, adults were captured by setting mist nets around the nest after the breeding period. The blood samples (approximately 10 μl) were collected by wing vein puncture, then preserved in 1 ml anhydrous ethanol in a 1.5-ml sterilized centrifuge tube and stored in liquid nitrogen in the field. Total genomic DNA was extracted using AxyPrep Blood Genomic DNA Miniprep Kit (Aisijin, Hangzhou, China) following the manufacturer’s instructions.

MtDNA amplification, sequencing, and analysis

A few samples did not amplify due to the limited volume of blood. Of the 117 samples, 116 were used for CR amplification and sequencing. We used the novel primers A629 (5′-GTATTCAACCGACCACTT-3′) and S19 (5′-CTCCGAGACCTACGACCTG-3′) to amplify 540 base pairs (bp) of the CR based on sequences obtained previously (GenBank accession no. MH670589 to MH670647). DNA amplification of CR was performed in 20-μl volumes containing 10 μl of 2× Taq PCR Master Mix (Sangon Biotech, Shanghai, China), 0.25 μl of each primer (10 μM), 8.5 μl H₂O, and 1 μl template DNA (approx. 50 ng/μl). The PCR cycling program was as follows: an initial denaturation step at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 53.7°C for 30 s, and extension at 72°C for 30 s; then a final extension at 72°C for 10 min. The PCR products were electrophoresed on 1% agarose gels and sequenced by Sangon in Shanghai, China.

Sequences were checked and aligned using Bioedit Sequence Alignment Editor (Hall, 1999) and Clustal X v1.83 (Jeanmougin et al., 1998). It is important to note that most analyses were performed only on populations with sufficient sample size (i.e., ALH, CG, GHT). KD and LB populations were not used for neutrality tests, mismatch analyses, and estimates of N_e due to possible bias due to small sample sizes (<10). Nevertheless, we still retain KD and LB populations in other analyses (i.e., genetic diversity, migration, and genetic structure), as these two populations are essential for understanding population connectivity and stratification in the study area. To evaluate genetic diversity, common statistical parameters, such as the number of polymorphic sites (S), number of haplotypes (H), average number of nucleotide differences (k), haplotype diversity (H_d), and nucleotide diversity (π), were calculated with DnaSP v6 (Rozas et al., 2017).

To clarify the genetic relationships among the five local populations, a statistical parsimony (i.e., TCS) network of mtDNA CR was constructed in PopART (Leigh and Bryant, 2015) using a 95% parsimony connection limit. The TCS method provides an empirical assessment of parsimony biases and is considered the most informative method for exploring the phylogeographic history of species with low divergence (Posada and Crandall, 2001).

The distribution of variation among and within local populations was assessed by Analysis of Molecular Variance and F-statistics in Arlequin v3.5 (Excoffier and Lischer, 2010) with 1,000 permutations to assess statistical significance of fixation indices.

Neutrality tests based on allele frequency distribution were carried out for both the local populations and the metapopulation with Arlequin to calculate the Tajima’s D and Fu’s F_s values and their significance. Tajima’s D and Fu’s F_s, as methods of testing neutral evolution, can be used to trace the demographic events (Fu, 1997; Liu et al., 2006). Tajima’s D test is more powerful for assessing the long-term history, while Fu’s F_s is sensitive to recent fluctuations (Depaulis et al., 2003). Arlequin was also used to analyze the mismatch distribution of pairwise nucleotide differences under the constant population size model (Rogers and Harpending, 1992). Subsequently, a mismatch distribution diagram was plotted in DnaSP and used to visualize the historical population dynamics. If an expansion event was detected, the formula $\tau =2\mu kt$ (Rogers and Harpending, 1992) was used to estimate the population expansion time (t), where τ (tau) is the expansion constant under the sudden expansion model calculated by Arlequin, μ is the mean mutation rate of nucleotides per generation, and k is the sequence length in nucleotides. We did not use spatial expansion model due to the lack of reliable information on population history. Due to the lack of available mutation rates, the average mutation rate of 8.72 × 10⁻⁸ substitutions per site used for the first domain of the mtDNA CR of the yellow warbler (Dendroica petechia, Passeriformes) and 1.156 × 10⁻⁸ substitutions per site calibrated for the coal tits (Periparus ater, Passeriformes) were applied to the above calculations to generate robust estimates (Milot et al., 2000; Cabanne et al., 2007; Pentzold et al., 2013; Tritsch et al., 2018).

Microsatellite amplification, genotyping, and analysis

For microsatellite analysis, a total of 15 polymorphic loci and 103 of the 117 samples were used in this study without multiplex amplification. We selected the primers used in the previous study for this species and its relatives. These primers represent most of the available loci of this species and are suitable for genetic polymorphism analysis through experimental verification. The specific information for each pair of microsatellite primers can be found in Supplementary Table 1. The PCR procedure and reaction volume previously described for CR amplification, as well as the annealing temperature optimized for each primer pair, was used to perform subsequent microsatellite amplification. The genotyping procedure was executed on ABI 3730xl DNA Analyzer platform. Final fragment sizes were calculated in GeneMarker v1.91 (SoftGenetics LLC, State College, PA, United States) and scored manually by Sangon Biotech (Shanghai, China).

We checked for deviations from the HWE resulting from heterozygote excess or deficiency at each microsatellite locus and overall with Genepop v4.7 (Raymond and Rousset, 1995). Fisher’s exact tests of linkage disequilibrium (LD) were also implemented in Genepop for all pairwise combinations of loci. The exact p-values were corrected by sequential Bonferroni tests. Micro-checker v2.2.1 (Oosterhout et al., 2010) was used to investigate the presence of possible null alleles, large allele dropout, or stuttering for each locus. Estimation of locus and population-based diversity indices including the number of alleles (N_a), allelic richness (A_r), observed heterozygosity (H_o), and expected heterozygosity (H_e) was performed in FSTAT v2.9.4 (Goudet, 1995). In addition, Cervus v3.0.7 (Kalinowski et al., 2007) was used to calculate the polymorphic information content (PIC), a parameter used to describe the informativeness of microsatellite loci.

Two complementary methods based on microsatellite data were used to detect potential genetic structure. First, Analysis of Molecular Variance and F-statistics were implemented in Arlequin, performing 1,000 and 3,000 permutations of significance tests, respectively, to evaluate the genetic similarity of local populations. Second, we used a Bayesian approach implemented in STRUCTURE v2.3.4 (Pritchard et al., 2000) to obtain reasonable group clustering and visualize the genetic connections between local populations. Genetic subdivision was evaluated by estimating the likelihood and sample composition of independent runs of subgroups (K = 1–10), assuming an admixture model with correlated allele frequencies. Each run was repeated 10 times, with 1,000,000 Markov chain Monte Carlo (MCMC) iterations and a burn-in period of 100,000 to check for convergence. The optimal K value was determined according to the ΔK statistical method (Evanno et al., 2005) as well as the likelihood values.

Despite the strong migratory capacity of birds, recent habitat fragmentation may still have cryptic effects on genetic communication between populations, resulting in a certain geographical distribution pattern of intraspecific variation. To test for a pattern of isolation by distance, a Mantel test of genetic (Codom genotypic genetic distance matrix) and geographic (Unit: km) distance between populations was carried out by using the software package GenAIEx v6.5 (Peakall and Smouse, 2006; Peakall and Smouse, 2012). 1,000 permutations were used to test for statistical significance.

Historical and recent gene flow tests were performed in two programs using a Bayesian approach and MCMC sampling implemented in Migrate-n v3.7.2 (Beerli, 2006) and BayesAss v3.0 (Wilson and Rannala, 2003), respectively. The two programs can produce comparable estimates of gene flow parameters (e.g., proportion of immigrants in a population per generation) at different time scales: Migrate-n produces an N_e and gene flow parameters among populations over approximately the past 4N_e generations, while BayesAss allows a recent (1–3 generations) estimation of migration rates (Samarasin et al., 2017). As recommended in the tutorials for Migrate-n, we ran three independent replicates with the Brownian motion mutation model. The search strategy was set to Bayesian inference. One long chain with 5,000 recorded steps, 500,000 sampled genealogies, and a sample increment of 100 were run to allow a wider search of the genealogy space. The number of discarded trees per chain (burn-in) was set to 5,000. Each MCMC run consisted of an analysis of adaptive heating with different temperatures (1; 1.5; 3; and 10,000) to assure convergence. BayesAss allows genotype frequencies to deviate from HWE or mutation-drift equilibrium (Wilson and Rannala, 2003). We performed three independent MCMC runs with 10,000,000 iterations and a sampling frequency of 100. The first 1,000,000 iterations were discarded as burn-in. According to the tutorial, we increased the proposal step size for the mixing parameters associated with proposed moves of the allele frequencies (0.7), inbreeding coefficients (1.0), and migration rates (0.3) to ensure an acceptance rate between 0.2 and 0.6. To ensure that the MCMC chains had converged, we checked the log probabilities for each run using Tracer v1.7.2 (Rambaut et al., 2018). The commonly used microsatellite mutation rate of 5 × 10⁻⁴ for vertebrates and the mutation-scaled formulae $M=m/\mu$ and $N_m=\theta \times M/4$ were used for the parameter transformation to produce comparable values for the migration rate (Chiucchi and Gibbs, 2010; Li et al., 2017; Speer et al., 2017; DeSilva and Dodd, 2020).

To obtain a global view of the demographic history of Jankowski’s Bunting, we estimated and compared the long-term and contemporary N_e for each local population as well as for the metapopulation. Since heterozygosity-based estimators are sensitive to the microsatellite mutation model, the long-term N_e in this study was calculated under two extreme mutation models, infinite allele model (IAM; Maruyama and Fuerst, 1985) and stepwise mutation model (SMM; Ohta and Kimura, 1973). Although the mutational pattern of microsatellites may vary from locus to locus, which is out of the scope of this study, the actual mutational process is generally intermediate between IAM and SMM (Bagley et al., 1999; Xu and Fu, 2004). The long-term N_e was estimated by the formula $N_e=\left({\left(\frac{1}{1-H_e}\right)}^2-1\right)/8\mu$ assuming SMM and the formula $N_e=\frac{H_e}{4\mu \left(1-H_e\right)}$ under IAM, respectively, where H_e is average expected heterozygosity and μ represents microsatellite mutation rate (Lehmann et al., 1998). Insufficient independent knowledge of mutation rate limits the accuracy of estimating N_e from multilocus data, so we used a range of microsatellite mutation rates for common taxa, including a typical value of 5 × 10⁻⁴ for vertebrates and a reliable range of 10⁻³–10⁻⁴ for tit species (Weber and Wong, 1993; Xu et al., 2000; Garza and Williamson, 2001; Hansson et al., 2014; Tritsch et al., 2018). The contemporary N_e was estimated in NeEstimator v2 (Do et al., 2014) under the random mating model. Two single-sample estimators based on linkage disequilibrium and molecular coancestry, respectively, and the standard temporal approach were used for the analyses (Waples, 2006; Nomura, 2008). Given the limited sample size and the fact that rare alleles can lead to an underestimation of LD-based values, we chose a minor allele frequency threshold of 0.01 to exclude single copy alleles. Compared with the single-sample estimator, the temporal estimator is less sensitive to migration and more accurate with long sampling interval (Luikart et al., 2010; Gilbert and Whitlock, 2015). The 1.7-year generation time of Jankowski’s Bunting was used for this estimation (Saether et al., 2005; Irwin et al., 2009). We assumed that there are three generations between the composite samples from 2012 to 2013 and the samples from 2019.

We used two methods based on allele frequency distribution and heterozygosity to detect bottleneck events. Under the IAM, SMM, and two-phase mutation model (TPM; Di Rienzo et al., 1994) in Bottleneck v1.2.02 (Piry et al., 1999), we tested for recent genetic bottlenecks (0.8–4.0 N_e generations), which did not require prior information on population size or genetic variation (Romane et al., 2010). TPM tests were run assuming 70% single-step mutations with 30% multi-step mutations. We performed the sign test and Wilcoxon sign-rank test to estimate heterozygosity excess or deficiency, the latter being most appropriate in the case of fewer than 20 microsatellite loci (Piry et al., 1999). Additionally, we checked for deviations in the allele frequency distribution from the normal L-shaped distribution (mode-shift). Rare alleles tend to be lost in a severely reduced population, resulting in fewer low-frequency alleles than intermediate-frequency alleles (Luikart et al., 1998). A mode-shift indicator, represented graphically, can detect the loss of rare alleles during a recent bottleneck, specifically over the last few dozen generations (Cornuet and Luikart, 1996; Luikart et al., 1998). All tests were based on 1,000 iterations.

To explore plausible scenarios of recent population dynamics in the Jankowski’s Bunting, ABC modeling was implemented in DIYABC v 2.1.0 (Cornuet et al., 2014) using a combined dataset of mtDNA and microsatellites. Due to the lack of prior information on the demographic history, we first defined seven possible scenarios (Scenario 1: constant N_e; Scenario 2: a recent bottleneck; Scenario 3: an old bottleneck; Scenario 4: a recent expansion; Scenario 5: an old expansion; Scenario 6: an old bottleneck followed by a recent expansion; Scenario 7: an old expansion followed by a recent bottleneck) for the metapopulation and then estimated the posterior probability of each scenario via two methods, direct estimation and logistic regression. Finally, we selected the best model for estimating historical demographic parameters. We adjusted the prior settings of some parameters before implementing the simulations. For genetic data and summary statistics, the motif size and allelic range for microsatellite loci were reset according to detailed locus information, and default settings were used for other parameters. All available parameters for microsatellites (e.g., mean number of alleles, mean genetic diversity, and mean size variance) and DNA sequences (e.g., number of haplotypes, mean of pairwise differences, and Tajima’s D) were used as summary statistics to generate the reference table. For the historical model of predefined scenarios, the Uniform distribution and adjusted intervals for time and N_e were applied to all parameters. The final parameter setting is shown in Supplementary Table 2. Kimura 2-parameter model was used as the mutation model of mtDNA CR. Following the program’s recommendations, 1,000,000 simulated datasets were applied to each scenario for robust results. Routine analyses, including pre-evaluation of scenarios and prior distributions, comparison of scenarios, estimation of parameters, and model checking, were performed to evaluate the confidence of the results. Parameters of interest were converted to calendar years based on an assumed 1.7-year generation time of Jankowski’s Bunting (Saether et al., 2005; Kondo et al., 2008; Irwin et al., 2009).

Results

Contemporary genetic diversity

The information of genetic diversity for all populations is shown in Table 1. In terms of the metapopulation, the mtDNA haplotype diversity (H_d) was 0.877, the nucleotide diversity (π) was 0.0046, and the average number of nucleotide differences (k) was 2.489. For the five local populations, the H_d, π, and k ranged from 0.833 to 0.871, 0.0028 to 0.0061, and 1.500 to 3.317, respectively. The GHT population had the highest π and k, 0.0061 and 3.317, respectively, which were higher than those of the CG population with the largest sample size.

For the microsatellite markers, all 15 loci were successfully amplified and bore high polymorphism (PIC > 0.5). We detected 191 alleles within 103 samples. Supplementary Table 3 shows the estimates for all parameters of interest, such as the observed heterozygosity (H_o) and expected heterozygosity (H_e). The number of alleles at each locus ranged from 4 to 27, with an average of 12.7. After sequential Bonferroni correction, loci Mme8 and LS2 significantly deviated from HWE (p < 0.01). Linkage disequilibrium tests of all pairs of loci revealed no significant p-values, suggesting that all loci were inherited independently. The null allele frequency exceeded 0.05 for four loci (Sosp13, Escμ04, Ecit2, and LS2). In isolated wild populations, loci affected by demographic bottlenecks, genetic drift, null allele, inbreeding, and natural selection are likely to deviate from HWE (Dąbrowski et al., 2015; Murphy et al., 2019). We retained these loci as these deviations were unlikely to significantly affect subsequent genetic analyses (Chapuis and Estoup, 2007; Murphy et al., 2019). The H_o and H_e of the five local populations ranged from 0.683 to 0.767 and 0.680 to 0.770, respectively, while those of the metapopulation reached 0.712 and 0.762 (Table 1). Although the small sample size of KD and LB populations may lead to biased estimations of diversity, in general, each local population of Jankowski’s Bunting showed similar and high genetic diversity.

Genetic similarity among populations

We defined 29 haplotypes based on the CR sequence data with 52 polymorphic sites. The statistical parsimony network showed a complicated topology, with most haplotypes radially clustered around three major haplotypes (Hap 1, 2, and 4; Figure 2). Except for Hap 27, which differed from the nearest potential haplotype by 21 nucleotide mutations, only 1–2 mutational steps occurred between most of the remaining haplotypes (Figure 2). All haplotypes were admixed without clear group clustering, suggesting genetic similarity among local populations.

The results of the genetic structure analyses, which were run to investigate the optimal population structure and eventual admixture patterns among the five local populations, are shown in Figure 3. The optimal K value was determined to be 3 by ΔK method (Figure 3B). No clear signals of subdivided populations were observed, either at the optimal K value (K = 3), or at the second (K = 2) or third (K = 6) most probable K values, indicating a high genetic homogeneity of the five local populations.

The test for population differentiation, as a clearer quantitative description, became a powerful complement to the STRUCTURE analysis. The Analysis of Molecular Variance and F-statistics results for the metapopulation are presented in Table 2. Regardless of the marker type, most of the genetic variation observed at the population level could be attributed to “within populations” rather than “among populations.” Variation among populations accounted for less than 5% of the total variation. The above results suggested that the five local populations present extremely limited genetic differentiation and should be considered as a whole.

The Mantel test revealed no significant correlation between genetic and geographic distance (r = −0.265, P (rxy-rand ≥ rxy-data) = 0.309) among the five local populations (Supplementary Figure 1), although the linear distance between the two most distant populations was more than 137 km.

Comparison of historical and contemporary gene flow

High gene flow is often associated with low differentiation among populations. Except for the low historical gene flow (N_m) of 0.178 from KD to GHT, we found relatively high values between the remaining population pairs (Table 3; Wright, 1978). According to the results of BayesAss, the proportion of individuals derived from other populations in each local population exceeded 0.25, implying considerable gene flow among populations (Table 3). Comparisons of historical and contemporary migration rates revealed a slight decrease in gene flow for most population pairs (Figure 4). Nevertheless, the distribution of the difference between the historical and contemporary migration rates was centered on zero (Supplementary Figure 2). Although immigration is secondary dominant over local retention in all habitat patches (local populations) as indicated by Table 3 and Figure 4, large-scale and asymmetric gene flow has occurred between local populations during the ancient and recent history. Migration rates among some populations may be underestimated due to the limited sample size, but overall, past landscape changes seem to not have resulted in severely limited population connectivity.

Demographic history

Historical N_e for the metapopulation, based on the common mutation rate range of 10⁻³–10⁻⁴ and the average H_e of 0.762, was estimated to range from 800 to 20,818 (Table 4), while contemporary N_e using linkage disequilibrium, molecular coancestry, and temporal method was estimated to be 144, 27, and 115, respectively. This suggested that contemporary N_e was greatly reduced compared with the historical values.

We reconstructed the demographic history of the Jankowski’s Bunting population by analyzing for recent population expansions and bottlenecks. For the ALH population, the neutrality test showed no signal of demographic expansion; however, significant negative values (p < 0.05) of Tajima’s D and Fu’s F_s were observed in the metapopulation, implying that it did not deviate from the null hypothesis, i.e., the metapopulation may have experienced a recent expansion (Table 5). Moreover, the unimodal pattern of mismatch distribution also supported the inference of a population expansion (Supplementary Figure 3). The bimodal pattern of mismatch distribution occurred only in the GHT population, while the unimodal pattern was found in all the other local populations as well as in the metapopulation. The expansion time of the metapopulation based on the two mutation rates was estimated to be 0.021 and 0.158 Ma, respectively (Table 5). Analysis of recent demographic bottleneck, based on allele frequencies and heterozygosity, are shown in Table 6 and Supplementary Figure 4. The IAM, TPM, and SMM under the sign test showed heterozygosity deficiencies in 1, 8, and 12 loci, respectively. Both the sign test and Wilcoxon sign-rank test under the IAM showed a significant heterozygosity excess in the metapopulation (p < 0.05). The sign test under the SMM showed a significant heterozygosity excess; however, the Wilcoxon sign-rank test with stronger statistical power did not show a significant value for population bottleneck. The outcome for TPM, a model thought to be most consistent with microsatellite evolution, supported the absence of recent bottleneck in the Jankowski’s Bunting population, and the null hypothesis that the population was in mutation-drift equilibrium. The distribution patterns of allele frequencies can also be used to distinguish demographic stability from bottleneck events. Mode shift analysis showed a normal L-shaped distribution of allele frequencies (Supplementary Figure 4), indicating the absence of population bottlenecks.

The DIYABC analysis produced reliable results (Figures 5A, B). The pre-evaluation of scenario and parameter prior showed that the summary statistics of the observed dataset fell within the core distribution range of the simulated summary statistics under all scenarios (Supplementary Figure 5). A limited overlap between the 95% confidence interval of posterior probability for Scenario 6 and those for Scenario 3 was observed under direct estimation. Nevertheless, the posterior probability for Scenario 6 under direct estimate and logistic regression was 73.20 and 68.45%, respectively, which were much higher than the other scenarios (Figure 5B). In addition, the Type I error rate for Scenario 6 was relatively high (0.77). If Scenario 6 had been frequently selected when it was not the true scenario for simulated dataset (high Type II error rate), we could have little confidence in its selection for the observed dataset. However, the Type II error rate for Scenario 6 was consistently low (0.01–0.09), with an average of 0.05 (Supplementary Table 4), which provides great confidence that when Scenario 6 is chosen it is likely to be the truest of the proposed scenarios. Scenarios 1 and 7 had higher Type I error rate than Scenario 6, showing that they have lower power of detection. Scenarios 2, 3, 4, and 5 all showed lower Type I error rates than Scenario 6 but consistently high Type II error rates; thus, it could not be accepted with confidence if it had been selected as the true scenario. To sum up, Scenario 6, i.e., an old bottleneck followed by a recent expansion, was considered as the most likely scenario for the demographic dynamics of Jankowski’s Bunting. Posterior distributions of the demographic parameters drawn for Scenario 6 along with model checking are shown in Supplementary Figures 6, 7, respectively. All parameter estimates showed high RMedAD values (>0.2) and cannot be considered fully reliable (Supplementary Table 5). Our ABC demographic analysis estimated the initial N_e of the Jankowski’s Bunting population under Scenario 6 to be 8,270 (95% CI: 244–26,800; N4, see Supplementary Table 5). Results suggested that there was a bottleneck 5,746 (95% CI: 1,217-14,059) years ago (t2), resulting in the N_e decrease to 4,470 (95% CI: 755–8,650; N3). Subsequently, the N_e remained relatively stable until the expansion 238 (95% CI: 82–335) years ago (t1). Although the estimated mean value of N1 (3,820) was slightly lower than that before expansion (N3), the actual value of N1 was still larger than that of N3 considering the high overlap of 95% confidence intervals and the priori constraint (N3 < N1) applied to the historical model.

Discussion

Unexpected high genetic diversity is found in threatened populations that have experienced recent devastating declines, prompting a conflict with the subconscious view under the neutral model of DNA evolution that small populations typically exhibit low genetic diversity (Kimura, 1983; Lippe et al., 2006; Torres-Florez et al., 2014; Zhang et al., 2015; Hoffmann et al., 2017). In wild populations, factors affecting genetic variation, such as population structure, demographic history, and inbreeding, often lead to the deviation of polymorphism patterns from neutral expectations and ultimately manifest as a lack of clear correspondence between genetic diversity and population size (Amos and Harwood, 1998). Although the Jankowski’s Bunting population is thought to have experienced a substantial decline over the past few decades, the present results reconfirm the high potential polymorphism in neutral regions of the genome in the recent past (Wang et al., 2011; Li et al., 2021; Huang et al., 2022). This “mismatch” pattern of population size and genetic diversity has been attributed to moderate historical dynamics, extensive gene flow between populations and time-lag effect of genetic response in many species, such as the copper redhorse (Moxostoma hubbsi), California sea lions (Zalophus californianus), Trinidadian guppies (Poecilia reticulata) and Chilean blue whales (Balaenoptera musculus; Lippe et al., 2006; Suárez et al., 2010; Torres-Florez et al., 2014; Fitzpatrick et al., 2016).

Demographic inferences from this study support our previous hypothesis (1) that large historical N_e and moderate demography contribute to rich contemporary genetic variation. The long-term N_e was estimated to range from 800 to 20,818, much higher than the contemporary estimate of 144 (LD method), indicating that Jankowski’s Bunting has experienced pronounced demographic decline. Nevertheless, no signals of population bottleneck were detected by classical estimators. Both the heterozygosity excess test and the mode-shift test prompted us to reject the hypothesis that the population has undergone a sharp decline. These suggested that the population decline may be moderate, and therefore less likely to be imprinted in the genome. Despite the lack of historical samples for genetic comparisons, the current high genetic diversity implied a limited genetic impact of population decline. A similar case is that the widespread New Zealand endemic Notiomystis cincta declined to a single remnant population following European colonization, however, high levels of historical population size, long generation times, random mating and absence of bottlenecks combined to maintain high genetic diversity (Brekke et al., 2011). We found an ancient demographic expansion based on the mismatch distribution of nucleotides and neutrality tests, with an estimated expansion time of 0.021 or 0.158 Ma. Uncertainty in mutation rate severely hindered the estimates of expansion time, but it is certain that this expansion occurred between the end of the penultimate glacial period (c. 0.17 Ma) and the last glacial period (c. 0.115 to c. 0.011 Ma). A similar pattern of population expansion since the penultimate glacial period was also found in the greater horseshoe bat (Rhinolophus ferrumequinum) in Northeast Asia, the black-spotted frog (Pelophylax nigromaculata) and the Chinese Hwamei (Leucodioptron canorum canorum) in East Asia (Zhang et al., 2008; Li et al., 2009; Liu et al., 2016). Among them, the Korean and Japanese populations of the greater horseshoe bat in Northeast Asia expanded at 0.12 Ma (Eemian interglaciation) and 0.15 Ma (Saale glaciation), respectively (Liu et al., 2016). Actually, the most pronounced pattern observed in many species was a severe reduction in N_e coinciding with the beginning of glacial period (Nadachowska-Brzyska et al., 2015). The different demographic pattern we observed may be related to the distinct geological history of Northeast Asia. The expansion time we inferred span the last interglacial period (c. 0.13 to 0.115 Ma) accompanied by global warming (Montoya, 2014). During this period, Northeast Asia was almost ice-free, and the relatively mild climate provided the necessary conditions for population expansion (Wei et al., 2006). This scenario may have shaped the demographic processes and geographic distribution patterns of native species, which differ from those of European and North American species (Wang et al., 2013). Existence of glacial refuges in Northeast Asia may explain the moderate demography of Jankowski’s Bunting (Ye et al., 2019; Lee et al., 2021). In addition, although the climate has gradually cooled since the last glacial period, the grasslands in Northeast China have expanded rapidly in the last glacial maximum, which may provide additional suitable habitat for Jankowski’s Bunting (Li X. et al., 2019). Shrubs and herbs, especially the genus Stipa, became the dominant vegetation in Inner Mongolia during this period (Sun, 1992), which may be related to the population stability and habitat preference of Jankowski’s Bunting. The stronger drought tolerance compared with forest birds may have improved the fitness of Jankowski’s Bunting during the late Pleistocene, which contributes to the similar demographic response with the Chinese Hwamei to the cold period (stadials; Li et al., 2009).

The ABC analysis provided new insights into our inferences on the recent demographic history. The optimal model revealed a historical bottleneck, but the confidence interval for the expansion time (t2) was rather wide. We focused our attention on population dynamics over the last hundred years since the Industrial Revolution. Historical records showed that the habitat of Jankowski’s Bunting has been severely damaged by fire, farmland reclamation, grazing, and development projects over the past few decades (Jiang et al., 2008; Wang et al., 2011). In addition, parasitism and nest predation, which seriously impede reproductive success, have also been observed in field investigations (Zhang et al., 2020). However, consistent with the classical bottleneck tests, we still found no recent declines during periods of dramatic habitat change (t1; 95% CI: 82–335 years ago). The constant N_e in the simulated scenarios is not always interpreted as population stability, as the population decline may be so slow that it cannot become a bottleneck. Given the significant difference between historical and contemporary N_e, the probability of population size remaining stable in the recent past appears to be quite low or, more likely, that the recent population decline could have occurred gradually. After the expansion of Jankowski’s Bunting population at t1, there may have been a slight decline that was not enough to become a bottleneck. Taking the above results together, we believe that a large historical N_e and moderate demographic decline without recent bottlenecks contribute to the formation and maintenance of high genetic diversity in Jankowski’s Bunting. A major concern with inferring demographic factors from genetic data, however, is the time-lag effect of genetic responses to environmental changes in some species (Bailleul et al., 2018; Aavik et al., 2019). The bottleneck-induced heterozygosity excess would persist for 0.2 to 4.0 N_e generations until the equilibrium between mutation and drift is re-established (Luikart and Cornuet, 1998). For the Jankowski’s Bunting population, with an estimated contemporary N_e of 144 (LD method) and a generation time of 1.7 years, a new equilibrium can be established after 28.8 to 576 generations, i.e., 48.96 to 979.2 years. Consequently, any loss of genetic variation may only be detectable after a longer period, which means that the existence of recent bottlenecks cannot be completely denied.

The possible consequences of habitat loss, degradation, and fragmentation caused by human disturbance are genetic isolation, loss of genetic variation, and the extinction of local populations (Lindsay et al., 2008; Segelbacher et al., 2008). Evolutionary processes over historical timescales (earlier than human influence), together with direct effects of habitat change, account for the observed genetic structure (Keyghobadi, 2007; Vandergast et al., 2007; Hansen et al., 2009; Pavlacky et al., 2009). Due to the lack of samples spanning large time scales, we only quantified and compared historical and contemporary gene flow to obtain a dynamic view of population connectivity and genetic structure. The absence of genetic structure, extremely low proportion of interpopulation variation, and high gene flow indicated good connectivity among local populations. Together, the five local populations constitute a large gene pool that may effectively buffer against loss of genetic variation. In addition, the fine-scale geographic analysis suggested that geographic distance is not the dominant factor affecting the current distribution pattern of genetic variation. Considering the high spatial dispersal capacity and wintering dynamics of Jankowski’s Bunting, we suggested that distance-dependent population differentiation caused by habitat fragmentation may be attenuated, or even counteracted, by gene flow between populations. Numerous traditional studies have attempted to reveal the impact of recent human-induced landscape changes by comparing historical and contemporary genetic connectivity (Chiucchi and Gibbs, 2010; Li et al., 2017; Austin et al., 2022). Although our results reveal a relatively high population connectivity, there is a concern that direct comparisons of gene flow across different timescales may lead to the underestimation of changes in genetic connectivity (Samarasin et al., 2017). Therefore, the population connectivity in Jankowski’s Bunting may have experienced a more severe decline over the recent past than is currently observed. Even so, such high contemporary gene flow implied a rather limited impact of habitat change on population connectivity.

During the process of long-term evolution, the combined effects of multiple evolutionary forces, including mutation, natural selection, and gene flow, can account for the current state of genetic variation in species (Jones et al., 2013; Hague and Routman, 2016). As an important complement to genetic resources, gene introgression from the sympatric related species Meadow Bunting (Emberiza cioides) may contribute to the genetic stability of the Jankowski’s Bunting population, although the genetic impact of introgressive hybridization may not be as significant as expected (Huang et al., 2022). In addition, the long-distance dispersal activity of Jankowski’s Bunting during the non-breeding season suggests that unknown parapatric (e.g., Mongolian) populations may also serve as potential genetic reservoirs. Currently, it is difficult to clearly measure and compare the relative effects of various evolutionary forces in the absence of sufficient information. Combined analysis of samples in different years may cause some confusion in the inference of recent demography, e.g., repeated capture of individuals in the same nest may lead to the overestimation of specific genotype frequency. However, compared with the evolution history, the time span of our sampling is obviously quite short, so we can almost regard it as a single time sampling. Our study failed to reveal details of the population history, which may be attributed to the fact that the limited number of molecular markers and sample size did not adequately resolve the information contained in the genome. Furthermore, the high RMedAD value, which represents the low accuracy of parameter estimates, may lead to the confusion of our understanding on true history. Despite these drawbacks, we believe that our study provides valuable insights to understand the recent history of Jankowski’s Bunting. Genome-based Markovian coalescent analysis and single nucleotide polymorphism markers generated by reduced representation sequencing may improve our ability to address these problems.

Although species with short generation times and small N_e are expected to be more sensitive to habitat fragmentation, the time-lag effect on genetic responses cannot be denied (Varvio et al., 1986; Keyghobadi, 2007). In addition, a human-induced bottleneck that occurs over a sufficient time span would result in a significant loss of heterozygosity, while recent anthropogenic effects may not be sufficient to cause observable genetic changes (Nei et al., 1975; Hughes, 2010). Therefore, recent or current population genetic characteristics should be treated with caution in conservation management, as observed levels and spatial patterns of genetic diversity may reflect historic rather than contemporary conditions. Our study highlighted the necessity of temporal genetic monitoring. The contemporary high genetic diversity does not equate to evolutionary health. Factors such as the current state of genetic variation, population history, life history, and ecological traits should be integrated into a systematic framework to generate a robust assessment of population health. In recent years, a series of grassland restoration projects have been implemented that, to some extent, have compensated for the loss of suitable habitat (Han et al., 2017). Even with the high genetic diversity and gene flow, we are not optimistic about the survival of Jankowski’s Bunting populations, considering inbreeding, genetic load, ecological specialization, and the small N_e that do not meet the widely accepted long-term genetic viability target (N_e ≥ 500) (Franklin, 1980; Jamieson and Allendorf, 2012; Laikre et al., 2016). Subsequent conservation efforts should remain focused on the protection of grassland habitats to promote population survival and continuation. Due to high genetic connectivity and absence of genetic structure, the five local populations in this study can be regarded as a genetic unit for regular census and genetic monitoring. Nonbreeding seasons and non-invasive methods, such as feces and feather sampling, are preferred in genetic sampling. Considering that the habitat of Jankowski’s Bunting is relatively fragmented, multiple small reserves should be established to achieve refined protection, but also to construct ecological corridors to reduce inbreeding. Furthermore, to explore the future population dynamics and implement reintroduction (ex situ conservation) programs, it is necessary to use systematic geographical research methods, such as ecological niche modeling, and a combination of historical and contemporary climate information to predict potentially suitable habitats.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Genbank Bank: OP435804-OP435919.

Ethics statement

The animal study was reviewed and approved by Ethics Committee of Jilin Agriculture University.

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 RY helped to collect samples and classify species. GF, WS, and DL helped to perform the microsatellite experiments. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (grant nos. 31601856 and 31670398).

Acknowledgments

The authors would like to thank everyone contributing with samples for this study. Xue Zhang and Yufei Cui are especially acknowledged for their support in the field. We thank Suzanne Leech, from Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript. Lastly, we thank the financial support from the National Natural Science Foundation of China. The content of the manuscript has previously appeared in a preprint online.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.996617/full#supplementary-material

SUPPLEMENTARY FIGURE 1 | Correlation between genetic (Codom genotypic genetic distance matrix) and geographic distances of five local populations inferred by Mantel test.

SUPPLEMENTARY FIGURE 2 | Distribution of differences between historical and contemporary migration rates for population comparisons (m_c-m_h). The contemporary migration rate (m_c) is calculated by BayesAss, while the historical migration rate (m_h) is derived from Migrate-n.

SUPPLEMENTARY FIGURE 3 | Mismatch distribution analyses of Jankowski’s Bunting for mtDNA control region under constant population size model. (A) Metapopulation; (B) ALH population; (C) CG population; (D) GHT population.

SUPPLEMENTARY FIGURE 4 | Mode shift analysis based on allele frequency distribution in Jankowski’s Bunting. (A) Metapopulation; (B) ALH population; (C) CG population; (D) GHT population.

SUPPLEMENTARY FIGURE 5 | Fit of observed data to models defined using ABC analysis. The principal component analysis was generated from 7,000,000 summary statistics from simulations, using prior distributions of historical parameters: open circles are simulated data sets, generated using prior distributions of historical parameters; closed yellow circle is the observed data set.

SUPPLEMENTARY FIGURE 6 | Posterior distributions of the demographic parameters drawn for Scenario 6. N1, recent effective population size; t1, recent population expansion/bottleneck time; t2, ancestral population expansion/bottleneck time; N3, effective population size before expansion; N4, effective population size before bottleneck.

SUPPLEMENTARY FIGURE 7 | Model checking by applying a PCA on the optimal scenario (Scenario 6) in DIYABC analysis. Summary statistics, based on 10,000 simulations, visualize the fit between the simulated and the observed datasets.

SUPPLEMENTARY TABLE 6 | Raw data of microsatellite genotype. The numbers in the table indicate the length of alleles.

Abbreviations

ABC, Approximate Bayesian computation; BSP, Bayesian skyline plot; CR, Control region; HWE, Hardy–Weinberg equilibrium; IAM, Infinite allele model; LD, Linkage disequilibrium; SMM, Stepwise mutation model; TPM, Two-phase mutation model.

References