Elevational divergence in the great tit complex revealed by major hemoglobin genes.

Gene flow and demographic history can play important roles in the adaptive genetic differentiation of species, which is rarely understood in the high-altitude adaptive evolution of birds. To elucidate genetic divergence of populations in the great tit complex (Parus major, P. minor and P. cinereus) at different elevations, we compared the genetic structure and gene flow in hemoglobin genes with neutral loci. Our results revealed the elevationally divergent structure of α A -globin gene, distinctive from that of the β A -globin gene and neutral loci. We further investigated gene flow patterns among the populations in the central-northern (> 1,000 m a.s.l.), south-eastern (< 1,000 m a.s.l.) regions and the Southwest Mountains (> 2,000 m a.s.l.) in China. The high-altitude (> 1,000 m a.s.l.) diverged α A -globin genetic structure coincided with higher α A -globin gene flow between highland populations, in contrast to restricted neutral gene flow concordant with the phylogeny. The higher α A -globin gene flow suggests the possibility of adaptive evolution during population divergence, contrary to the lower α A -globin gene flow homogenized by neutral loci during population expansion. In concordance with patterns of historical gene flow, genotypic and allelic profiles provide distinctive patterns of fixation in different high-altitude populations. The fixation of alleles at contrasting elevations may primarily due to highland standing variants α A 49Asn/72Asn/108Ala originating from the south-western population. Our findings demonstrate a pattern of genetic divergence with gene flow in major hemoglobin genes depending on population demographic history.

Widespread species in hierarchical population structure can provide opportunities for investigating the potential influence of demographic history in shaping genetic variation arising in response to environmental change. A comprehensive example is a fitness-correlated morphological gene in the adaptive evolution of marine and freshwater three-spine stickleback Gasterosteus aculeatus responding to gene flow and demographic history (Schluter et al. 2004; Colosimo et al. 2005; Schluter and Conte 2009; Raeymaekers et al. 2014). For high-altitude adaptive evolution, genetic variation of oxygen-carrying hemoglobin (Hb) has long been used to track natural selection in native vertebrates under hypoxia stress (Perutz 1983; Monge and Leon-Velarde 1991; Scott and Milsom 2006; Storz 2007; Weber 2007). However, intraspecific genetic divergence in Hb genes may not necessarily cause adaptive changes in Hb in birds (McCracken et al. 2009a; Munoz-Fuentes et al. 2013; Cheviron et al. 2014; Galen et al. 2015; Natarajan et al. 2015). A hypothesis proposed that highland adaptation after colonization from lowlands following divergence based on the Hb gene flow from lowland to highland in Andean waterbirds Anatidae (McCracken et al. 2009b; Wilson et al. 2012). The hypothesis was subsequently supported by highland genotypes derived from ancestral lowland genotypes in functional experiments (Natarajan et al. 2015). For avian species at broad elevations, whether the genetic structure and gene flow in Hb genes can reveal intraspecific adaptive differentiation in elevations is still an open question.

To investigate the role of gene flow and demographic history in elevational genetic divergence, we performed population genetic analyses on major Hb genes and unlinked neutral loci in a widespread species complex, the great tit. The HbA tetramer, which is the major Hb isoform expressed in adult birds, comprises 2 α-chains and 2 β-chains encoded by α- and β-globin gene, respectively (Mairbaurl and Weber 2012; Opazo et al. 2015). It has been revealed that the adaptive genetic variants mostly presented in the β-globin gene in high-altitude Andean birds (Projecto-Garcia et al. 2013; Galen et al. 2015; Natarajan et al. 2015). In contrast, rare causative variants were found in the α-globin gene of cinnamon teals Anas cyanoptera and nightjars Hydropsalis ( Wilson et al. 2012; Natarajan et al. 2015; Kumar et al. 2017). The great tit complex, mainly comprising Parus major, P. minor, and P. cinereus, is known as a superspecies encircling the Qinghai–Tibet Plateau (QTP) and Central Asian desert (Kvist et al. 2003; Päckert et al. 2005; Kvist et al. 2007; Gill and Donsker 2016). This superspecies is widely distributed from sea level to above 4,000 m a.s.l. in Eurasia, with the highest habitats located in the eastern Himalaya (Gosler et al. 2017). A monophyletic clade of P. minor in eastern Himalaya (P. m. tibetanus and P. m. subtibetanus, clade B) has been further identified to compose the 5-clade mtDNA phylogeny of the great tit complex (Zhao et al. 2012). The demographic history of P. minor revealed historically stable growth in the eastern Himalaya populations before the last glacial maximum (LGM) and rapid population expansion in the East Asia during the LGM (Zhao et al. 2012; Qu et al. 2015).

