Changes in Genetic Diversity and Demographic History of Jankowski’s Bunting (Emberiza Jankowskii), as Revealed by Mitochondrial DNA and Microsatellite Markers

The genetic variation and distribution of a population depend largely on the demographic history. For instance, populations that have recently experienced shrinkage usually have a lower genetic diversity. However, some endangered species with a narrow distribution have a high genetic diversity resulting from large historical population sizes and long generation times. In addition, very recent population bottlenecks may not be reected in the population’s genetic information. In this study, we used a mitochondrial DNA marker and 15 microsatellite markers to reveal the genetic diversity, recent changes, inbreeding, and demographic history of a Jankowski’s bunting (Emberiza jankowskii) population in eastern Inner Mongolia. The results show that the genetic diversity of the population remained at a relatively stable and high level until recently. Severe population shrinkage did not result in a considerable lack of genetic variation because of the large historical population size and relatively short periods of human disturbance. In addition, introgression and gene ow among populations compensate for the loss of genetic variation to some extent. Considering the current small effective population size and the existence of inbreeding, we recommend that habitat protection be continued to maximize the genetic diversity of the Jankowski’s bunting population.


Introduction
Genetic diversity and variation are closely related to a species' environmental adaptation and evolutionary potential [1] . Species with high genetic diversity and abundant variation often show strong adaptations to environmental changes [2] . Given the possibility of inbreeding effects and genetic drift, a small population size and narrow distribution usually indicate low genetic diversity [3] , and similar associations have been demonstrated in many previous studies [4,5] . In addition, a population bottleneck is usually accompanied by the loss of genetic variation and the decline of genetic diversity [6] . However, some recent studies have shown that the association is not inevitable and that a large historically effective population size (N e ) and long generation time can act as buffers against the loss of genetic variation [7,8] . Therefore, it is necessary to study species from multiple perspectives, such as population history and species characteristics, to deepen our understanding of the current genetic variation patterns in a population.
Owing to excessive land-use intensity [9] , mainly caused by agricultural production, overgrazing, and other human activities [10] , grassland habitats in Inner Mongolia are threatened by fragmentation and degradation [11] , which greatly affects the survival and reproduction of grassland birds. Jankowski's bunting, Emberiza jankowskii, a resident bird of Eastern Inner Mongolia, has a limited range of several hundred kilometers because of habitat specialization [12] . Historically, Jankowski's bunting was a locally common bird in Northeast China but has undergone a rapid population decline during the last 50 years [13] . Jankowski's bunting is listed as an 'endangered' species on the IUCN Red List of Threatened Species (BirdLife International, 2018). It has not been established whether the population lacks genetic variation as a result of bottlenecks and whether the long-term survival of this species is at risk.
Previous studies on Jankowski's bunting have mainly focused on census, distribution prediction, gut microbiota, reproduction, and nest biology, while population genetic information is very limited [10,11,23] .
Low genetic diversity has been revealed in Jankowski's bunting populations by studying the major histocompatibility complex (MHC class II B) genes, which is an adaptive genetic marker [24] . MHC genes, which are associated with immune responses, are thought to maintain relatively high polymorphisms in natural populations because of heterozygous advantages and frequency-dependent selection [25,26] . The relatively low levels of MHC polymorphism in Jankowski's bunting compared with the common meadow bunting E. cioides indicates that the effects of demographic processes in shaping MHC variation exceed those of natural selection, as well as indicating that the species has a small N e and low resistance to pathogens and parasites. Nevertheless, because of the lack of N e information, the role of demographic processes in shaping genetic variation is not yet understood. The genetic information provided by the adaptive marker is in uenced by both genetic drift and natural selection, which cannot clearly describe the long-term historical population processes [27] . The evolutionary history of populations over long time scales needs to be inferred using neutral markers [7] .
Mitochondrial DNA and nuclear microsatellite markers are considered powerful tools for genetic diversity and demographic history research [28] . Mitochondrial DNA is widely used for studying species and lower taxonomic categories because of the high mutation rate, maternal inheritance, and absence of recombination events [29] . The control region (CR), also known as the D-loop, is considered to have one of the highest evolutionary rates, higher than that of protein-coding genes [30] . Microsatellite markers are especially suitable for endangered species, having advantages such as codominance, selective neutrality, random distribution, low sample requirement, and interspeci c generality [31,32] . Evolutionary forces act differently on the various molecular markers [33] , and a single marker cannot fully re ect the evolutionary process of a population; however, the conjoint analysis of multiple sites can solve this problem [34] . To obtain an accurate view of the current genetic diversity and evolutionary history of the population, it is necessary to compare the genetic variation information obtained using mtDNA and microsatellite markers, respectively [28,31] . Moreover, although single-sample-based population genetics studies provide a snapshot of the current patterns of genetic variation, they lack the ability to reveal population historical processes compared with temporal sampling methods [35,36] . In the absence of historical samples, it is necessary to monitor the short-term changes in genetic diversity, as this helps to attain an accurate picture of the current population status and to predict recent population dynamics, and, at the same time, it is critical for assessing current conservation efforts and developing further targeted management [37] .
In this study, we used one mtDNA fragment (partial mitochondrial CR) and 15 nuclear microsatellite loci as molecular markers to detect recent changes in genetic diversity and inbreeding in a Jankowski's bunting population in eastern Inner Mongolia. The analysis of changes in genetic diversity was based on data obtained from annual sampling. We also reconstructed the demographic history of the species to reveal temporal changes in the N e and to explain the current pattern of genetic variation from a historical perspective. In this report, we discuss the effectiveness of recent protective measures and provide reasonable suggestions.

