- Research article
- Open Access
Deep south-north genetic divergence in Godlewski’s bunting (Emberiza godlewskii) related to uplift of the Qinghai-Tibet Plateau and habitat preferences
BMC Evolutionary Biology volume 19, Article number: 161 (2019)
Geological events and climatic changes played important roles in shaping population differentiation and distribution within species. In China, populations in many species have contracted and expanded responding to environmental changes with the uplift of the Qinghai-Tibet Plateau (QTP) and glacial cycles during Pleistocene. In this study, we analysed the population structure of Godlewski’s Bunting, Emberiza godlewskii, to determine the effects of major historical events, geographic barriers and past climatic changes on phylogenetic divergence and historical demographic dynamics of this species.
A phylogeny based on concatenated mitochondrial and nuclear DNA datasets show two (northern and southern) clades approximately diverged 3.26 million years ago (Ma). The West Qinling Mountains serve as a dividing line between the two lineages. Both lineages experienced a recent demographic expansion during interglacial periods (marine isotope stages (MISs) 2–6). Bayesian skyline plots and the results of ecological niche modelling suggested a more intensive expansion of the northern lineage during the late Pleistocene, whereas the southern lineage was comparatively mild in population growth.
Our results provide insights into the distribution patterns of avian taxa and the possible mechanisms for a south and north divergence model in China. The deep divergence may have been shaped by the uplift of the QTP. Habitat preferences might have facilitated the lineage divergence for E. godlewskii. Moreover, the West Qinling Mountains act as a dividing line between the two lineages, indicating a novel phylogeographic pattern of organisms in China. The difference in population expansion mode between two lineages resulted from different effects caused by the climate of the LGM and the subsequent habitat changes accompanying the arrival of a colder climate in northern and southern regions of China.
The phylogeographic patterns of avian taxa in China are diverse . Among them, the south-north divergence is not quite common for widely-distributed birds. For example, studies on generalist species such as the Black-billed Magpie, Pica pica, which inhabit various open and semi-open habitats including short-grass areas but depend on trees and tall bushes, did not find major divergence between northern and southern populations . In contrast, the azure-winged magpie, Cyanopica cyanus, which requires well-wooded habitat, shows a north-south genetic division in eastern China . Multiple factors contribute to the observed phylogeographic patterns of divergence in China, including environmental features , geographical barriers , glacial periods , habitat preferences of organisms . It is important to combine geo-climatical and ecological components in interpreting phylogeographical patterns we observed from empirical studies.
Zoogeographical studies have shown that the Himalayan Mountains and the southwestern mountain systems in China are important avian biodiversity hotspots in the Northern Hemisphere  and the Qinling Mountains – Huai River has been regarded as an ecological and zoogeographical boundary between the Oriental and Palaearctic realms in China . However, the roles of mountain system in southwest China are not substantially illustrated (but see literatures by Song et al.  and Qu et.al. ), and few study to demonstrated the barrier effects of Qinling Mountains-Huai River line on lineage diversification in birds but questioned by Song et al. . Therefore, it is necessary to evaluate the phylogeographical impacts of these geographic features and their related geological events on birds in central and south China.
Old World buntings (genus Emberiza) are widely distributed across the Palearctic, from East Asia and the Himalayas to the Middle East and Africa . The species occupy a variety of breeding habitats. Most prefer semi-open to open taiga forest habitats, steppes and forest edges to very dense forests [12, 13]. Godlewski’s Bunting, Emberiza godlewskii, shows distribution and habitat preferences typical of the genus. The species tends to select bushy and rocky hill slopes that are often near forests, thickets, ravines, and farm fields . In China, Godlewski’s Bunting is mainly distributed in the foothills of the Tien Shan range and on the western rim of the Tarim Basin. The species is also distributed throughout Gansu and Qinghai Provinces, northern and western Sichuan Province, Tibet (mainly on the northern side of the Himalayan range), northern and western Yunnan Province, Guizhou and Guangxi Provinces, and the mountains from the north of Shanxi-Beijing-Hebei to the south of the Qinling range in China [12, 14]. Five [15, 16] or six subspecies  have been recognized according to geographical distribution and morphological characteristics. Cheng et al.  proposed the existence of a subspecies styani based on differences in morphology and allopatric distribution from other subspecies. In our research, we included six subspecies and considered styani have a special status.
Although numerous species diversified in the Old World, there are few studies on the phylogeography of Emberiza. Päckert et al.  studied the phylogeny of buntings and proposed that elevational segregation between alpine and lowland breeding habitats might be responsible for the northern and southern lineage separation between E. g. yunnanensis and its northern conspecifics. Nevertheless, due to limited geographic coverage and a small sample size, the genetic structure and population history of E. godlewskii is not clear. In addition, previous studies failed to provide information on detailed demographic changes during glacial-interglacial cycle and distribution shift of the species caused by suitable habitats change during the last interglacial (LIG) and Last Glacial Maximum (LGM).
In this study, we incorporate an intensive sampling coverage and integrative approaches combining phylogeography with an isolation-with-migration (IM) model  to study the range-wide phylogeographic pattern of E. godlewskii, and investigate demographic dynamics with hypothetical alternative scenarios. Our specific objectives are as follows: (1) To identify the factors driving the phylogeny of E. godlewskii across its distribution range, locate the divergence boundary and determine whether apparent gene flow and migration existed between the major lineages; (2) To reconstruct a detailed demographic history of the major lineages with approximate Bayesian computation (ABC) modeling  in order to evaluate the effects of Quaternary climatic changes, use ecological niche models (ENMs)  to explore whether genetic differentiation within regions occurred and whether the potential distribution of the bunting contracted into refugia during the LGM in response to suitable habitat changes during Pleistocene glacial-interglacial cycles, and examine the fragmentation (pre-Pleistocene ‘vicariance’ model) between the major populations in which genetic differentiation occurred according to a null model of isolation without migration.
Sample collection and laboratory work
We collected 190 blood samples from 26 sites covering the major distribution of Godlewski’s Bunting in China (see Fig. 1, Additional file 1: Table S1). Individuals were caught by mist net and blood samples were extracted from the vein under wing. All individuals were processed in the field and released immediately thereafter. The samples were preserved in anhydrous ethanol at − 20 °C and stored in the Zoological Laboratory of Lanzhou University.
Total genomic DNA was extracted using a TIANamp Genomic DNA Kit (Tiangen, Beijing, China). Cytochrome oxidase I (COI), cytochrome b (Cytb), the control region (CR) and two nuclear gene introns (beta-fibrinogen intron 5, FIB; muscle-specific kinase intron, MUSK) were amplified by polymerase chain reaction (PCR). The primer pairs and PCR conditions with optimized annealing temperatures are described in Additional file 1: Table S2. The PCR products were purified with a StarPrep Gel Extraction Kit and examined with 1% agarose gel electrophoresis. Sequencing was performed by the Sangon Biotech Company in Shanghai. The original sequences were assembled with SeqMan software (DNASTAR). The consensus sequences were aligned using ClustalW and manually proof read in MEGA.
Genetic diversity and phylogenetic structure analyses
We calculated the genetic parameters of the population, including the number of segregating sites (S), the number of haplotypes (Nhap), haplotype diversity (Hd), and nucleotide diversity (Pi) in DnaSP 5.0. The heterozygous sequences of the two nuclear genes were phased using the program PHASE  with the default parameters.
We constructed the intraspecific phylogeny based on the concatenated mitochondrial haplotypes using a constant population size coalescent model in BEAST 1.8 . The optimum model setting was determined by the results of MrModeltest 2.3 . Because the sequence substitutions were expected to be constant at the intraspecific level, we applied the strict molecular clock model . We ran 100 million generations for the mitochondrial DNA (mtDNA) data set of the Markov chain Monte Carlo (MCMC) chain, with a sampling frequency of 10000 generations. Convergence of the posterior distributions of the parameters was evaluated by monitoring the effective sample size (ESS > 200) and trace plots in Tracer1.5 . A maximum-credibility tree, which represents the topology of maximum posterior probability, was calculated in TreeAnnotator 1.8  after discarding the first 10% of trees as burn-in.
To test the validity of the major lineages revealed in the mtDNA, a species tree was estimated by *BEAST  based on mtDNA and nucler DNA (nuDNA) of 101 individuals. The *BEAST analysis was implemented in BEAST 2.0 . We unlinked the substitution models across the five genes (Cytb, COI, CR, FIB and MUSK) and set the substitution parameters for each locus according to the MrModeltest results. We selected a piecewise linear and constant root model as the prior species tree and employed default molecular clock settings (strict clock model, rate of the first locus = 1.0, estimated rates of the remaining loci). The MCMC chains operated for 500 million generations. The sampling frequency was every 50000 generations, and the first 10% was discarded as burn-in. Convergence of the posterior distribution parameters was examined by monitoring the ESS (> 200) and trace plots in Tracer1.5 . The results were visualized with DensiTree software . We also reconstructed phylogenetic trees with MrBayes  and RAxML  based on concatenated fragments (mtDNA and nuDNA). The best-fitting model for the BI and RAxML analyses was determined separately for each gene at the five loci by MrModeltest 2.3 . We included sequences of Emberiza rutila and Emberiza leucocephalos the outgroups. We conducted an analysis of molecular variance (AMOVA) based on FST to further determine where the major divergence occurred. We used AMOVA based on the concatenated mitochondrial set in Arlequin 3.5  to diagnose population structure. As the sample sizes of some the populations, such as Lanzhou (LAN), Lianhuashan (LHS), Changdu (CD), Mangkang (MK), Deqin (DQ) and Derong (DR) were too small, we combined the LAN and LHS populations into LAN/LHS, the Hebei (HB) and Beijing (BJ) populations into HB/BJ, the CD and MK populations into CD/MK, the Lhasa (LS) and Xiongse (XS) populations into LS/XS, and the DQ and DR populations into DQ/DR according to the corrected average pairwise differences [pi XY- (pi X + pi Y)/2] and geographic conditions. The significance of each of the various components in the AMOVA was tested with 10000 permutations. Moreover, to verify the consistency of topological structure, we conducted clustering analyses of the median-joining networks with mtDNA to infer the intraspecific relationship among haplotypes using NETWORK 5.0 software .
Divergence dating and historical demographic reconstruction
According to the phylogenetic results, we calculated the divergence time and reconstructed the demographic dynamics of the major lineages. Divergence times were estimated in BEAST 1.8  based on the concatenated mitochondrial data set. We accepted the substitution rate of 0.01035 per site per million years (molecular clock for the passerine mitochondrial Cytb gene)  and modified it for the combined mtDNA sequences. We calculated the average overall mean p distances of Cytb and the combined mtDNA. The average overall mean p distances of Cytb and the combined mtDNA were 0.02100 and 0.02140, respectively, which implied a substitution rate of combined mtDNA 1.0190 times greater than that of Cytb. Therefore, we obtained a substitution rate of 0.01055 per site per million years for the combined mtDNA sequences . The analyses were implemented with a constant growth as a coalescent constant tree prior and a strict molecular clock. The MCMC chains were run for 100 million generations with sampling every 10000 generations. We used Tracer 1.5 to examine the posterior distribution and ESSs and then summed the posterior probabilities of each parameter after discarding the first 10% of the samples as burn-in.
We calculated Tajima’s D , Fu’s Fs  and Ramos-Onsins’s and Rozas’s (R2) statistic  for each lineage to assess the neutral evolution of the mitochondrial genes. In Arlequin 3.5, we tested for historical demographic expansion using mismatch distribution analyses, and statistical significance was tested with two important parameters: the sum of squared deviations (SSD) and the raggedness index (r) between the observed and expected values . The P values for these parameters were obtained from 1000 bootstrap samples.
To explore whether the major phylogenetic lineages experienced historical fluctuation, we conducted Bayesian skyline plots (BSPs) using mtDNA in BEAST 1.8 to reconstruct the changes in effective population size through time since the most recent common ancestor (TMRCA). We ran the analysis for 100 million generations or more until the ESSs were larger than 200, with sampling every 10000 generations after discarding the first 10% of samples as burn-in. Population curves of demographic history through time were reconstructed in Tracer 1.5.
To explore the demographic scenarios of divergence and expansion, we also applied approximate Bayesian computation (ABC) analysis in DIY-ABC, version 2.1.0.  with the mtDNA and nuDNA data sets. We adopted the divergence with isolation model based on the phylogenetic structure, which clearly displayed two lineages diverged from a common ancestor, and tested the possible scenarios of expansion for each lineage (Additional file 1: Figure S1). For the prior settings for parameters (Additional file 1: Table S3) in DIY-ABC, we used two groups: group 1, including 1.00E-9 to 1.00E-7 per site per year uniform priors, and group 2, including 1.00E-11 to 1.00E-8 per site per year uniform priors, for mtDNA and nuDNA respectively. We tested the following five scenarios of demographic changes for each lineage: a constant population size since divergence (scenario 1); recent expansion (scenario 2); old expansion (scenario 3); expansion-shrinkage (scenario 4); and expansion-shrinkage-expansion (scenario 5)  (Additional file 1: Figure S2a). To select the best model that could explain the genetic polymorphism, we simulated 1 000 000 multilocus genetic data sets for each scenario. The 1% of the simulated data sets closest to the observed data was used to estimate the relative posterior probability (with 95% confidence intervals (CIs)) of each scenario via logistic regression and posterior parameter distributions.
To better estimate gene flow and divergence time between the major lineages, we used IMa2 based on the IM model  with nuDNA and mtDNA data. In the analysis, demographic parameters, including divergence time (t), effective population size of each extant population (q0 and q1) and the ancestor (q2), and migration rates (m0 and m1) between two groups, scaled by the mutation rate (μ), were estimated. The maximum priors for the parameters used for IMa analyses were t = 0.48, q0 = q1 = q2 = 0.084, and m0 = m1 = 8.3. The inheritance scalars were set to 0.25 and 1 for mtDNA and nuDNA respectively. We performed multiple runs to ensure that posterior probability distributions converged by monitoring ESS values and trend lines along MCMC chains. In the simulations, we conducted each IMa simulation for 100000 steps with a burn-in of 100000 steps. The mtDNA and FIB loci were assumed to follow a Hasegawa-Kishino-Yano (HKY) model, and the MUSK locus was assumed to follow an infinite sites (IS) model of mutation, with a mean mutation rate value of 1.6E-8, 1.035E-8, 2.4E-8, 1.2E-9 and 1.9E-9 and confidence intervals of 1.0E-8 to 3.0E-8, 1.0E-8 to 3.0E-8, 1.0E-8 to 4.0E-8, 7.0E-10 to 2.0E-9 and 1.0E-9 to 3.0E-9 for COI, Cytb, CR, MUSK and FIB, respectively.
Distribution changes with ecological niche modelling
To explore whether genetic differentiation within regions occurred in response to distribution changes of suitable habitat during Pleistocene glacial-interglacial cycles, we employed ENMs to predict the distribution of the E. godlewskii with 19 bio-climate variables from WorldClim  during four periods, namely, the present time (1960--1990 AD), the Mid-Holocene (MIH, 0.006 Ma), the LGM (0.021--0.018 Ma) and the LIG (0.14--0.12 Ma). We collected 766 occurrence records from the Global Biodiversity Information Facility (GBIF) (www.gbif.com), eBird, and our sampling sites, and after quality control steps to reduce spatial auto-correlation and sampling bias, we obtained 253 locations. We performed ENM analyses of the two different lineages based on the phylogenetic results, with 126 sites for the northern lineage and 127 sites for the southern lineage, with at least one of the longitude and latitude differences greater than 0.1° between every pair of locations. We constructed the model by randomly selecting 80% of the occurrence data, and the remaining 20% of the data were used to test the accuracy of the model in MaxEnt v3.3.3 k. The parameters were set to the default convergence threshold (10--5), 2000 maximum iterations and 10 replicates. MaxEnt was used to test the model by calculating the average area over ten replications under the receiver operating characteristic curve (AUC) and the binominal probabilities indicating the predictive ability of the model. Then, we reclassified the results into binary states under the 10% logistic threshold, indicating an unsuitable area, and higher than the threshold, indicating a suitable area.
Genetic diversity and molecular evolution
We obtained 1103 base pairs (bp), 1323 bp, and 1094 bp of Cytb, COI and CR, respectively, for 190 individuals. For the nuclear gene introns, we obtained 509 bp and 558 bp of FIB and MUSK, repectively, for 101 individuals. All sequences were deposited in GenBank (see Additional file 1: Tables S4 and S5). The Cytb, COI, CR, and COI + Cytb+CR data sets contained 101, 103, 100 and 304 variable sites, defining 72, 75, 89 and 144 haplotypes, respectively. MUSK and FIB contained 16 and 62 variable sites, defining 18 and 70 haplotypes, respectively. CR had the highest Hd (0.968) and Pi (0.024) values among the mitochondrial genes, whereas FIB had larger Hd (0.960) and pi (0.018) values than MUSK. Neither Tajima’s D test nor Fu’s Fs indicated significant deviations from neutrality at the five loci, with the exception of a significantly negative Tajima’s D value (P < 0.05) for the MUSK gene (Table 1).
Intraspecific phylogeny and geographic divergence
The best molecular evolutionary model of the combined mtDNA was GTR + I + G. The best-fitting models of COI, Cytb, CR, MUSK and FIB were GTR + I, HKY + I, GTR + I + G, HKY and HKY + I, respectively. The intraspecific phylogeny based on mtDNA strongly supported deep divergence of the northern (clade A) and southern (clade B) clades in China (Fig. 2). There were no subclades corresponding to geographical subspecies. The RAxML, BI and *BEAST analyses recovered intraspecific phylogenies supported by high posterior probabilities and bootstrap values and that were consistent with the phylogenies from the mtDNA data (Fig. 3). The median-joining network of the mtDNA was topologically congruent with the trees (Figs. 2 and 3).
The results of the AMOVA based on the concatenated mtDNA showed that the maximum percentage (93.35%) of molecular variance resulted from between groups when the whole population localities were divided into northern and southern lineages. The variation percentages among the localities within group and within the populations were 0.89 and 5.75%, respectively (Table 2), and all of the P values were significant.
Divergence dating and historical demographic reconstruction
The divergence time between clade A and clade B was dated to 3.26 (95% highest posterior density (HPD): 2.53–4.08) Ma, which was approximately in the late Pliocene (Table 3). The TMRCA of clade A was dated to 0.32 Ma (95% HPD: 0.21–0.38 Ma), and the TMRCA of clade B was dated to 0.27 Ma (95% HPD: 0.18–0.32 Ma). The divergence time estimated by IMa2 also showed divergence between the two lineages during approximately the late Pliocene. The maximum likelihood estimates of 2 Nm from north to south and vice versa were 0.091 and 0.077 migrants per generation, respectively (Table 4).
The BSPs showed different demographic trends between the two lineages: The northern lineage showed rapid expansion, whereas the trend curve of the southern lineage was shallower (Fig. 4). Past population dynamics indicated that the rapid population growth of the northern lineage occurred 0.05 Ma whereas a mild expansion of the southern lineage occurred 0.12 Ma. The results of ABC (Table 5) for the divergence with isolation model indicated that the effective population sizes of the two groups were 6.97E+ 05 (95% CI: 4.15e+ 05 to 9.42e+ 05) and 6.23E+ 05 (95% CI: 3.41e+ 05 to 9.11e+ 05) and divergence time was 1.06E+ 06 (95% CI: 3.26e+ 05 to 4.42e+ 06) years ago. The logistic regression of the 1% of simulated data sets are similar to the observed data showed, the recent expansion (scenario 2) of both the northern and southern lineages fit best (Additional file 1: Figure S2b), indicating that both the northern and southern lineages experienced recent expansion (Fig. 5, Additional file 1: Table S6).Tajima’s D and Fu’s Fs values showed significant deviation from neutral evolution for both the northern (P values for Tajima’s D and Fu’s Fs: 0.001and 0.000, respectively) and southern clades (P values for Tajima’s D and Fu’s Fs: 0.003 and 0.000, respectively), and the values of R2 were small and positive (0.0309 and 0.0310, respectively) for both the northern and southern clades, indicating an excess of low-frequency mutations relative to expectations under the standard neutral model. The mismatch distributions of the northern (P values for the SSD and r: 0.430 and 0.370, respectively) and southern clades (P values for SSD and r: 0.740 and 0.880, respectively) fit the population expansion model, although the northern mismatch curve had two peaks. These results supported the model of historical demographic expansion (Table 6 and Fig. 6).
The ENM results (Fig. 7) showed that the LGM distribution was not consistent with the LIG distribution; in particular, the LGM distribution was different from the LIG and MIH/CURRENT distributions, suggesting that the distribution of the southern lineage experienced partial contraction and expansion during the LGM (25,000–15,000 years ago) and post-LGM, while the northern lineage contracted into North China from the LIG to the LGM and weakly expanded to the east and west after LGM.
Deep north-south lineage divergence in E. godlewskii
The phylogenies showed that the E. godlewskii populations were divided into northern and southern clades. Each clade was further differentiated into subclades, but did not fully correspond to the subspecies. These findings indicate that subspecies of E. godlewskii cannot be genetically classified into simple groups and need re-evaluation.
As previously indicated by Päckert et al. , E. godlewskii was deeply divided into two lineages in China. The northern lineage was widely distributed in the northwest (Gansu, Xinjiang, Ningxia and Shaanxi) and north (Shanxi, Hebei and Beijing) of China, which are arid and semi-arid areas of the temperate zone. The southern lineage was distributed in the southeast of the Qinghai-Tibet Plateau (QTP) and the southwest (Yunnan, Guizhou, Guangxi and Sichuan) of China, corresponding to semi-humid areas of the subtropical zone. The belt in the narrow range from LHS to SP seems to act as the line dividing between the two lineages. The fragmentation of species distributions seems an important cause of allopatric divergence. Our ENM results indicated that there is a geographic boundary between the northern and southern clades. Moreover, IMa2 analysis confirmed a nearly null model of isolation without gene flow between the diverging populations, a gap between northern and southern populations. According to the divergence time estimated by the BEAST, the genetic differentiation might have been formed before the Pleistocene in a ‘vicariance’ model.
We propose that the deep divergence may be attributed to the uplift of the QTP, and a boundary corresponds to the West Qinling Orogenic Belt, located approximately 96°-106° E and 33°-37° N. The subduction collision and orogeny of this belt have been dated to the Silurian, but rapid uplift occurred in the late Pliocene . Phylogenetic studies have revealed that some divergent events were consistent with the QTP uplift. Bao et al.  showed that the Tibetan partridge, Perdix hodgsoniae, diverged from its ancestor approximately 3.6 Ma, corresponding to the early intensive uplift of the QTP. The changes in climate, ecology, and habitats following the uplift of the QTP during the Pleistocene were likely the most important factors driving the pattern of phylogeography in Asian birds . The difference in demographic dynamics and distribution range shifts between north and south lineages indicates that the specific habitat requirements of E. godlewskii was another important factor driving its expansion and shaping its distribution. The patterns of the south-north lineages are at least partly explained by ecological differences, especially in habitat preference.
The distribution pattern of E. godlewskii may be related to three geological events that occurred at different timescales in China. The first event was the early intensive uplift of the QTP, which may have been the major driver of the deep north-south split at 3.26 Ma, corresponding to the late Pliocene to early Pleistocene. The uplift of the QTP led to the geomorphological process of major downcutting in the Hengduan Mountain System , including the Yunnan-Guizhou Plateau , and shearing basins and graben valleys formed in the Southwest Mountains from the late Tertiary to the early Quaternary [46, 47]. These rugged basins or valleys at the edge of forests would have become refugia for E. godlewskii during glaciation. The second event consisted of ecological and climatic changes. During the uplift and the following climate fluctuations, the southeastern flanks of the Himalayan and Hengduan Mountains were part of the monsoon realm . These areas were influenced by the southwestern monsoon winds that travelled deep into the continent through the north-south-oriented mountain ranges in Yunnan from the Indian Ocean . Most of the areas were covered by subtropical mixed forest comprising evergreen broadleaf and broadleaf deciduous woodlands . The dense forests separated the lowland and alpine-shrubland habitats at the forest edge into patches. At this time, in the intermontane basins of the Yunnan-Guizhou Plateau, the climate was dry and warm and the dominant vegetation in the late Pliocene was shrub meadow . These conditions met the habitat requirements of E. godlewskii. The third event was the gradual shift to a drier and colder climate on the QTP, with the development of glaciers accompanying the uplift . At this time, E. godlewskii was constricted into two lineages and dispersed to northern and southwestern China along both flanks of the western Qinling Mountains, which underwent extensive uplift along with the QTP and formed new basins for settlement . The lineages refuged in the montane basins and gorges during the Quaternary glaciation. The species was separated into two lineages by the West Qinling Mountains, which is consistent with the traditional zoogeographical boundary between the Oriental and Palaearctic realms in China . However, the observed phylogeographic pattern of the two lineages is at least partly explained by ecological dissimilarities within southern China (southern lineage) and northern China (northern lineage), especially in habitat preference.
Divergence dating and historical demographic reconstruction
Another important question need to be addressed is whether apparent gene flow corresponding to secondary contact and admixture between the two lineages exits in a continuous time period. However, the results from IMa2 indicated no significant historical gene flow between the southern and northern lineages. Therefore, we simulated the demographic expansion of major lineages of this species based on a null model of isolation without migration (gene flow) between diverging populations with DIY-ABC, and the results indicated that both lineages underwent recent demographic expansions. The BSP indicated that the population experienced an expansion earlier than the LGM, predominantly during the interglacial periods (MISs 2–6) . The northern and southern lineages expanded 0.05 Ma and 0.12 Ma, respectively, during the last interglacial period . The demographic history analysis indicated that the distribution of this bird expanded before the LGM. This time is in contrast to the post-LGM expansion observed in most European bird populations and reflects the different timescales and frequencies of glaciations that occurred in Europe and Asia during the Pleistocene . In fact, a warmer climate in the LIG provided opportunities for each lineage to expand their distribution range. However, the ENM showed different modes of distribution range shift in the two lineages, which may have resulted from different effects caused by the climate of the LGM and the subsequent habitat changes accompanying the arrival of a colder climate in northern and southern regions in China. For the southern lineage, the LGM distribution was contracted compared to the LIG, MIH and CURRENT. This likely relate to the landform heterogeneity in the Southwest Mountains of China [1, 54]. This coherent with BSP result of a mild expansion in the southern lineage. While for the northern lineage, a convergence was occurred and which contracted together into some refugia during the LGM, followed by partial expansion in the MIH after the LGM. This pattern may be related to the relative homogeneity of landforms in northern China, where the bird populations were more likely to be affected by the colder climate in the LGM.
The diverse local-topography and heterogonous climates in China during the Pleistocene can be expected to have varied effects on species-distributions and historical demography. As revealed by Su et al. , the climate during Pleistocene displayed an overall warm and humid trend, with a series of warming and cooling cycles. The vegetation in the Jiuquan Basin included spruces, pine, sagebrush and Solanaceae species approximately 0.04–0.012 Ma. Shi et al.  described a warm stage during the interstage of the last glaciation (0.06–0.03 Ma), when the temperature was approximately 3–4 °C higher than it is today. In the MIS3a stage (0.035 Ma), when the temperature in the Tsaidam Basin was 2 °C higher than today, the alpine and sub-alpine meadows turned into alpine coniferous forests and meadows in Zoige , and the distribution of vegetation expanded towards northern and western China . These findings indicate that the forest developed and shrub meadow habitats arose near forest edges in northwestern China. Therefore, E. godlewskii rapidly expanded along with the suitable habitat in northern China in a short period of time, as indicated by the BSP. However, the glaciers in southwestern China retreated during MIS5 (0.11–0.071 Ma) , and a series of prolonged, mild interglacial periods developed , corresponding to the unique warming stage, MIS5e (0.125 Ma), in the LIG period . Pollen studies have indicated that both alpine meadow and steppe had spread to the eastern part of the plateau during the glaciations . Therefore, it is possible that this species expanded its distribution range during the warmer MIS5 period of the mild interglacial period, which resulted in vegetation similar to the present-day flora of East Asia .
With a detailed phylogeographic structure of Godlewski’s Bunting, E. godlewskii, we provide new insights into the potential mechanisms of genetic divergence and historical dynamics of avian taxa in China. A south-north lineage divergence in E. godlewskii was identified and inferred related to geographic barrier of the West Qinling Mountains. The observed phylogeographic pattern may also have been driven by series geological events like uplift of QTP and subsequent environmental changes in glaciation. In particular, the two lineages showed different models of expansion in response to habitat changes during the Pleistocene glacial-interglacial cycles. Our findings may inform future investigations and potential taxonomic revisions.
Availability of data and materials
All sequences of Godlewski’s Bunting generated during the current study are available in the NCBI GenBank database under accession numbers MG194539-MG194728 and MG228461-MG229050. The data sets supporting the results and conclusions of this article are available in Additional file 1: Tables S3 and S4.
Analysis of molecular variance
Bayesian skyline plots
Cytochrome oxidase I
Effective sample size
Beta-fibrinogen intron 5
Last Glacial Maximum
Million years ago
Markov chain Monte Carlo
Marine isotope stages
Muscle-specific kinase intron
Number of haplotypes
Polymerase chain reaction
Ramos-Onsins and Rozas
Number of segregating sites
Sum of squared deviations
Time of the most recent common ancestor
Lei F, Qu Y, Song G. Species diversification and phylogeographical patterns of birds in response to the uplift of the Qinghai-Tibet plateau and quaternary glaciations. Curr Zool. 2014;60:149–61.
Song G, Zhang R, Per A, Martin I, Cai T, Qu Y, et al. Complete taxon sampling of the avian genus pica (magpies) reveals ancient relictual populations and synchronous late-pleistocene demographic expansion across the northern hemisphere. J Avian Biol. 2017;e01612:1–14.
Zhang R, Song G, Qu Y, Alström P, Ramos R, Xing X, Ericson PGP, Fjeldså J, Wang H, Yang X, et al. Comparative phylogeography of two widespread magpies: importance of habitat preference and breeding behavior on genetic structure in China. Mol Phylogenet Evol. 2012;65:562–72.
Lange R, Durka W, Holzhauer SIJ, Wolters V, DiekÖTter TIM. Differential threshold effects of habitat fragmentation on gene flow in two widespread species of bush crickets. Mol Ecol. 2010;19:4936–48.
Song G, Zhang R, Qu Y, Wang Z, Dong L, Kristin A, Alström P, Ericson PGP, Lambert DM, Fjeldså J, et al. A zoogeographical boundary between the Palaearctic and Sino-Japanese realms documented by consistent north/south phylogeographical divergences in three woodland birds in eastern China. J Biogeogr. 2016;43:2099–112.
Zhao N, Dai C, Wang W, Zhang R, Qu Y, Song G, Chen K, Yang X, Zou F, Lei F. Pleistocene climate changes shaped the divergence and demography of Asian populations of the great tit Parus major: evidence from phylogeographic analysis and ecological niche models. J Avian Biol. 2012;43:297–310.
Roselaar CS, Sluys R, Aliabadian M, Mekenkamp PGM. Geographic patterns in the distribution of Palearctic songbirds. J Ornithol. 2007;148:271–80.
Tan Z. The calculation and simulation of Chinese south–north demarcation based on GIS (in Chinese). Lanzhou: Master Thesis, Lanzhou University; 2011.
Song G, Qu Y, Yin Z, Li S, Liu N, Lei F. Phylogeography of thealcippenmorrisonia (aves: timaliidae): long population history beyond late pleistocene glaciations. BMC Evol Biol. 2009;1:143.
Qu Y, Song G, Gao B, Quan Q, Ericson PGP, Lei F, et al. The influence of geological events on the endemism of east asian birds studied through comparative phylogeography. J Biogeography. 2015;1:179–92.
Päckert M, Sun Y-H, Strutzenberger P, Valchuk O, Tietze DT, Martens J. Phylogenetic relationships of endemic bunting species (Aves, Passeriformes, Emberizidae, Emberiza koslowi) from the eastern Qinghai-Tibet plateau. Vertebrate Zool. 2015;65:135–50.
Byers C, Olsson U, Curson J. Buntings and sparrows : a guide to the buntings and north American sparrows. East Sussex: Pica Press; 1995.
Rose C. Family Emberizidae (buntings and New World sparrows) - species accounts Emberiza. In: Hoyo J, Elliot A, Sargatal J, editors. Handbook of the birds of the world, volume 16, tanagers to New World blackbirds. Barcelona: Lynx Edicions; 2011. p. 508–36.
Fu T, Song Y, Gao W. Fauna Sinica, Aves Vol. 14. Passeriformes: Ploceidae and Fringillidae. Beijing: Science Press; 1998.
Handbook of the Birds of the World, Family Old World Buntings (Emberizidae). Barcelona,1989. http://www.hbw.com/species/godlewskis-bunting-emberiza-godlewskii. Accessed 18 Aug 2018.
Gill F, Donsker D, editors. IOC World Bird List (v8.2); 2018. https://doi.org/10.14344/IOC.ML.8.2.
Cheng T, Qian Y, Tan Y, Zheng B, Guan G, Li G, Min Z, Chen F, Zhao T, Shi D. The avifauna of Qinling mountain. Beijing: Science Press; 1973.
Beaumont MA, Zhang W, Balding DJ. Approximate Bayesian computation in population genetics. Genetics. 2002;162:2025–35.
Nielsen R, Wakeley J. Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics. 2001;158:885–96.
Guisan A, Thuiller W. Predicting species distribution: offering more than simple habitat models. Ecol Lett. 2005;9:993–1009.
Stephens M, Smith NJ, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001;68:978–89.
Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian Phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.
Nylander J. MrModeltest v2. Program distributed by the author. Evolutionary Biology Centre. Uppsala: Uppsala University; 2009.
Drummond AJ, Ho SYW, Phillips MJ, Rambaut A. Relaxed Phylogenetics and dating with confidence. PLoS Biol. 2006;4:e88.
Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:1–8.
Heled J, Drummond AJ. Bayesian inference of species trees from multilocus data. Mol Biol Evol. 2010;27:570–80.
Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu C-H, Xie D, Suchard MA, Rambaut A, Drummond AJ. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10:e1003537.
Bouckaert RR. DensiTree: making sense of sets of phylogenetic trees. Bioinformatics. 2010;26:1372–3.
Huelsenbeck JP, Ronquist F. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics. 2001;17:754–5.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.
Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and windows. Mol Ecol Resour. 2010;10:564–7.
Bandelt HJ, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16:37–48.
Weir JT, Schluter D. Calibrating the avian molecular clock. Mol Ecol. 2008;17:2321–8.
Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.
Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:915–25.
Ramosonsins SE, Rozas J. Statistical properties of new neutrality tests against population growth. Mol Biol Evol. 2002;19:2092–100.
Rogers AR, Harpending H. Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992;9:552–69.
Cornuet JM, Pudlo P, Veyssier J, Dehne-Garcia A, Gautier M, Leblois R, Marin JM, Estoup A. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics. 2014;30:1187–9.
Ren G, Mateo RG, Liu J. Genetic consequences of Quaternary climatic oscillations in the Himalayas Primula tibetica as a case study based on restriction site-associated DNA sequencing. New Phytol. 2017;213:1500–12.
Hey J. Isolation with migration models for more than two populations. Mol Biol Evol. 2010;27:905–20.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2010;25:1965–78.
Feng Y, Cao X, Zhang E, Hu Y, Pan X, Yang J, Jia Q, Li W. Tectonic evolution framework and nature of the west Qinling Orogenic Belt. Northwest Geol. 2003;36:1–10.
Bao X, Liu N, Qu J, Wang X, An B, Wen L, Song S. The phylogenetic position and speciation dynamics of the genus Perdix (Phasianidae, Galliformes). Mol Phylogenet Evol. 2010;56:840–7.
Chen F. Hengduan event: An important tectonic event of the late cenozoic in eastern asia. J Mt Res. 1992;4:195–204.
Zhu H, Chen Y, Pu P, Wang S, Zhuang D. Environments and sedimentation of fault lakes, Yunnan province. Beijing: Science Press; 1989.
Wu G. The modes and mechanism of quaternary fault movement in Lijiang-Dali area, northwestern Yunnan and their influence on environment. Quaternary Sciences. 1992;3:265–76.
Shi X. Late cenozoic uplift of the yulong snow mountain (5596m), se tibetan plateau, caused by erosion and tectonic forcing. Quat Sci. 2008;28:222–31.
Zhang R. Geological events ad mammalian distribution in China. Acta Zool Sin. 2002;48:141–53.
Zhou Y, Lu X, Xu J, Zhang H, Jiang T. China climate change impact report: Yunnan province. Beijing: China Meteorological Press; 2011.
Xiao XY, Shen J, Xiao HF. Paleovegetation and paleoclimate of the Heqing Basin during 2.780–1.802 MaB.P. in Yunnan Province, China (in Chinese). Quat Sci. 2007;27:417–26.
Wu Y, Cui Z, Liu G, Ge D, Yin J, Xu Q, Pang Q. Quaternary geomorphological evolution of the Kunlun pass area and uplift of the Qinghai-Xizang (Tibet) plateau. Geomorphology. 2001;36:203–16.
Wang K, Zhang J, Lei F. Geographical distribution pattern and its temporal and spatial variation of birds breeding in zoogeographical subregion of China. Acta Zootaxon Sin. 2010;35:145–57.
Shi Y, Zheng B, Su Z. Glaciations, glacial-interglacial cycles and environmental changes in the quaternary. In: Shi YF, editor. Glaciers and their environments in China. Beijing: Science Press; 2000. p. 320–55.
Shi Y. Characteristics of late Quaternary monsoonal glaciation on the Tibetan plateau and in East Asia. Quat Int. 2002;97:79–91.
Su J, Wu Y, Li Q, Zhang Y, Wen X. Environmental evolution of the Jiuquan Basin and its relation with the uplift of the Qilian Mountains since the quaternary. Acta Geosci Sin. 2005;26:443.
Shi Y. Uplift of the Qinghai-Xizang (Tibetan) plateau and east asia environmental change during late cenozoic. Acta Geograph Sin. 1999;54:20–8.
Zhao J, Shi Y, Jie W. Comparison between quaternary glaciations in China and the marine oxygen isotope stage (MIS):An improved schema. Acta Geograph Sin. 2011;66:867–84.
Zhou S, Wang X, Wang J, Xu L. A preliminary study on timing of the oldest Pleistocene glaciation in Qinghai–Tibetan plateau. Quat Int. 2006;44:154–5.
Yuan D, Cheng H, Edwards RL, Dykoski CA, Kelly MJ, Zhang M, Qing J, Lin Y, Wang Y, Wu J, et al. Timing, duration, and transitions of the last interglacial Asian monsoon. Science. 2004;304:575–8.
Liu J, Yu G, Chen X. Palaeoclimate simulation of 21 ka for the Tibetan plateau and eastern Asia. Clim Dyn. 2002;19:575–83.
We thank Bo Du, Sen Song, Lixun Zhang, Aiwu Jiang, Wei Zhao, Qingshan Zhao, Mengmeng Guan, Ying Wang, Hongxing Niu, Yinghong Hao and Yuehua Sun for their assistance with field sampling, Guangpeng Ren for assistance with the approximate Bayesian computation analysis in DIY-ABC and Yanhua Qu for assistance with statistics. The authors dedicate this manuscript to their esteemed professor Naifa Liu who departed in February 2019.
This study was supported by grants from the National Natural Science Foundation of China (Grant Nos. 41071031, 31672296 and 31471991). The funders had no role in the design of the study; collection, analysis and interpretation of data; and writing of the manuscript.
Ethics approval and consent to participate
Animal studies, experiments and data collection for this study were conducted following national laws and the guidelines of Lanzhou University’s Institutional Care and Use Committee. All processes in this study were approved by Lanzhou University’s Institutional Care and Use Committee.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1 The sampling localities for Emberiza godlewskii. Table S2 The PCR primers used in this study. Table S3 Descriptions of prior settings for all parameters used in DIY-ABC. Table S4 Genbank Accession numbers for individuals of mitochondrial genes used in this study. Table S5 Genbank Accession numbers for individuals of nuclear genes used in this study. Table S6 Estimations of posterior distributions of parameters revealed by DIY-ABC for the best scenario of demographic history of northern and southern population, respectively. Figure S1 Demographic expansion and divergence models of Emberiza godlewskii. ABC results showing demographic history. Effective population sizes are marked in different colours and times of events (not to scale) are indicated. Figure S2 (a) Schematic representation of ABC modelling of changes in population sizes. Na, ancestral population size; N1, current population size; Ne and Nb, population sizes between Na and N1; t2, old expansion time; tb, bottleneck time; t1, recent expansion time. (b) Posterior probabilities obtained by logistic regression of 1% of the closest simulated datasets for the northern and southern lineages, respectively, from left to right. (DOCX 279 kb)
About this article
- Emberiza godlewskii
- Uplift of the Qinghai-Tibet plateau
- West Qinling Mountains
- Historical demographic dynamics
- Population expansion