In the present study, we conducted population genetic analyses for exploring the role of gene flow and demographic history in Hb genetic differentiation corresponding to changes in elevation in great tit populations. This study mainly focused on 1) comparing the difference in genetic structure of Hb genes with neutral loci; 2) examining gene flow among local populations at hierarchical elevations of P. minor; 3) identifying genotypes and profiles along elevational gradients and 4) elevational variation of allelic frequencies in Hb genes across parallel transects.

Materials and Methods

Sampling, sequencing and gametic phasing

To determine the genetic structure of all populations, we employed a total of 159 samples from all 5 clades of the 3 species (P. minor, P. cinereus and P. major) within the great tit complex (Supplementary Table S1). All the collections were carried out with permission from the State Forestry Administration of the People’s Republic of China and in accordance with the National Wildlife Conservation Law of China. To explore local population differentiation in response to changes in elevation, we analyzed a sub-dataset of P. minor populations. According to the genetic structure partitioned by elevation (Figure 1), we sampled highland (1,000–4,282 m a.s.l.) and lowland individuals (7–210 m a.s.l.) in 3 local populations of P. minor (N = 66) (Supplementary Table S1). The 3 populations were located in the Southwest Mountain region (highland Clade B, HB, N = 23), the alpine region of central China (highland Clade A, HA, N = 20), and the plain of southeast China (lowland Clade A, LA, N = 23). The HB and HA populations were adjacent to the southeast and the east of the QTP, respectively (Figure 2A).

Figure 1

Distribution of samples and unrooted haplotype networks in samples of the Great Tit complex. (A) Samples collected in species maps. Colorful shadows represent the distribution of each clade in neutral phylogeny, corresponding to the nodes in mtDNA dendrogram (Zhao et al. 2012). (B) In networks, each circle represents for a haplotype and circle size for the individual numbers with a shared haplotype. Pied colors represent for different phylogenetic clades in areas proportional to haplotype frequencies. An elevational differentiation presents in α-globin gene in the samples above 1,000 m a.s.l of clades A, B and E (bold nodes in dash-line rectangles, Supplementary Table S1). The presence of non-synonymous substitution sites are labeled in Hb networks.

Figure 2

Gene flow diagrams and networks in P. minor populations. (A) The collection localities of 3 local populations of clade A and B in dash circles (Figure 1). The sampling localities of HB, HA, and LA populations are in black, gray and white. The direction of gene flow between elevational differentiated populations is in red and white arrows for α-globin gene and neutral loci, respectively. The strength of gene flow is classified into relatively strong as solid and weak as dash lines compared in α-globin gene or neutral loci. (B-C) The genetic structures of major Hb genes and neutral reference loci (mtDNA and nucleotide introns) were illustrated in networks. Black, gray and white colors are the same individual of the 3 populations in panel A. The presence of nonsynonymous substitution sites are labeled in Hb networks (B).

We amplified and sequenced α- and β-globin genes and 6 unlinked reference loci (2 mtDNA and 4 autosomal introns) for all samples. The reference loci are widely used neutral markers in general phylogenetic studies, including the mitochondrial control region (CR, 539 bp), NADH dehydrogenase subunit 2 (ND2, 737 bp), beta-fibrinogen intron 5 (Fib, 281 bp), transforming growth factor beta 2 intron 5 (TGFB2, 334 bp), myoglobin intron 2 (Myo, 556 bp) and ornithine decarboxylase intron 7 (ODC, 284 bp). We designed subunit-specific primers for Hb genes (660 bp of α-globin gene and 1,274 bp of β-globin gene). After DNA extraction from venous blood or pectoral muscle with a DNeasy Tissue kit (QIAGEN, Hilden, Germany), we amplified the 8 genes with KD-Plus DNA polymerase (TransGen Biotech, Beijing, China). The general PCR thermal profile was 5 min at 94°C, 35 cycles of 40 s at 94 °C, 40 s at 52–64 °C, 1 min at 72 °C, and 10 min at 72 °C (with varied annealing temperatures, as shown in Supplementary Table S2). For each gene, we sequenced the PCR products on the ABI3730 sequencing system (BGI, Beijing, China).