Changes in genetic diversity
We rst analyzed all the datasets together, followed by the annual data separately, to understand the short-term changes in genetic diversity. We sequenced 549 bp of the CR and, after excluding 132 sites with alignment gaps or missing data, obtained the nal 417 bp of the CR sequence to be used in the genetic analysis. A summary of the CR genetic diversity is presented in Table 1. We de ned 23 haplotypes based on the CR sequence data with 33 polymorphic sites. The haplotype diversity (H d ) was 0.839, the nucleotide diversity (π) was 0.0045, and the average number of nucleotide differences (k) was 1.867. Escµ03, Ecit2, and LS2), of these, the null allele frequency at Mme8 and LS2 reached 0.343 and 0.299, respectively, showing extremely signi cant deviations from the Hardy-Weinberg equilibrium (p < 0.05). In addition, these two loci had low heterozygosity values and were, therefore, removed from subsequent analyses. Considering the in uence of inbreeding, Sosp01, Escµ01, Escµ03, and Ecit2 loci were preserved. In fact, inbreeding, which was evidenced in subsequent analyses, led to the presence of homozygous excess loci and substantially affected the frequency of null alleles. In general, genetic polymorphism at each informative locus is considered low (PIC ≤ 0.25), medium (0.25 < PIC < 0.5), or high (PIC ≥ 0.5) [38] . The polymorphic information content (PIC) values of all loci, except Mme12, were greater than 0.5, with an average of 0.744, suggesting high polymorphism. The mean expected heterozygosity was 0.774, and it varied from 0.434 (Mme12 locus) to 0.952 (Mme2 locus), which were higher than the observed heterozygosity (0.733).
We placed particular emphasis on π and H e , which were not signi cantly correlated with sample size.
There was no substantial decrease in π or H e among the samples from 2012, 2013, 2017, 2018, or 2019 (Tables 1 and 3). In summary, the present study revealed high genetic polymorphism in a Jankowski's bunting population based on both mitochondrial DNA and nuclear microsatellite markers. Although the Jankowski's bunting population has recently diminished, a high level of genetic variation has been preserved.

Inbreeding and effective population size
Relatively close genetic distances usually imply a high interindividual similarity, indicating low intrapopulation differentiation [39] . Under the Kimura 2-parameter model, the mean genetic distance between individuals in the population based on CR data was 0.005 ± 0.001. The genetic distances between pairs of sequences indicated there was a high genetic similarity between individuals within the population and low levels of intra-population differentiation.
An estimated inbreeding coe cient (F is ) of 0.054 indicated the presence of inbreeding within the population. According to the average expected heterozygosity of 0.774, the long-term N e was estimated to be 4645. We used the linkage disequilibrium method to estimate the contemporary N e . The sampling scheme adopted in this study was in accordance with the requirements of the linkage disequilibrium method (> 90 individuals, > 6 loci) [40,41] . This method is less affected by population decline and can avoid the underestimation of N e that occurs in the temporal method because of a small generation interval [42,43] . Because the presence of rare alleles leads to an overestimation of the N e value, we set the minimum allele frequency used in this study at 0.02 to exclude single-copy alleles. On the basis of the linkage disequilibrium method, the N e was estimated to be 106 (95%CI = 94-120). This clearly shows that, compared with the historical population, the contemporary N e was greatly reduced.

Demographic reconstruction
We reconstructed the demographic history of the Jankowski's bunting population by analyzing for recent population expansions and bottlenecks. A recent demographic expansion was inferred from the neutral test and mismatch distribution analysis. In the neutral test results, the Tajima's D and Fu's F s values (Table 1), based on CR data, were signi cantly negative (p < 0.05), indicating demographic expansion has occurred. In the nucleotide mismatch analysis, the Harpending's raggedness index (Hri) and sum of squared deviation (SSD) values implied that the population did not deviate from the null hypothesis, i.e., the population may have experienced a recent expansion. Moreover, the unimodal pattern of mismatch distribution also supported the inference of a recent population expansion (Fig. 1). The population expansion time was estimated at 0.15 Ma.
The bottleneck analysis using allele frequency and heterozygosity are depicted in Table 4 and Fig. 2. The in nite allele model (IAM), two-phase mutation model (TPM), and stepwise mutation model (SMM) under the sign test showed de ciencies in 7, 9, and 18 loci, respectively. Under the IAM, the Wilcoxon sign-rank test revealed a signi cant deviation from the mutation drift equilibrium (p < 0.01), showing there was excess heterozygosity. The SMM showed signi cant heterozygosity de ciency in both tests as a result of its over-sensitivity to microsatellite mutations (p < 0.01). The mutation pattern of microsatellites is considered to be more consistent with the TPM [44] . The outcome for TPM supported the absence of any bottleneck in the E. jankowskii population, and the null hypothesis that the population was in mutationdrift equilibrium was accepted under the TPM. Allele frequency distribution patterns can also be used to distinguish population stability from bottleneck events. The normal L-shaped distribution of allelic frequency (Fig. 2) was revealed by mode-shift analysis, indicating the absence of any recent bottlenecks.