For heterozygous autosomal genes, we inferred the gametic phases using PHASE v2.1 (Stephens et al. 2001). The PHASE algorithm was run 5 times using random seeds to obtain stable results. For the α-globin gene (660 bp) and 4 autosomal introns mentioned above, we inferred the gametic phases of the complete sequences. For the β-globin gene (1,274 bp), we inferred the gametic phases for the fragments spanning from exon 1 to exon 2 and partial intron 2 (453 bp), which had no recombination detected in DnaSP v5.10.1 (Librado and Rozas 2009). All substitutions that occurred in the exon 3 (126 bp) of the β-globin gene were synonymous. The sequences with high posterior probabilities (>95%) and the consensus sequence of the cloned samples unresolved in PHASE proceeded into downstream analyses. All sequences are accessible under accession numbers MF321902-MF322504 in GenBank, NCBI, NIH (Supplementary Table S2).

Comparison of genealogical structures between hb genes and neutral loci

To understand the genealogical discordance between Hb genes and neutral loci, we reconstructed unrooted networks of haplotypes for each gene in local and all samples by using the median-joining algorithm in NETWORK 4.6.1 (Bandelt et al. 1999). Networks based on the coding sequence (CDS) of the Hb genes were also reconstructed to exclude the effects of introns in Hb genes. For qualitative analysis of the genetic structures of local populations, we calculated the posterior probabilities of assigned individuals into highland and lowland populations of P. minor in STRUCTURE 2.3.4 (Pritchard et al. 2000). We performed STRUCTURE analyses on 4 different sub-datasets: all nuclear introns, the introns plus α-globin gene, the introns plus β-globin gene, and the introns plus α- and β-globin genes. The analyses were performed for 1 or 2 population models (K = 1 or 2) without prior sampling and conducted with admixture model and independent allele frequencies (α = 1, λ = 1). The algorithm was run for 100,000 iterations with a burn-in of 25,000 iterations.

We calculated polymorphism and diversity of all loci in each population and in all 3 populations of P. minor combined in DnaSP v5.10.1 (Pritchard et al. 2000). Tajima’s D (Tajima 1989) and Fu’s Fs (Fu and Li 1993; Fu 1997) statistics were computed to test deviation from neutrality. The significant deviation from neutrality is an outlier of 95% confidence interval in a null model of neutral expectation in coalescent simulation based on segregated sites and recombination rate per gene. We calculated pairwise FST values for each of 8 loci and the CDS of Hb genes in ARLEQUIN 3.5.

Genetic diversity and gene flow in local populations

We used a coalescent-based “isolation-with-migration” model implemented in IMa (Hey and Nielsen 2007) to estimate inter-population gene flow in the sampled P. minor populations. We focused on the gene flow strength in the same dataset (α-globin or reference loci) and asymmetric gene flow ratio scaled by migration rates (m* = 2 Nm) i.e., locus-specific effective population size rather than mutation rates. We thus conducted analyses for the α-globin gene and joint reference loci (2 mtDNA and 4 autosomal introns, Supplementary Table S3). We detected intra-locus recombination with the 4-gamete test and retained the longest non-recombining block of sequences in each locus under the IMa assumption of the absence of recombination within a locus (Hey and Nielsen 2004). According to the biasing effects of substitution models, we fitted the Hasegawa–Kishino–Yano (HKY) and Infinite Sites (IS) mutation models for mtDNA and autosomal DNA regardless of calculated models (Strasburg and Rieseberg 2010). Final parameter estimations were scaled to effective population size multiplied by mutation rate (θ = 4 Nμ) with estimated migration rates (m = m/μ) transformed to 2 Nm (M = θ × m/2) (Hey and Nielsen 2004). Given inheritance scalars of 0.25 (mtDNA) and 1.00 (autosomal DNA), we initially performed IMa runs with flat priors for each parameter. On the basis of preliminary runs, we tuned upper bounds for priors that encompassed the overall posterior distributions. We finally ran the analyses in 6 Markov-coupled chains with geometric heating (g1 = 0.8; g2 = 0.9), discarding at least 20,000 samples as “burn-in”. Each chain was continued until the effective sample size (ESS) was >100. To ensure the convergence of parameter estimates, 2 replicated runs in each analysis were conducted by using different random seeds.