Discussion
In recent years, the genetic variation of many endangered birds has been revealed and utilized effectively in conservation management [45,46] . In general, threatened species usually exhibit low genetic diversity because of their small population size and narrow distribution. However, it seems that this cannot be veri ed in the studied Jankowski's bunting population. According to the present data, the nding of nearly neutral variation suggested there is a high potential of polymorphism in the neutral region at the whole genome level. The haplotype (H d = 0.839) and nucleotides (π = 0.0045) in the CR of the bunting population were more diverse than those of the endangered Elliot's pheasant Syrmaticus ellioti (H d = 0.584, π = 0.0015), and even higher than those of the domesticated/commercial turkey M. gallopavo (H d = 0.569, π = 0.0016), a species of 'least concern' [47,48] . Interestingly, the expected heterozygosity values based on microsatellite analyses were found to be high in a few endangered birds, such as the blackcapped vireo Vireo atricapilla (H e = 0.810) [49] , which is slightly above that of Jankowski's bunting (H e = 0.774). This characteristic of high genetic polymorphism is clearly unusual in endangered birds, as these populations are thought to have undergone substantial contraction over the past few decades.
Despite the lack of census records, Jankowski's bunting is believed to have a large historical population, which is the basis for the high genetic diversity [13] . Historical records show that the habitat of Jankowski's bunting has been severely damaged by re, farmland reclamation, grazing, and development projects [10,11] . In addition, parasitism and nest predation, which seriously impede reproductive success, have also been observed in our eld investigations [50] , but there was no considerable decline in population genetic diversity over nearly a decade in this study. The demographic decline may have only slightly affected the genetic variation, most of which was retained. This scenario was also found in a comparison of historical and contemporary samples of the endangered black-capped vireo V. atricapilla, which has experienced more severe census shrinkage over the last century than Jankowski's bunting. E. cioides, a species of least concern, can be found within the distribution range of Jankowski's bunting. The two species are phenotypically and genetically similar and are clustered within a highly supportive branch of the phylogenetic tree [51] . The presence of identical or similar alleles between the two species has been found in previous studies [24] . Therefore, we hypothesized that gene introgression has compensated for the loss of genetic diversity during population decline to some extent. Gene ow between populations leads to the convergence of genetic composition and the sharing of different alleles among populations. The population studied is an aggregate population composed of several small subpopulations. Close interindividual genetic distances within a population suggest there is su cient gene ow between the smaller populations. Reduction of the amount of suitable habitat available promotes the fusion of multiple populations, which is one of the compensation mechanisms for genetic variation within populations.
In addition, recent habitat changes may not be enough to affect bird species with higher dispersal abilities. Reportedly, when a population decreases, the population heterozygosity remains stable for 0.2 to 4.0 generations until a new equilibrium is re-established [52] . In the Jankowski's bunting population, with an estimated contemporary N e of 106 individuals and a generation time of 1.7 years, a new equilibrium could be established after 21 to 424 generations, or 36 to 721 years. Consequently, a loss of genetic variation may only be detectable after a longer period.
Current patterns of genetic variation provide a snapshot of the demographic history. The estimated contemporary N e (106) was substantially less than the estimated long-term N e (4645). The sharp contrast between the historical and current N e suggests that the population has experienced a gradual decline, but we have almost failed to detect recent bottlenecks. Both the heterozygosity excess test and the modeshift test prompted us to reject the hypothesis that the population had undergone a sharp decline.
Therefore, the reduction occurred slowly or was not severe enough to be detectable. Another possible reason for the high genetic diversity and absence of recent bottlenecks is that the population decline caused by habitat fragmentation and degradation is rather recent. It began only a few decades ago, thus, there has not been enough time for a loss of genetic variation. Conversely, the nucleotide diversity was relatively low compared with the high haplotype diversity, indicating the presence of recent population expansion. This was con rmed by the results of the neutral test and mismatch distribution analysis.
The quaternary climate uctuation is considered to have greatly in uenced the geographical distribution, differentiation, and demographic history of global species. Previous studies have shown that few large glaciers formed in Northeast Asia during the alternation of the last glacial and interglacial periods [53] , and the relatively mild climate provided the necessary conditions for population expansion. In this research, population expansion was estimated to occur approximately 0.15 Ma, during the penultimate glacial period (about 0.17 Ma). An unusual pre-last interglacial (LIG) expansion may be attributed to the lack of glaciers in eastern Inner Mongolia, which provided a refuge for Jankowski's bunting. Palynology and paleobotanical records show that the average annual temperature during the LIG was higher than that of modern times [54,55] . In addition, herbaceous plants, especially of the genus Stipa, have become the dominant vegetation in Inner Mongolia [56] , which may be related to the particular habitat preference of Jankowski's bunting. We speculated that the population continued to expand during the LIG, accumulating a wealth of genetic variation. Similar interglacial expansion patterns have been discovered in several East Asian species, such as the Chinese hwamei Leucodioptron canorum canorum [57] and the green-backed tit Parus monticolus [19] . Subsequently, the global climate cooled and dried leading up to the last glacial maximum (LGM; 0.021 − 0.018 Ma). Many European and North American species experienced population shrinkage during the LGM, followed by expansion [58,59] ; however, we did not detect any sharp population uctuations. Glacial and interglacial shelters played a key role in the long evolutionary history of the species. The development of the modern climate towards drought has led to the deserti cation and fragmentation of grasslands, which further affect the survival of Jankowski's bunting. However, the impact of human interference is undeniable. The short-term changes in genetic diversity provide no evidence of a substantial recent loss of genetic variation. The absence of considerable changes may be attributed to the short sampling interval, but it may also indicate that the population has remained relatively stable in the recent past. Our study failed to reveal details of the population evolution of the species. To explore the future dynamics of populations and implement reintroduction programs, it is necessary to use systematic geographical research methods, such as ecological niche modeling, and combine historical and contemporary climate information to predict potentially suitable habitats.
Although high genetic diversity has been observed in some endangered species with small population sizes and narrow distributions, this does not mean that high genetic diversity is unnecessary for the longterm survival of the population [60] . The absence of a considerable recent decline in genetic diversity suggests that conservation efforts have been effective. In recent years, a series of grassland restoration projects have been implemented, which, to some extent, have compensated for the loss of suitable habitat. A small population size is often associated with inbreeding depression [61] . Considering the small effective population size and the existence of inbreeding, the survival of the Jankowski's bunting population is still under threat. Subsequent conservation efforts should remain focused on the protection of grassland habitats to preserve population survival and continuation and to maximize genetic diversity.

Materials And Methods
Ethics approval and consent to participate 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 nal extension at 72°C for 10 min.
The PCR products were electrophoresed on 1% agarose gels and sequenced by Sangon in Shanghai, China.

SSR ampli cation and genotyping
A total of 15 polymorphic SSR loci and 165 of the 200 samples were used in this study. All the primers for ampli cation were obtained from previously published papers. The speci c sequence and information for each microsatellite primer can be found as Supplementary Table S1 online. We optimized the reaction volumes and PCR protocols for each primer pair. Microsatellite loci were scored by Sangon Biotech (Shanghai, China). Final fragment sizes were calculated in GeneMarker v1.91 (SoftGenetics LLC, State College, PA, USA) and scored manually.
mtDNA data analysis Unreliable sequences at both ends were deleted by comparing the sequencing peak diagrams using Bioedit Sequence Alignment Editor [62] . Clustal X v1.83 was used for multiple sequence alignment [63] .
Standard genetic statistics were performed with universal software. To evaluate genetic diversity, common diversity 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 [64] .
To assess the level of inbreeding in the population, the inbreeding coe cient (F is ) was computed on Arlequin v3.5 [65] . The mean genetic distance between individuals within the population was calculated using MEGA X [66] under a Kimura 2-parameter model. All parameters were evaluated following the bootstrap method and used 1000 sampling repeats to obtain exact p-values.
Neutral tests based on allele frequency distributions were carried out with Arlequin v3.5 to calculate the Tajima's D and Fu's F s values and their signi cance. Tajima's D and Fu's F s , as methods of testing neutral evolution, can be used to trace the demographic history [67,68] . Tajima's D test is more powerful for assessing the longer history of the population, while Fu's F s is sensitive to recent uctuations [69] . Arlequin v3.5 was also used to analyze the mismatch distribution based on the sudden expansion model [70] . The Harpending's raggedness index (Hri) and sum of squared deviation (SSD), which are parameters describing the degree of observed and expected curve tting, were used to evaluate the validity of the model. Subsequently, the mismatch distribution diagram was plotted in DnaSP v6 and was used to visualize the historical dynamics of the population. If an expansion event was detected, the formula τ = 2µkt [70] was used to estimate the population expansion time (t), in which τ(tau) is the expansion constant, µ is the mean mutation rate of nucleotide, and k is the sequence length of the CR. The value of τ was calculated via Arlequin v3.5.

SSR data analysis
The genotype and its frequency distribution were determined according to genotyping results. We tested for deviations from the Hardy-Weinberg Equilibrium (HWE) resulting from heterozygote excess or de ciency at each microsatellite locus and over all loci with Genepop v4.7 [71] . To detect whether the locus was genetically independent, Fisher's exact test of linkage disequilibrium (LD) was implemented in Genepop v4.7 for all pairwise combinations of loci. The exact p-values obtained by the Markov chain method were corrected by sequential Bonferroni tests. Micro-checker v2.2.1 [72] was used to evaluate possible null alleles, large allele dropout, or stuttering for each locus.
Estimations of the number of alleles (N a ), allelic richness (A r ), observed heterozygosity (H o ), and expected heterozygosity(H e ) for each studied locus and overall loci were performed in FSTAT v2.9.4 [73] . Cervus was used to calculate the value of the polymorphic information content (PIC), which is a parameter used to describe the genetic information richness contained in the site [74] .
The calculation of the inbreeding coe cient (F is ) based on genotype data was implemented in Genepop v4.7. The long-term effective population size (N e ) was estimated using the formula N e = (1/[1 − H E ] 2 − 1)/8µ, where H E is the average heterozygosity and µ is the mutation rate of 5 × 10 − 4 commonly used for microsatellite loci [75] . We also estimated the contemporary N e based on the random mating model by NeEstimater v2 [76] . The linkage disequilibrium method and the heterozygosity excess method were based on a single sample, while the standard temporal method was based on multiple samples. Considering the characteristics of the samples used in this study, the rst estimators were used for analysis.
We used two methods based on allele frequency distribution and heterozygosity to detect the demographic bottleneck events. Under the in nite allele model (IAM) [77] , stepwise mutation model (SMM) [78] , and the two-phase mutation model (TPM) [44] , recent population genetic bottlenecks were tested by Bottleneck v1.2.02 [79] . TPM tests were run assuming 70% single-step mutations with 30% multi-step mutations. In general, the sign test, standardized differences test, and a Wilcoxon sign-rank test are the three most commonly used methods to determine signi cantly excess heterozygosity. We performed the sign test and Wilcoxon sign-rank test to analyze heterozygosity excess or de ciency, as they were statistically capable and appropriate for this study. Using this method to infer bottleneck events does not require the prior acquisition of population size and genetic variation information. Additionally, we also checked for deviations of the allele frequency distribution from the normal L-shaped distribution (mode shift). Rare alleles tend to be lost in populations experiencing bottlenecks, resulting in fewer lowfrequency alleles than intermediate-frequency alleles. A mode-shift indicator, represented graphically, can detect the loss of rare alleles in a bottleneck population. All tests were based on 1000 iterations.

Data Availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.