Genotypic profiles of hb genes and elevational clines of hb alleles

To characterize possible genotypes differentiated between highland and lowland populations, we estimated genotypic frequencies in Hardy–Weinberg equilibrium (HWE) for haplotypes and genotypes of the 2 Hb genes in ARLEQUIN 3.5 (Excoffier and Lischer 2010). To explore the threshold of elevational genetic variation, we conducted a maximum likelihood estimation of Hb allele frequencies by assigning individuals of P. minor populations into 2 parallel elevational transects: T1 in the south and T2 in the north (Supplementary Figure S2). We used the cline shapes to acquire the profile of highland replacements and a modified Markov Chain Monte Carlo (MCMC) algorithm to determine support limits in ClineFit 2.0b (Porter et al. 1997). The 4-parameter model for the cline center (c), cline width (w), minimum allele frequency (P) and maximum allele frequency (P) at the ends of the cline shapes was estimated. We set a burn-in of 800 parameter tries per annealing step and 1000 tries for the second run that are combined in sampling; subsequently, 2000 replicates were saved in the sampling parameter space, with 20 replicates between each save. To assess the coincidence of the cline center and the concordance of the width, we combined the sampled parameter space of replicated runs according to the maximum 2 LnL limits and subsequently generated support plots (not shown) and sigmoidal cline graphics in Mathematica v10. The FST for each non-synonymous substitution between highland and lowland populations was calculated using locus-by-locus AMOVA in ARLEQUIN 3.5 (Excoffier and Lischer 2010).

Results

The population genetic structure along elevations

For both the great tit superspecies and P. minor, significant genetic divergence was found in Hb genes, especially for the α-globin gene (Figures 1 and 2). The divergence of the α-globin gene partitioned the superspecies into 2 elevational clades: the monophyletic clade of high-altitude individuals from the phylogeographic clade A, B, C, and E, and the other clade including the clade D and low-altitude individuals from the clade A (Figure 1). Similarly, the network of the α-globin gene in P. minor populations presented the monophyletic highland clade including HA and HB, which was diverged from the lowland clade-A population (LA) (Figure 2). The β-globin gene, however, did not exhibit any clear genetic structure. The single β21Ala/Thr substitution presented in the southwest mountains was not fixed in frequency (Figure 2 and 4). Therefore, in further analyses of genetic variance, we mainly used the local population dataset to explore a potential genetic response to the altitudinal gradient in P. minor.

Figure 4

Genotype frequencies in amino acid alleles of major Hb genes. FST represents for differentiation of each allele between highland and lowland populations, *P < 0.05. The striking heterozygous highland and lowland genotypes involve α22Asp/Gly in LA population and β21Ala/Thr in HB population, in contrast to low-level admixed in α35/57. The more homozygous α49Asn/Ser-72Asn/His-108Ala/Val is an admixed genotype only in HA population at lower frequency while completely fixed in HB vs. LA population.

The genetic structure in differences of elevation, which could be largely due to differences in the α-globin gene, was further confirmed by STRUCTURE analysis on local populations (Figure 3). The different assignment to lowland or highland (>1,100 m a.s.l.) populations was most dramatic in average posterior probability of the α-globin (6.1% versus 88.2%, ±8.0% and 18.7% standard deviation [SD]) and less so of joint α- and β-globin genes (8.6% versus 82.6%, ±10.9% and 22.0% SD). In contrast, the assignments showed much less divergence in response to elevation either for the single β-globin (43.6% and 55.3%) or for nuclear introns (55.7% and 44.8%).

Figure 3

Posterior probability of individual assignment to highland and lowland populations in STRUCTURE 2.3. Black, gray, and white represent individuals from population HB, HA and LA (Fig. 2). (A) 4 introns of reference loci, (B) 4 reference introns and α-globin gene, (C) 4 reference introns and β-globin gene, (D) 4 reference introns and 2 hemoglobin genes. The higher divergence occurs in the dataset including α-globin gene but lower with β-globin gene add-in.

Gene flow among elevationally differentiated populations of P. Minor

The coalescent estimations of IMa showed unbalanced downslope gene flow (from highland to lowland) in the joint reference loci (assumed as neutral gene flow) among local P. minor populations (Table 1; Figure 2A). The highest gene flow was found in the downslope (HA to LA) neutral loci and in 16.1, 7.3 and 18.5 times higher level than the upslope (LA to HA) and the gene flow between 2 clades (HB-LA and HB-HA) (Table 1). For the α-globin gene, in contrast, downslope gene flow between clades (HB-HA) was 23.9 times higher than upslope gene flow within the clade A (LA-HA). The lower gene flow of neutral loci is between 2 clades (HB versus HA/LA), while the lower α-globin gene flow is between populations differed by elevations (LA versus HB/HA). Relative to the restrained neutral gene flow, the much higher gene flow in α-globin from HB into HA is in accordance with the elevational structure of local populations.

Table 1

Coalescent estimation of population genetic parameters among the P. minor populations using IMa

Loci	Populations	θH	θL	mH	mL	mH*	mL*	mL*/mH*
	(High versus lowland)
All neutralloci	HA versus LA	0.43	1.68	2.50	10.36	0.54	8.69	16.09
		(0.18–1.18)	(0.78–4.82)	(0.68–41.37)	(3.24–45.80)
	HB versus LA	1.14	1.09	0.00	2.18	0.00	1.19	NA
		(0.66–2.03)	(0.64–2.00)	(0.02–4.64)	(0.25–5.51)
	HB versus HA	1.11	0.53	0.00	1.76	0.00	0.47	NA
		(0.67–1.93)	(0.28–1.14)	(0.01–2.81)	(0.42–5.62)
αA–globingene	HA versus LA	0.43	0.41	1.03	0.01	0.22	0.00	NA
		(0.16–2.17)	(0.13–1.93)	(0.27–8.06)	(0.02–5.60)
	HB versus LA	0.49	0.60	0.00	0.00	0.00	0.00	NA
		(0.13–2.07)	(0.23–2.26)	(0.01–2.28)	(0.01–2.10)
	HB versus HA	0.26	1.20	0.01	8.73	0.00	5.25	NA
		(0.10–1.57)	(0.52–31.37)	(0.12–19.41)	(2.73–51.57)

The neutral parameter θ = 4 Nμ, N is the effective population size of each population, μ is the mutation rate. mH, the migration rate from lowland to highland (upslope); mL, migration rate from highland to lowland (downslope). The upslope gene flow between HB and HA is from the relative lowland east QTP (HA) to the highland west QTP (HB), the opposite is downslope from HB to HA (Supplementary Figure S1 and see Appendix S3). The 95% highest posterior density (95% HPD) are shown in parentheses. mH* and mL*, were scaled m by effective population size by m × θ/2 of upslope and downslope gene flow. The mL*/mH* ratio is an index of the asymmetric gene flow. NA is an ineffective ratio of zero m values.

Genetic divergence and genotypic profiles in major Hb genes of P. Minor

We identified nucleotide polymorphisms of all genes in the P. minor populations (Supplementary Table S3). The haplotype diversity was the lowest in ND2 and β-globin gene (H< 0.5) of HA and LA populations compared with other loci (average 0.726 ± 0.21 SD). The nucleotide diversity of both Hb genes and ND2 of clade A exhibited the lowest levels (π < 0.002, average 0.0052 ± 0.0050 SD). Significant negative Tajima’s D and Fu’s F values were present in different populations of β-globin and reference loci (Supplementary Table S3). In α-globin gene of all populations, insignificant negative values were predominant, indicating lower than expected neutrality. Significant negative values predominated among nuclear loci of the LA population and in β-globin for both HA and LA populations, but for mtDNA in the HB population. In P. minor populations of clade A, significantly greater divergence was indicated by fixation indices (FST) in the α-globin gene in contrast to other loci (FST* > 0.5 versus < 0.1, Table 2). The significant inter-clade divergence of mtDNA was consistent with genetic structure. The pairwise FST values of the α-globin gene showed significant differentiation among populations at different elevations (Table 2).

Table 2

Pairwise FST values for nucleotide haplotypes of the P. minor populations

Pairwise	αA	αA-cds	βA	βA-cds	CR	ND2	TGFB2	Fib	Myo	ODC
HA versus HB	0.15	0.20	0.07	0.06	0.81	0.95	0.15	0.26	0.02	0.31
HA versus LA	0.60	0.67	0.00	0.00	0.10	0.00	0.02	0.07	0.05	0.00
HB versus LA	0.76	0.87	0.10	0.13	0.69	0.95	0.10	0.15	0.10	0.29

Significant values (P < 0.05) are shown in bold. Sample sizes were shown in Figure 1.

All the Hb genotypes were profiled in HWE analyses (Supplementary Table S4). The genotype of each Hb gene exhibited differences in frequency between highland and lowland populations. The 6-allele genotype in the α-globin gene was completely fixed in all individuals in the HB population (the same genotype in 46 chromosomes). Accordingly, both the observed genotypic heterozygosity of the α-globin gene in HB was lower than that in the other 2 populations (Ho/He = 0.753, 0.924 and 0.861 in HB, HA, and LA). The β-globin genotype of most individuals was homozygous, and a unique substitution of β21Ala/Thr occurred in HB populations at a low frequency (23.5%). From the allele-specific analysis, α57Gly/Ala (FST = 0.004) and β21Ala/Thr (FST* = 0.168, *P < 0.05) exhibited the lowest differentiation between highland and lowland populations (Figure 4). The α-globin genotypes of 2 distant alleles were fixed according to elevation (FST* > 0.9): α35Ala/72Asn was predominant in highland populations (HA and HB). Moreover, α49Asn/Ser and α108Ala/Val also exhibited high levels of elevational differentiation (FST* = 0.845).

Elevational variation in allele frequencies across parallel transects

Specific α-globin allelic profiles in local populations of P. minor were further determined by cline-fitting analyses along parallel elevational gradients (Figure 5 and Table 3). The highland alleles shifted along altitudinal gradients coincide at a threshold of 1,100 m a.s.l. (±150 m SD) that is consistent with the genetic structure of superspecies sorted by elevation. The center fittings of cline shapes were similar across parallel transects at or above 1,000 m. The average estimated cline center of T1 curves was 1028.4 (817–1,145) m a.s.l., and the cline width was 471 (59–958) m (Table 3). The fully consistent cline centers (1,060, 962–1,152 m a.s.l.) and widths (446, 14–813 m) occurred in α49/72/108 of the T1 transects, whereas the approximately coincident centers of α22/35/108 were in T2. The average center of T2 was slightly higher (1,131.3 m a.s.l.) and wider (707 m) than that of the T1 transect. In the more elevationally fixed alleles (α49/72/108 of T1), higher differentiation (FST = 1.00, P < 0.05) was apparent across elevational transects. Accordingly, the α49-72-108 across T1 exhibited the highest differentiation (FST = 1.00), and α22 across T1 had the lowest differentiation (FST = 0.36; Table 3). According to the genetic structure, the unique β-globin replacement without an elevation pattern could not be used as the marker to reveal the elevation pattern in local populations (Figure 4 and Supplementary Figure S3).

Figure 5

Elevational cline fitting of α-globin allele frequencies in parallel transects of P. minor. Highland allele frequencies at a given elevation are in red points, the support limit of cline center is in gray shade. The shift at ca 1, 100 m a.s.l. appears to be a threshold of most high-altitude Hb alleles in the P. minor populations. The coincident cline varied elevational shifts occurred in α49/72/108 in T2 transect and α22/35 in T1 transect. The genotypic frequency of α49Asn/Ser-72Asn/His-108Ala/Val in T1 transect and α72Asn/His in both transects are fixed in high/low elevations outside the cline shift regions.

Table 3

Maximum likelihood of elevational clines of αA-globin allele frequencies across parallel transects in P. minor

Alleles	Transect	Centre/m	Width/m	PL	PR	FST(*P < 0.05)
αA22	T1	1145	59	0.52	1.00	0.36*
		(1066–1262)	(0–606)	(0.41–0.63)	(0.96–1.00)
	T2	1106	719	0.25	1.00	0.68*
		(440–1415)	(4–2275)	(0.00–0.44)	(0.92–1.00)
αA35	T1	817	958	0.00	1.00	0.97*
		(666–1048)	(385–1387)	(0.00–0.15)	(0.95–1.00)
	T2	1027	764	0.08	1.00	0.86*
		(632–1269)	(17–1713)	(0.00–0.24)	(0.92–1.00)
αA49	T1	1060	446	0.00	1.00	1.00*
		(962–1152)	(14–813)	(0.00–0.08)	(0.96–1.00)
	T2	861	31	0.00	0.67	0.60*
		(181–1203)	(0–1539)	(0.00–0.09	(0.52–0.82)
αA72	T1	1060	446	0.00	1.00	1.00*
		(962–1152)	(14–813)	(0.00–0.08)	(0.96–1.00)
	T2	1261	638	0.00	1.00	0.84*
		(1081–1449)	(329–1207)	(0.00–0.08)	(0.93–1.00)
αA108	T1	1060	446	0.00	1.00	1.00*
		(962–1152)	(14–813)	(0.00–0.08)	(0.96–1.00)
	T2	1047	56	0.00	0.68	0.60*
		(846–1339)	(0.25–1272)	(0.00–0.07)	(0.53–0.83)

Notes: The parenthesized ranges are parameter support limits up to lnLmax–lnL < 2 likelihood.

The Hb allelic variations in the great tit complex samples comprised the same 6 α-globin substitutions and 1 β-globin substitution as in the P. minor populations. In accordance with local populations, the allele frequencies of α49/72/108 exhibited relatively greater fixed elevational variation and α72 exhibited complete elevational fixation (Supplementary Figure S3). And frequencies of the same allele in highland populations exhibited different results in parallel transects that α49/108 was flexible in T2 (Figure 5 and Supplementary Figure S3). In contrast, the α22/35/57 lacked differentiation along the altitudinal gradient for samples from all clades (Figure 1). The β21Thr occurred at a lower frequency from lowland to highland populations except for fixation at the peak elevation (< 36% at lower than 3,000 m a.s.l.).

Discussion

Elevational differentiation of genetic structure in αA-globin gene

We found an elevationally differentiated structure in the α-globin gene, which is distinctive from the mtDNA phylogeny and the unstructured networks in nuclear loci. This topological mismatch between functional gene and neutral genes might be induced by incomplete lineage sorting or introgressive hybridization between sister lineages, as in sister Anatidae species (Natarajan et al. 2015). The elevationally differentiated genetic structure of the α-globin gene is also distinctive from the unstructured β-globin. This result differs from most Andean birds that elevational differentiated β-globin gene is contributed to the conspecific adaptive divergence (McCracken et al. 2009b; Bulgarella et al. 2012; Galen et al. 2015; Natarajan et al. 2015; but see Wilson et al. 2012). The high- and low-altitude divergence in α-globin gene only occurred in cinnamon teals as similar manner as β-globin gene divergence in other waterfowl (McCracken et al. 2009b; Bulgarella et al. 2012; Wilson et al. 2012). The adaptive divergence in Hb genes of Andean waterbirds consistently contributed to the dispersal of ancestral genotypes in lowland populations by colonization into the highlands. In conclusion, the discrepancy in genetic structures between the 2 major Hb genes might be resulted from adaptive evolution in the α-globin gene of great tits.

Gene flow patterns among local populations of P. minor

From coalescent gene flow estimation of gene flow in P. minor, we found higher than expected α-globin gene flow between clades with restricted neutral gene flow. In contrast to background gene flow revealed by neutral loci, the inter-clade downslope gene flow of α-globin gene is nearly absent between elevation-diverged populations within the clade. Similar patterns exhibited in the restricted β-globin gene flow in the adaptive elevational divergence with highland-derived replacements in Andean waterbirds (McCracken et al. 2009b; Wilson et al. 2012; Natarajan et al. 2015). The α-globin gene flow from subtropical highlands (HB) into temperate lowlands (HA) of P. minor is in a direction similar to the β-globin gene flow from tropical highlands into temperate lowland of south-American waterfowl (McCracken et al. 2009b; Munoz-Fuentes et al. 2013). The restricted α-globin gene flow in the same clade of P. minor may also have been homogenized by high gene flow of neutral loci contingent on population history, as in the pattern shown for major Hb genes in Andean waterfowl (McCracken et al. 2009b). Our results suggest the role of α-globin gene flow in driving genetic divergence along elevational gradients in population history of a widespread species.

Given the demographic history of the great tit (Zhao et al. 2012; Qu et al. 2015), we propose a plausible route for historical population expansion. The α-globin genetic variants may have existed in the high-altitude eastern Himalayan population and dispersed to the lower habitats in central China earlier before the inter-clade divergence in neutral loci. The strong downslope neutral gene flow, in contrast, implies the historical rapid population expansion was from the central highland into the southeast lowland population after inter-clade divergence. Thus, the restricted α-globin gene flow within the same clade may be a result of rapid population expansion in the lowland population, offsetting the possible countervailing selection on high-altitude genetic variants after inter-clade divergence.

Elevational variation of alleles depends on demographic history

The genotypic profile of major Hb genes in HWE coincides with the genetic structure of elevational segregation in the α-globin gene. However, the haplotypic Ho/He ratios of the α-globin gene in the HB and LA populations deviated significantly from equilibrium. This result suggests subunit-specific allelic profiles. From elevational cline fitting, we identified the high/lowland fixed α49/72/108 genotypes that coincide with possible standing variants α49Asn/72Asn/108Ala, originating from an ancestral population in the eastern Himalaya (Tietze and Borthakur 2012; Qu et al. 2015). From the ancestral states reconstruction of Hb genes in Paridae, the highland α-globin haplotype was inferred to be the standing genetic variant in P. minor (Zhu et al., unpublished data). In contrast, the lowland fixed β21Thr is more likely a result of newly derived mutations in highland eastern Himalaya populations and fixed after the rapid expansion. These results implied that standing genetic variation of the α-globin gene might have been fixed in clade B before the inter-clade divergence. The heterozygous α49/72/108 genotypes in clade A, in contrast, may be a mixed result of countervailing selection on standing high-altitude variants and homogenizing neutral gene flow, as the case in Andean waterfowl (McCracken et al. 2009b; Wilson et al. 2012; Natarajan et al. 2015;). Therefore, we suggest that the α-globin gene is likely under historical adaptive selection rather than the β-globin gene in P. minor and the great tit complex.

It is known that adaptive divergence in a wild species might not be determined by phenotypic changes of Hb , whereas the fixation of genetic variants of functional genes eventually depends on demographic history (Olson-Manning et al. 2012; Raeymaekers et al. 2014; Ferchaud and Hansen 2016). The fixation of standing genetic variants in highland populations (HB or HA) was inferred for most alleles (8/10, except for α49/108 in T2) (Zhu et al., unpublished data) that might be homogenized by neutral gene flow. Similarly, much narrower shifts but unfixed highland α49Asn/108Ala alleles in HA and wider cline shifts (α22/35/72 of T2) may have been caused by neutral gene flow in clade A. In agreement with these results, the α22/35 were absent from altitudinal clines in the full sampling, whereas only the completely diverged α72Asn/His suggested that an ancestral standing polymorphism was fixed in all highland populations of the great tit complex. Therefore, although the same elevational differentiation was detected in the genetic structure of the superspecies, specific allele profiles of in P. minor populations were affected by local homogenizing gene flow. In a word, differed from the case of standing Hb variants originated in low-altitude Andean waterfowl (Natarajan et al. 2015), the recurrently replaced alleles in deeply diverged lineages of Eastern Himalaya and Central Asia can be interpreted as the fixation of standing high-altitude alleles in East-Himalaya.

In summary, our study highlights the necessity of integrating gene flow and demographic history into exploring potential adaptive evolution in the genetic divergence of widespread species. This preliminary investigation on genetic variants in a widespread passerine superspecies provides an ideal model to further explore the adaptive evolution of Hb genes in wild birds.

Data archiving

All DNA sequences are accessible under GenBank accession numbers (MF321902–MF322504), and sample information for all individuals is found in Table S1 of the Supplementary file.

Click here for additional data file.