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

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 b -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.) A A 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 expan- sion. In concordance with patterns of historical gene flow, genotypic and allelic profiles provide dis- tinctive patterns of fixation in different high-altitude populations. The fixation of alleles at contrast- ing 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. Key words: elevational divergence, gene flow, genetic structure, hemoglobin gene. Widespread species in hierarchical population structure can provide natural selection in native vertebrates under hypoxia stress (Perutz opportunities for investigating the potential influence of demo- 1983; Monge and Leon-Velarde 1991; Scott and Milsom 2006; graphic history in shaping genetic variation arising in response to Storz 2007; Weber 2007). However, intraspecific genetic divergence environmental change. A comprehensive example is a fitness- in Hb genes may not necessarily cause adaptive changes in Hb in correlated morphological gene in the adaptive evolution of marine birds (McCracken et al. 2009a; Munoz-Fuentes et al. 2013; and freshwater three-spine stickleback Gasterosteus aculeatus re- Cheviron et al. 2014; Galen et al. 2015; Natarajan et al. 2015). A sponding to gene flow and demographic history (Schluter et al. hypothesis proposed that highland adaptation after colonization 2004; Colosimo et al. 2005; Schluter and Conte 2009; Raeymaekers from lowlands following divergence based on the Hb gene flow et al. 2014). For high-altitude adaptive evolution, genetic variation from lowland to highland in Andean waterbirds Anatidae of oxygen-carrying hemoglobin (Hb) has long been used to track (McCracken et al. 2009b; Wilson et al. 2012). The hypothesis was V C The Author (2017). Published by Oxford University Press. 1 This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 2 Current Zoology, 2017, Vol. 0, No. 0 A A subsequently supported by highland genotypes derived from ances- We amplified and sequenced a - and b -globin genes and 6 un- tral lowland genotypes in functional experiments (Natarajan et al. linked reference loci (2 mtDNA and 4 autosomal introns) for all 2015). For avian species at broad elevations, whether the genetic samples. The reference loci are widely used neutral markers in gen- structure and gene flow in Hb genes can reveal intraspecific adaptive eral phylogenetic studies, including the mitochondrial control region differentiation in elevations is still an open question. (CR, 539 bp), NADH dehydrogenase subunit 2 (ND2, 737 bp), To investigate the role of gene flow and demographic history in beta-fibrinogen intron 5 (Fib, 281 bp), transforming growth factor elevational genetic divergence, we performed population genetic ana- beta 2 intron 5 (TGFB2, 334 bp), myoglobin intron 2 (Myo, 556 bp) lyses on major Hb genes and unlinked neutral loci in a widespread and ornithine decarboxylase intron 7 (ODC, 284 bp). We designed species complex, the great tit. The HbA tetramer, which is the major subunit-specific primers for Hb genes (660 bp of a -globin gene and Hb isoform expressed in adult birds, comprises 2 a-chains and 2 1,274 bp of b -globin gene). After DNA extraction from venous blood A A b-chains encoded by a -and b -globin gene, respectively (Mairbaurl or pectoral muscle with a DNeasy Tissue kit (QIAGEN, Hilden, and Weber 2012; Opazo et al. 2015). It has been revealed that the Germany), we amplified the 8 genes with KD-Plus DNA polymerase adaptive genetic variants mostly presented in the b -globin gene in (TransGen Biotech, Beijing, China). The general PCR thermal profile high-altitude Andean birds (Projecto-Garcia et al. 2013; Galen et al. was 5 min at 94 C, 35 cycles of 40 s at 94 C, 40 s at 52–64 C, 1 min at 2015; Natarajan et al. 2015). In contrast, rare causative variants were 72 C, and 10 min at 72 C (with varied annealing temperatures, as found in the a -globin gene of cinnamon teals Anas cyanoptera and shown in Supplementary Table S2). For each gene, we sequenced the nightjars Hydropsalis ( Wilson et al. 2012; Natarajan et al. 2015; PCR products on the ABI3730 sequencing system (BGI, Beijing, China). Kumar et al. 2017). The great tit complex, mainly comprising Parus For heterozygous autosomal genes, we inferred the gametic major, P. minor,and P. cinereus, is known as a superspecies encircling phases using PHASE v2.1 (Stephens et al. 2001). The PHASE algo- the Qinghai–Tibet Plateau (QTP) and Central Asian desert (Kvist rithm was run 5 times using random seeds to obtain stable results. et al. 2003; P€ ackert et al. 2005; Kvist et al. 2007; Gill and Donsker For the a -globin gene (660 bp) and 4 autosomal introns mentioned 2016). This superspecies is widely distributed from sea level to above above, we inferred the gametic phases of the complete sequences. 4,000 m a.s.l. in Eurasia, with the highest habitats located in the east- For the b -globin gene (1,274 bp), we inferred the gametic phases ern Himalaya (Gosler et al. 2017). A monophyletic clade of P. minor for the fragments spanning from exon 1 to exon 2 and partial intron in eastern Himalaya (P. m. tibetanus and P. m. subtibetanus, clade B) 2 (453 bp), which had no recombination detected in DnaSP v5.10.1 has been further identified to compose the 5-clade mtDNA phylogeny (Librado and Rozas 2009). All substitutions that occurred in the of the great tit complex (Zhao et al. 2012). The demographic history exon 3 (126 bp) of the b -globin gene were synonymous. The se- of P. minor revealed historically stable growth in the eastern quences with high posterior probabilities (>95%) and the consensus Himalaya populations before the last glacial maximum (LGM) and sequence of the cloned samples unresolved in PHASE proceeded into rapid population expansion in the East Asia during the LGM (Zhao downstream analyses. All sequences are accessible under accession et al. 2012; Qu et al. 2015). numbers MF321902-MF322504 in GenBank, NCBI, NIH In the present study, we conducted population genetic analyses (Supplementary Table S2). for exploring the role of gene flow and demographic history in Hb genetic differentiation corresponding to changes in elevation in great Comparison of genealogical structures between hb tit populations. This study mainly focused on 1) comparing the dif- genes and neutral loci ference in genetic structure of Hb genes with neutral loci; 2) examin- To understand the genealogical discordance between Hb genes and ing gene flow among local populations at hierarchical elevations of neutral loci, we reconstructed unrooted networks of haplotypes for P. minor; 3) identifying genotypes and profiles along elevational gra- each gene in local and all samples by using the median-joining algo- dients and 4) elevational variation of allelic frequencies in Hb genes rithm in NETWORK 4.6.1 (Bandelt et al. 1999). Networks based across parallel transects. on the coding sequence (CDS) of the Hb genes were also recon- structed to exclude the effects of introns in Hb genes. For qualitative analysis of the genetic structures of local populations, we calculated Materials and Methods the posterior probabilities of assigned individuals into highland and Sampling, sequencing and gametic phasing lowland populations of P. minor in STRUCTURE 2.3.4 (Pritchard To determine the genetic structure of all populations, we employed et al. 2000). We performed STRUCTURE analyses on 4 different a total of 159 samples from all 5 clades of the 3 species (P. minor, sub-datasets: all nuclear introns, the introns plus a -globin gene, the A A A P. cinereus and P. major) within the great tit complex introns plus b -globin gene, and the introns plus a - and b -globin (Supplementary Table S1). All the collections were carried out with genes. The analyses were performed for 1 or 2 population models permission from the State Forestry Administration of the People’s (K ¼ 1 or 2) without prior sampling and conducted with admixture Republic of China and in accordance with the National Wildlife model and independent allele frequencies (a ¼ 1, k ¼ 1). The algo- Conservation Law of China. To explore local population differenti- rithm was run for 100,000 iterations with a burn-in of 25,000 ation in response to changes in elevation, we analyzed a sub-dataset iterations. of P. minor populations. According to the genetic structure parti- We calculated polymorphism and diversity of all loci in each tioned by elevation (Figure 1), we sampled highland (1,000–4,282 m population and in all 3 populations of P. minor combined in DnaSP a.s.l.) and lowland individuals (7–210 m a.s.l.) in 3 local populations v5.10.1 (Pritchard et al. 2000). Tajima’s D (Tajima 1989) and Fu’s of P. minor (N ¼ 66) (Supplementary Table S1). The 3 populations Fs (Fu and Li 1993; Fu 1997) statistics were computed to test devi- were located in the Southwest Mountain region (highland Clade B, ation from neutrality. The significant deviation from neutrality is an HB, N ¼ 23), the alpine region of central China (highland Clade A, outlier of 95% confidence interval in a null model of neutral expect- HA, N ¼ 20), and the plain of southeast China (lowland Clade A, ation in coalescent simulation based on segregated sites and recom- LA, N ¼ 23). The HB and HA populations were adjacent to the bination rate per gene. We calculated pairwise F values for each of ST southeast and the east of the QTP, respectively (Figure 2A). 8 loci and the CDS of Hb genes in ARLEQUIN 3.5. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 ZHu et al.  Elevational divergence in major hemoglobin genes of great tits 3 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 pro- portional to haplotype frequencies. An elevational differentiation presents in a -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. scalars of 0.25 (mtDNA) and 1.00 (autosomal DNA), we initially Genetic diversity and gene flow in local populations performed IMa runs with flat priors for each parameter. On the We used a coalescent-based “isolation-with-migration” model im- basis of preliminary runs, we tuned upper bounds for priors that en- plemented in IMa (Hey and Nielsen 2007) to estimate inter- compassed the overall posterior distributions. We finally ran the population gene flow in the sampled P. minor populations. We analyses in 6 Markov-coupled chains with geometric heating focused on the gene flow strength in the same dataset (a -globin or (g1 ¼ 0.8; g2 ¼ 0.9), discarding at least 20,000 samples as “burn- reference loci) and asymmetric gene flow ratio scaled by migration in”. Each chain was continued until the effective sample size rates (m* ¼ 2 Nm) i.e., locus-specific effective population size rather (ESS) was >100. To ensure the convergence of parameter estimates, than mutation rates. We thus conducted analyses for the a -globin 2 replicated runs in each analysis were conducted by using different gene and joint reference loci (2 mtDNA and 4 autosomal introns, random seeds. 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). Genotypic profiles of hb genes and elevational clines of According to the biasing effects of substitution models, we fitted the hb alleles Hasegawa–Kishino–Yano (HKY) and Infinite Sites (IS) mutation To characterize possible genotypes differentiated between highland models for mtDNA and autosomal DNA regardless of calculated and lowland populations, we estimated genotypic frequencies in models (Strasburg and Rieseberg 2010). Final parameter estimations Hardy–Weinberg equilibrium (HWE) for haplotypes and genotypes of were scaled to effective population size multiplied by mutation rate the 2 Hb genes in ARLEQUIN 3.5 (Excoffier and Lischer 2010). To (h ¼ 4Nl) with estimated migration rates (m ¼ m/l) transformed to explore the threshold of elevational genetic variation, we conducted a 2Nm (M ¼ h  m/2) (Hey and Nielsen 2004). Given inheritance maximum likelihood estimation of Hb allele frequencies by assigning Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 4 Current Zoology, 2017, Vol. 0, No. 0 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 a -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 a -globin gene or neutral loci. (B-C) The genetic structures of major Hb genes and neutral reference loci (mtDNA and nucleotide introns) were illus- trated 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). individuals of P. minor populations into 2 parallel elevational tran- Results sects: T1 in the south and T2 in the north (Supplementary Figure S2). The population genetic structure along elevations We used the cline shapes to acquire the profile of highland replace- For both the great tit superspecies and P. minor, significant genetic ments and a modified Markov Chain Monte Carlo (MCMC) algo- divergence was found in Hb genes, especially for the a -globin gene rithm to determine support limits in ClineFit 2.0b (Porter et al. 1997). (Figures 1 and 2). The divergence of the a -globin gene partitioned The 4-parameter model for the cline center (c), cline width (w), min- the superspecies into 2 elevational clades: the monophyletic clade of imum allele frequency (P ) and maximum allele frequency (P )atthe L R high-altitude individuals from the phylogeographic clade A, B, C, ends of the cline shapes was estimated. We set a burn-in of 800 par- and E, and the other clade including the clade D and low-altitude in- ameter tries per annealing step and 1000 tries for the second run that dividuals from the clade A (Figure 1). Similarly, the network of the are combined in sampling; subsequently, 2000 replicates were saved a -globin gene in P. minor populations presented the monophyletic in the sampling parameter space, with 20 replicates between each highland clade including HA and HB, which was diverged from the save. To assess the coincidence of the cline center and the concordance lowland clade-A population (LA) (Figure 2). The b -globin gene, of the width, we combined the sampled parameter space of replicated however, did not exhibit any clear genetic structure. The single runs according to the maximum 2 LnL limits and subsequently gener- b21Ala/Thr substitution presented in the southwest mountains was ated support plots (not shown) and sigmoidal cline graphics in not fixed in frequency (Figure 2 and 4). Therefore, in further ana- Mathematica v10. The F for each non-synonymous substitution be- ST lyses of genetic variance, we mainly used the local population data- tween highland and lowland populations was calculated using locus- set to explore a potential genetic response to the altitudinal gradient by-locus AMOVA in ARLEQUIN 3.5 (Excoffier and Lischer 2010). in P. minor. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 ZHu et al.  Elevational divergence in major hemoglobin genes of great tits 5 Figure 3. Posterior probability of individual assignment to highland and lowland populations in STRUCTURE 2.3. Black, gray, and white represent individuals A A from population HB, HA and LA (Fig. 2). (A) 4 introns of reference loci, (B) 4 reference introns and a -globin gene, (C) 4 reference introns and b -globin gene, (D) 4 A A reference introns and 2 hemoglobin genes. The higher divergence occurs in the dataset including a -globin gene but lower with b -globin gene add-in. The genetic structure in differences of elevation, which could be populations compared with other loci (average 0.7266 0.21 SD). The largely due to differences in the a -globin gene, was further confirmed nucleotide diversity of both Hb genes and ND2 of clade A exhibited by STRUCTURE analysis on local populations (Figure 3). The different the lowest levels (p< 0.002, average 0.00526 0.0050 SD). Significant assignment to lowland or highland (>1,100 m a.s.l.) populations was negative Tajima’s D and Fu’s F values were present in different popu- A A most dramatic in average posterior probability of the a -globin (6.1% lations of b -globin and reference loci (Supplementary Table S3). In versus 88.2%,68.0% and 18.7% standard deviation [SD]) and less so a -globin gene of all populations, insignificant negative values were A A of joint a -and b -globin genes (8.6% versus 82.6%,610.9% and predominant, indicating lower than expected neutrality. Significant 22.0% SD). In contrast, the assignments showed much less divergence negative values predominated among nuclear loci of the LA popula- A A in response to elevation either for the single b -globin (43.6% and tion and in b -globin for both HA and LA populations, but for 55.3%) or for nuclear introns (55.7% and 44.8%). mtDNA in the HB population. In P. minor populations of clade A, sig- nificantly greater divergence was indicated by fixation indices (F )in ST the a -globin gene in contrast to other loci (F *> 0.5 versus< 0.1, ST Gene flow among elevationally differentiated Table 2). The significant inter-clade divergence of mtDNA was con- populations of P. Minor sistent with genetic structure. The pairwise F values of the a -globin ST The coalescent estimations of IMa showed unbalanced downslope gene showed significant differentiation among populations at different gene flow (from highland to lowland) in the joint reference loci elevations (Table 2). (assumed as neutral gene flow) among local P. minor populations All the Hb genotypes were profiled in HWE analyses (Table 1; Figure 2A). The highest gene flow was found in the down- (Supplementary Table S4). The genotype of each Hb gene exhibited slope (HA to LA) neutral loci and in 16.1, 7.3 and 18.5 times higher differences in frequency between highland and lowland populations. level than the upslope (LA to HA) and the gene flow between 2 clades The 6-allele genotype in the a -globin gene was completely fixed in (HB-LA and HB-HA) (Table 1). For the a -globin gene, in contrast, all individuals in the HB population (the same genotype in 46 downslope gene flow between clades (HB-HA) was 23.9 times higher chromosomes). Accordingly, both the observed genotypic heterozy- than upslope gene flow within the clade A (LA-HA). The lower gene gosity of the a -globin gene in HB was lower than that in the other flow of neutral loci is between 2 clades (HB versus HA/LA), while the 2 populations (H /H ¼ 0.753, 0.924 and 0.861 in HB, HA, and o e lower a -globin gene flow is between populations differed by eleva- LA). The b -globin genotype of most individuals was homozygous, tions (LA versus HB/HA). Relative to the restrained neutral gene flow, and a unique substitution of b21Ala/Thr occurred in HB popula- the much higher gene flow in a -globin from HB into HA is in accord- tions at a low frequency (23.5%). From the allele-specific analysis, ance with the elevational structure of local populations. a 57Gly/Ala (F ¼ 0.004) and b21Ala/Thr (F * ¼ 0.168, ST ST *P< 0.05) exhibited the lowest differentiation between highland Genetic divergence and genotypic profiles in major Hb and lowland populations (Figure 4). The a -globin genotypes of 2 genes of P. Minor distant alleles were fixed according to elevation (F * > 0.9): ST We identified nucleotide polymorphisms of all genes in the P. minor a35Ala/72Asn was predominant in highland populations (HA and A A populations (Supplementary Table S3). The haplotype diversity was HB). Moreover, a 49Asn/Ser and a 108Ala/Val also exhibited high the lowest in ND2 and b -globin gene (H < 0.5) of HA and LA levels of elevational differentiation (F * ¼ 0.845). d ST Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 6 Current Zoology, 2017, Vol. 0, No. 0 Table 1. Coalescent estimation of population genetic parameters among the P. minor populations using IMa Loci Populations h h m m m * m * m */m * H L H L H L L H (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 h ¼ 4Nl, N is the effective population size of each population, l is the mutation rate. m , the migration rate from lowland to highland (upslope); m , 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. m * and m *, were scaled m by effective population size by m  h/2 of upslope and downslope gene flow. The m */m * H L L H ratio is an index of the asymmetric gene flow. NA is an ineffective ratio of zero m values. ST Table 2. Pairwise F values for nucleotide haplotypes of the P. minor populations A A A A Pairwise a a -cds b b -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. Elevational variation in allele frequencies across complete elevational fixation (Supplementary Figure S3). And fre- quencies of the same allele in highland populations exhibited differ- parallel transects ent results in parallel transects that a 49/108 was flexible in T2 Specific a -globin allelic profiles in local populations of P. minor (Figure 5 and Supplementary Figure S3). In contrast, the a 22/35/57 were further determined by cline-fitting analyses along parallel ele- lacked differentiation along the altitudinal gradient for samples vational gradients (Figure 5 and Table 3). The highland alleles from all clades (Figure 1). The b 21Thr occurred at a lower fre- shifted along altitudinal gradients coincide at a threshold of 1,100 m quency from lowland to highland populations except for fixation at a.s.l. (6150 m SD) that is consistent with the genetic structure of the peak elevation (< 36% at lower than 3,000 m a.s.l.). superspecies sorted by elevation. The center fittings of cline shapes were similar across parallel transects at or above 1,000 m. The aver- age 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 Discussion consistent cline centers (1,060, 962–1,152 m a.s.l.) and widths (446, Elevational differentiation of genetic structure in 14–813 m) occurred in a 49/72/108 of the T1 transects, whereas A A the approximately coincident centers of a 22/35/108 were in T2. a -globin gene The average center of T2 was slightly higher (1,131.3 m a.s.l.) and We found an elevationally differentiated structure in the a -globin wider (707 m) than that of the T1 transect. In the more elevationally gene, which is distinctive from the mtDNA phylogeny and the un- fixed alleles (a49/72/108 of T1), higher differentiation (F ¼ 1.00, structured networks in nuclear loci. This topological mismatch be- ST P< 0.05) was apparent across elevational transects. Accordingly, tween functional gene and neutral genes might be induced by the a 49-72-108 across T1 exhibited the highest differentiation incomplete lineage sorting or introgressive hybridization between (F ¼ 1.00), and a 22 across T1 had the lowest differentiation sister lineages, as in sister Anatidae species (Natarajan et al. 2015). ST (F ¼ 0.36; Table 3). According to the genetic structure, the unique The elevationally differentiated genetic structure of the a -globin ST A A b -globin replacement without an elevation pattern could not be gene is also distinctive from the unstructured b -globin. This result used as the marker to reveal the elevation pattern in local popula- differs from most Andean birds that elevational differentiated tions (Figure 4 and Supplementary Figure S3). b -globin gene is contributed to the conspecific adaptive divergence The Hb allelic variations in the great tit complex samples com- (McCracken et al. 2009b; Bulgarella et al. 2012; Galen et al. 2015; A A prised the same 6 a -globin substitutions and 1 b -globin substitu- Natarajan et al. 2015; but see Wilson et al. 2012). The high- and tion as in the P. minor populations. In accordance with local low-altitude divergence in a -globin gene only occurred in cinna- A A populations, the allele frequencies of a 49/72/108 exhibited rela- mon teals as similar manner as b -globin gene divergence in other tively greater fixed elevational variation and a 72 exhibited waterfowl (McCracken et al. 2009b; Bulgarella et al. 2012; Wilson Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 ZHu et al.  Elevational divergence in major hemoglobin genes of great tits 7 Figure 4. Genotype frequencies in amino acid alleles of major Hb genes. F represents for differentiation of each allele between highland and lowland popula- ST A A tions, *P< 0.05. The striking heterozygous highland and lowland genotypes involve a 22Asp/Gly in LA population and b 21Ala/Thr in HB population, in contrast to low-level admixed in a35/57. The more homozygous a 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. et al. 2012). The adaptive divergence in Hb genes of Andean water- population expansion. The a -globin genetic variants may have birds consistently contributed to the dispersal of ancestral genotypes existed in the high-altitude eastern Himalayan population and dis- in lowland populations by colonization into the highlands. In con- persed to the lower habitats in central China earlier before the inter- clusion, the discrepancy in genetic structures between the 2 major clade divergence in neutral loci. The strong downslope neutral gene Hb genes might be resulted from adaptive evolution in the a -globin flow, in contrast, implies the historical rapid population expansion gene of great tits. was from the central highland into the southeast lowland population after inter-clade divergence. Thus, the restricted a -globin gene flow within the same clade may be a result of rapid population expansion Gene flow patterns among local populations of P. minor in the lowland population, offsetting the possible countervailing se- From coalescent gene flow estimation of gene flow in P. minor,we lection on high-altitude genetic variants after inter-clade divergence. found higher than expected a -globin gene flow between clades with restricted neutral gene flow. In contrast to background gene Elevational variation of alleles depends on flow revealed by neutral loci, the inter-clade downslope gene flow of a -globin gene is nearly absent between elevation-diverged popula- demographic history tions within the clade. Similar patterns exhibited in the restricted The genotypic profile of major Hb genes in HWE coincides with the A A b -globin gene flow in the adaptive elevational divergence with genetic structure of elevational segregation in the a -globin gene. highland-derived replacements in Andean waterbirds (McCracken However, the haplotypic H /H ratios of the a -globin gene in the o e et al. 2009b; Wilson et al. 2012; Natarajan et al. 2015). The a -glo- HB and LA populations deviated significantly from equilibrium. bin gene flow from subtropical highlands (HB) into temperate low- This result suggests subunit-specific allelic profiles. From elevational lands (HA) of P. minor is in a direction similar to the b -globin gene cline fitting, we identified the high/lowland fixed a49/72/108 geno- flow from tropical highlands into temperate lowland of south- types that coincide with possible standing variants a49Asn/72Asn/ American waterfowl (McCracken et al. 2009b; Munoz-Fuentes 108Ala, originating from an ancestral population in the eastern et al. 2013). The restricted a -globin gene flow in the same clade of Himalaya (Tietze and Borthakur 2012; Qu et al. 2015). From the P. minor may also have been homogenized by high gene flow of neu- ancestral states reconstruction of Hb genes in Paridae, the highland tral loci contingent on population history, as in the pattern shown a -globin haplotype was inferred to be the standing genetic variant for major Hb genes in Andean waterfowl (McCracken et al. 2009b). in P. minor (Zhu et al., unpublished data). In contrast, the lowland Our results suggest the role of a -globin gene flow in driving genetic fixed b21Thr is more likely a result of newly derived mutations in divergence along elevational gradients in population history of a highland eastern Himalaya populations and fixed after the rapid ex- widespread species. pansion. These results implied that standing genetic variation of the Given the demographic history of the great tit (Zhao et al. 2012; a -globin gene might have been fixed in clade B before the inter- Qu et al. 2015), we propose a plausible route for historical clade divergence. The heterozygous a49/72/108 genotypes in clade Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 8 Current Zoology, 2017, Vol. 0, No. 0 Figure 5. Elevational cline fitting of a -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 popula- tions. The coincident cline varied elevational shifts occurred in a49/72/108 in T2 transect and a22/35 in T1 transect. The genotypic frequency of a49Asn/Ser-72Asn/ His-108Ala/Val in T1 transect and a72Asn/His in both transects are fixed in high/low elevations outside the cline shift regions. Table 3. Maximum likelihood of elevational clines of aA-globin allele frequencies across parallel transects in P. minor Alleles Transect Centre/m Width/m P P F (*P < 0.05) L R ST a 22 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) a 35 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) a 49 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) a 72 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) a 108 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 lnL –lnL< 2 likelihood. max Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 ZHu et al.  Elevational divergence in major hemoglobin genes of great tits 9 A, in contrast, may be a mixed result of countervailing selection on Author contributions standing high-altitude variants and homogenizing neutral gene flow, X.Z., Y.G., and F.L. conceived the ideas; G.S. and Y.Q. collected samples; as the case in Andean waterfowl (McCracken et al. 2009b; Wilson Y.G. and Z.X. performed the experiments and analyzed the data; X.Z., Y.G., et al. 2012; Natarajan et al. 2015;). Therefore, we suggest that the G.S., G.D., and F.L. led the writing; X.Z., Y.G., G.S., Y.Q., D.G. and F.L. a -globin gene is likely under historical adaptive selection rather discussed the results and determined the conclusions. All authors agreed on than the b -globin gene in P. minor and the great tit complex. the main conclusions of the article. It is known that adaptive divergence in a wild species might not be determined by phenotypic changes of Hb , whereas the fixation of Supplementary material genetic variants of functional genes eventually depends on demo- graphic history (Olson-Manning et al. 2012; Raeymaekers et al. 2014; Supplementary material can be found at http://www.cz.oxfordjournals.org/. Ferchaud and Hansen 2016). The fixation of standing genetic vari- ants in highland populations (HB or HA) was inferred for most alleles References (8/10, except for a49/108 in T2) (Zhu et al., unpublished data) that might be homogenized by neutral gene flow. Similarly, much nar- Bandelt H-J, Forster P, Ro ¨ hl A, 1999. Median–joining networks for inferring rower shifts but unfixed highland a49Asn/108Ala alleles in HA and intraspecific phylogenies. Mol Biol Evol 16:37–48. Bulgarella M, Peters JL, Kopuchian C, Valqui T, Wilson RE et al., 2012. wider cline shifts (a22/35/72 of T2) may have been caused by neutral Multilocus coalescent analysis of haemoglobin differentiation between low– gene flow in clade A. In agreement with these results, the a22/35 were and high–altitude populations of crested ducks Lophonetta specularioides. absent from altitudinal clines in the full sampling, whereas only the Mol Ecol 21:350–368. completely diverged a72Asn/His suggested that an ancestral standing Cheviron ZA, Natarajan C, Projecto-Garcia J, Eddy DK, Jones J et al., 2014. polymorphism was fixed in all highland populations of the great tit Integrating evolutionary and functional tests of adaptive hypotheses: a case complex. Therefore, although the same elevational differentiation was study of altitudinal differentiation in hemoglobin function in an Andean detected in the genetic structure of the superspecies, specific allele pro- sparrow Zonotrichia capensis. Mol Biol Evol 31:2948–2962. files of in P. minor populations were affected by local homogenizing Colosimo PF, Hosemann KE, Balabhadra S, Villarreal G, Dickson M et al., gene flow. In a word, differed from the case of standing Hb variants 2005. Widespread parallel evolution in sticklebacks by repeated fixation of originated in low-altitude Andean waterfowl (Natarajan et al. 2015), ectodysplasin alleles. Science 307:1928–1933. Excoffier L, Lischer HEL, 2010. Arlequin suite ver 3.5: a new series of pro- the recurrently replaced alleles in deeply diverged lineages of Eastern grams to perform population genetics analyses under Linux and Windows. Himalaya and Central Asia can be interpreted as the fixation of stand- Mol Ecol Resour 10:564–567. ing high-altitude alleles in East-Himalaya. Ferchaud AL, Hansen MM, 2016. The impact of selection, gene flow and In summary, our study highlights the necessity of integrating demographic history on heterogeneous genomic divergence: three-spine gene flow and demographic history into exploring potential adap- sticklebacks in divergent environments. Mol Ecol 25:238–259. tive evolution in the genetic divergence of widespread species. This Fu YX, 1997. Statistical tests of neutrality of mutations against population preliminary investigation on genetic variants in a widespread passer- growth, hitchhiking and background selection. Genetics 147:915–925. ine superspecies provides an ideal model to further explore the adap- Fu YX, Li WH, 1993. Statistical tests of neutrality of mutations. Genetics 133: tive evolution of Hb genes in wild birds. 693–709. Galen SC, Natarajan C, Moriyama H, Weber RE, Fago A et al., 2015. Contribution of a mutational hot spot to hemoglobin adaptation in high- altitude Andean house wrens. Proc Natl Acad Sci USA 112: 13958–13963. Gill F, Donsker D, 2016. Ioc world bird list (v 6.4) [Cited 2016 May 21] Data archiving Available from: http://www.worldbirdnames.org/ All DNA sequences are accessible under GenBank accession num- Gosler A, Clement P, Christie DA, 2017. Great tit Parus major. In: del Hoyo J, bers (MF321902–MF322504), and sample information for all indi- Elliott A, Sargatal J, Christie DA, Juana d, editors. Hand Book of the Birds viduals is found in Table S1 of the Supplementary file. of the World Alive. Barcelona: Lynx Edicions. Hey J, Nielsen R, 2004. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. Persimilis. Genetics 167: 747–760. Acknowledgments Hey J, Nielsen R, 2007. Integration within the Felsenstein equation for im- proved Markov chain Monte Carlo methods in population genetics. Proc We thank Shimiao Shao for providing comments on the early draft, Yalin Natl Acad Sci USA 104:2785–2790. Chen and Dezhi Zhang for revising the latest version of this manuscript. We Kumar A, Natarajan C, Moriyama H, Witt CC, Weber RE et al., 2017. also thank Per Alstro ¨ m for valuable comments on this work. We thank Stability-mediated epistasis restricts accessible mutational pathways in the Tianlong Cai for providing the elevation maps and Ruiying Zhang for provid- functional evolution of avian hemoglobin. Mol Biol Evol 34:1240–1251. ing sample information. Kvist L, Arbabi T, P€ ackert M, Orell M, Martens J, 2007. Population differenti- ation in the marginal populations of the great tit (Paridae: Parus major). Biol J Linn Soc 90:201–210. Kvist L, Martens J, Higuchi H, Nazarenko AA, Valchuk OP et al., 2003. Funding Evolution and genetic structure of the great tit Parus major complex. Proc R This work was funded by the Strategic Priority Research Programme of the Soc Lond B Biol Sci 270:1447–1454. Chinese Academy of Sciences (XDB13020300), State Key Programme of Librado P, Rozas J, 2009. Dnasp v5: a software for comprehensive analysis of Natural Science Foundation of China (NSFC 31330073), a grant from the DNA polymorphism data. Bioinformatics 25:1451–1452. Ministry of Science and Technology of China (2014FY210200) to F.L., and a Mairbaurl H, Weber RE, 2012. Oxygen transport by hemoglobin. Compr grant from the NSFC programme (J1210002) to Y.G. and X.J.Z. G.S. is sup- Physiol 2:1463–1489. ported by the NSFC fund (31471991) and a grant (Y229YX5105) from the McCracken KG, Barger CP, Bulgarella M, Johnson KP, Kuhner MK et al., Key Laboratory of the Zoological Systematics and Evolution of the Chinese 2009a. Signatures of high-altitude adaptation in the major hemoglobin of Academy of Sciences. five species of Andean dabbling ducks. Am Nat 174:631–650. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 10 Current Zoology, 2017, Vol. 0, No. 0 McCracken KG, Bulgarella M, Johnson KP, Kuhner MK, Trucco J et al., Qu Y, Tian S, Han N, Zhao H, Gao B et al., 2015. Genetic responses to sea- 2009b. Gene flow in the face of countervailing selection: adaptation to high- sonal variation in altitudinal stress: whole-genome resequencing of great tit altitude hypoxia in the bA hemoglobin subunit of yellow-billed pintails in in eastern Himalayas. Sci Rep 5:14256. the Andes. Mol Biol Evol 26: 815–827. Raeymaekers JAM, Konijnendijk N, Larmuseau MHD, Hellemans B, De Monge C, Leon-Velarde F, 1991. Physiological adaptation to high altitude: Meester L et al., 2014. A gene with major phenotypic effects as a target for oxygen transport in mammals and birds. Physiol Rev 71: 1135–1172. selection vs. homogenizing gene flow. Mol Ecol 23:162–181. Munoz-Fuentes V, Cortazar-Chinarro M, Lozano-Jaramillo M, McCracken Schluter D, Clifford EA, Nemethy M, McKinnon JS, 2004. Parallel evolution KG, 2013. Stepwise colonization of the Andes by ruddy ducks and the evo- and inheritance of quantitative traits. Am Nat 163:809–822. lution of novel beta-globin variants. Mol Ecol 22: 1231–1249. Schluter D, Conte GL, 2009. Genetics and ecological speciation. Proc Natl Natarajan C, Projecto-Garcia J, Moriyama H, Weber RE, Munoz-Fuentes V Acad Sci USA 106 (1 Suppl): 9955–9962. et al., 2015. Convergent evolution of hemoglobin function in high-altitude Scott GR, Milsom WK, 2006. Flying high: a theoretical analysis of the factors Andean waterfowl involves limited parallelism at the molecular sequence limiting exercise performance in birds at altitude. Respir Physiol Neurobiol level. PLoS Genet 11:e1005681. 154:284–301. Olson-Manning CF, Wagner MR, Mitchell-Olds T, 2012. Adaptive evolution: Stephens M, Smith NJ, Donnelly P, 2001. A new statistical method for haplo- Evaluating empirical support for theoretical predictions. Nat Rev Genet type reconstruction from population data. Am J Hum Genet 68:978–989. 13:867–877. Storz JF, 2007. Hemoglobin function and physiological adaptation to hypoxia Opazo JC, Hoffmann FG, Natarajan C, Witt CC, Berenbrink M et al., 2015. in high-altitude mammals. J Mammal 88:8. Gene turnover in the avian globin gene families and evolutionary changes in Strasburg JL, Rieseberg LH, 2010. How robust are ‘isolation with migration’ hemoglobin isoform expression. Mol Biol Evol 32:871–887. analyses to violations of the IM model? A simulation study. Mol Biol Evol P€ ackert M, Martens J, Eck S, Nazarenko AA, Valchuk OP et al., 2005. The great 27:297–310. tit Parus major: a misclassified ring species. Biol J Linn Soc 86:153–174. Tietze DT, Borthakur U, 2012. Historical biogeography of tits (Aves: Paridae, Perutz MF, 1983. Species adaptation in a protein molecule. Mol Biol Evol Remizidae). Org Divers Evol (12):433–444. 1:1–28. Weber RE, 2007. High-ltitude adaptations in vertebrate hemoglobins. Respir Porter AH, Wenger R, Geiger H, Scholl A, Shapiro AM, 1997. The Pontiac Physiol Neurobiol 158: 132–142. daplidice - edusa hybrid zone in northwestern Italy. Evolution Wilson RE, Peters JL, McCracken KG, 2012. Genetic and phenotypic diver- 51:1561–1573. gence between low– and high–altitude populations of two recently diverged Pritchard JK, Stephens M, Donnelly P, 2000. Inference of population structure cinnamon teal subspecies. Evolution 67: 170–184. using multilocus genotype data. Genetics 155:945–959. Zhao N, Dai C, Wang W, Zhang R, Qu Y et al., 2012. Pleistocene climate Projecto-Garcia J, Natarajan C, Moriyama H, Weber RE, Fago A et al., 2013. changes shaped the divergence and demography of Asian populations of the Repeated elevational transitions in hemoglobin function during the evolution great tit Parus major: evidence from phylogeographic analysis and ecolo- of Andean hummingbirds. Proc Natl Acad Sci U S A (51):20669–20674. gical niche models. J Avian Biol 43: 297–310. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Current Zoology Oxford University Press

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

Free
10 pages

Loading next page...
 
/lp/ou_press/elevational-divergence-in-the-great-tit-complex-revealed-by-major-kYmJqRNqkf
Publisher
Editorial Office
Copyright
© The Author (2017). Published by Oxford University Press.
ISSN
1674-5507
eISSN
2396-9814
D.O.I.
10.1093/cz/zox042
Publisher site
See Article on Publisher Site

Abstract

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 b -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.) A A 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 expan- sion. In concordance with patterns of historical gene flow, genotypic and allelic profiles provide dis- tinctive patterns of fixation in different high-altitude populations. The fixation of alleles at contrast- ing 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. Key words: elevational divergence, gene flow, genetic structure, hemoglobin gene. Widespread species in hierarchical population structure can provide natural selection in native vertebrates under hypoxia stress (Perutz opportunities for investigating the potential influence of demo- 1983; Monge and Leon-Velarde 1991; Scott and Milsom 2006; graphic history in shaping genetic variation arising in response to Storz 2007; Weber 2007). However, intraspecific genetic divergence environmental change. A comprehensive example is a fitness- in Hb genes may not necessarily cause adaptive changes in Hb in correlated morphological gene in the adaptive evolution of marine birds (McCracken et al. 2009a; Munoz-Fuentes et al. 2013; and freshwater three-spine stickleback Gasterosteus aculeatus re- Cheviron et al. 2014; Galen et al. 2015; Natarajan et al. 2015). A sponding to gene flow and demographic history (Schluter et al. hypothesis proposed that highland adaptation after colonization 2004; Colosimo et al. 2005; Schluter and Conte 2009; Raeymaekers from lowlands following divergence based on the Hb gene flow et al. 2014). For high-altitude adaptive evolution, genetic variation from lowland to highland in Andean waterbirds Anatidae of oxygen-carrying hemoglobin (Hb) has long been used to track (McCracken et al. 2009b; Wilson et al. 2012). The hypothesis was V C The Author (2017). Published by Oxford University Press. 1 This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 2 Current Zoology, 2017, Vol. 0, No. 0 A A subsequently supported by highland genotypes derived from ances- We amplified and sequenced a - and b -globin genes and 6 un- tral lowland genotypes in functional experiments (Natarajan et al. linked reference loci (2 mtDNA and 4 autosomal introns) for all 2015). For avian species at broad elevations, whether the genetic samples. The reference loci are widely used neutral markers in gen- structure and gene flow in Hb genes can reveal intraspecific adaptive eral phylogenetic studies, including the mitochondrial control region differentiation in elevations is still an open question. (CR, 539 bp), NADH dehydrogenase subunit 2 (ND2, 737 bp), To investigate the role of gene flow and demographic history in beta-fibrinogen intron 5 (Fib, 281 bp), transforming growth factor elevational genetic divergence, we performed population genetic ana- beta 2 intron 5 (TGFB2, 334 bp), myoglobin intron 2 (Myo, 556 bp) lyses on major Hb genes and unlinked neutral loci in a widespread and ornithine decarboxylase intron 7 (ODC, 284 bp). We designed species complex, the great tit. The HbA tetramer, which is the major subunit-specific primers for Hb genes (660 bp of a -globin gene and Hb isoform expressed in adult birds, comprises 2 a-chains and 2 1,274 bp of b -globin gene). After DNA extraction from venous blood A A b-chains encoded by a -and b -globin gene, respectively (Mairbaurl or pectoral muscle with a DNeasy Tissue kit (QIAGEN, Hilden, and Weber 2012; Opazo et al. 2015). It has been revealed that the Germany), we amplified the 8 genes with KD-Plus DNA polymerase adaptive genetic variants mostly presented in the b -globin gene in (TransGen Biotech, Beijing, China). The general PCR thermal profile high-altitude Andean birds (Projecto-Garcia et al. 2013; Galen et al. was 5 min at 94 C, 35 cycles of 40 s at 94 C, 40 s at 52–64 C, 1 min at 2015; Natarajan et al. 2015). In contrast, rare causative variants were 72 C, and 10 min at 72 C (with varied annealing temperatures, as found in the a -globin gene of cinnamon teals Anas cyanoptera and shown in Supplementary Table S2). For each gene, we sequenced the nightjars Hydropsalis ( Wilson et al. 2012; Natarajan et al. 2015; PCR products on the ABI3730 sequencing system (BGI, Beijing, China). Kumar et al. 2017). The great tit complex, mainly comprising Parus For heterozygous autosomal genes, we inferred the gametic major, P. minor,and P. cinereus, is known as a superspecies encircling phases using PHASE v2.1 (Stephens et al. 2001). The PHASE algo- the Qinghai–Tibet Plateau (QTP) and Central Asian desert (Kvist rithm was run 5 times using random seeds to obtain stable results. et al. 2003; P€ ackert et al. 2005; Kvist et al. 2007; Gill and Donsker For the a -globin gene (660 bp) and 4 autosomal introns mentioned 2016). This superspecies is widely distributed from sea level to above above, we inferred the gametic phases of the complete sequences. 4,000 m a.s.l. in Eurasia, with the highest habitats located in the east- For the b -globin gene (1,274 bp), we inferred the gametic phases ern Himalaya (Gosler et al. 2017). A monophyletic clade of P. minor for the fragments spanning from exon 1 to exon 2 and partial intron in eastern Himalaya (P. m. tibetanus and P. m. subtibetanus, clade B) 2 (453 bp), which had no recombination detected in DnaSP v5.10.1 has been further identified to compose the 5-clade mtDNA phylogeny (Librado and Rozas 2009). All substitutions that occurred in the of the great tit complex (Zhao et al. 2012). The demographic history exon 3 (126 bp) of the b -globin gene were synonymous. The se- of P. minor revealed historically stable growth in the eastern quences with high posterior probabilities (>95%) and the consensus Himalaya populations before the last glacial maximum (LGM) and sequence of the cloned samples unresolved in PHASE proceeded into rapid population expansion in the East Asia during the LGM (Zhao downstream analyses. All sequences are accessible under accession et al. 2012; Qu et al. 2015). numbers MF321902-MF322504 in GenBank, NCBI, NIH In the present study, we conducted population genetic analyses (Supplementary Table S2). for exploring the role of gene flow and demographic history in Hb genetic differentiation corresponding to changes in elevation in great Comparison of genealogical structures between hb tit populations. This study mainly focused on 1) comparing the dif- genes and neutral loci ference in genetic structure of Hb genes with neutral loci; 2) examin- To understand the genealogical discordance between Hb genes and ing gene flow among local populations at hierarchical elevations of neutral loci, we reconstructed unrooted networks of haplotypes for P. minor; 3) identifying genotypes and profiles along elevational gra- each gene in local and all samples by using the median-joining algo- dients and 4) elevational variation of allelic frequencies in Hb genes rithm in NETWORK 4.6.1 (Bandelt et al. 1999). Networks based across parallel transects. on the coding sequence (CDS) of the Hb genes were also recon- structed to exclude the effects of introns in Hb genes. For qualitative analysis of the genetic structures of local populations, we calculated Materials and Methods the posterior probabilities of assigned individuals into highland and Sampling, sequencing and gametic phasing lowland populations of P. minor in STRUCTURE 2.3.4 (Pritchard To determine the genetic structure of all populations, we employed et al. 2000). We performed STRUCTURE analyses on 4 different a total of 159 samples from all 5 clades of the 3 species (P. minor, sub-datasets: all nuclear introns, the introns plus a -globin gene, the A A A P. cinereus and P. major) within the great tit complex introns plus b -globin gene, and the introns plus a - and b -globin (Supplementary Table S1). All the collections were carried out with genes. The analyses were performed for 1 or 2 population models permission from the State Forestry Administration of the People’s (K ¼ 1 or 2) without prior sampling and conducted with admixture Republic of China and in accordance with the National Wildlife model and independent allele frequencies (a ¼ 1, k ¼ 1). The algo- Conservation Law of China. To explore local population differenti- rithm was run for 100,000 iterations with a burn-in of 25,000 ation in response to changes in elevation, we analyzed a sub-dataset iterations. of P. minor populations. According to the genetic structure parti- We calculated polymorphism and diversity of all loci in each tioned by elevation (Figure 1), we sampled highland (1,000–4,282 m population and in all 3 populations of P. minor combined in DnaSP a.s.l.) and lowland individuals (7–210 m a.s.l.) in 3 local populations v5.10.1 (Pritchard et al. 2000). Tajima’s D (Tajima 1989) and Fu’s of P. minor (N ¼ 66) (Supplementary Table S1). The 3 populations Fs (Fu and Li 1993; Fu 1997) statistics were computed to test devi- were located in the Southwest Mountain region (highland Clade B, ation from neutrality. The significant deviation from neutrality is an HB, N ¼ 23), the alpine region of central China (highland Clade A, outlier of 95% confidence interval in a null model of neutral expect- HA, N ¼ 20), and the plain of southeast China (lowland Clade A, ation in coalescent simulation based on segregated sites and recom- LA, N ¼ 23). The HB and HA populations were adjacent to the bination rate per gene. We calculated pairwise F values for each of ST southeast and the east of the QTP, respectively (Figure 2A). 8 loci and the CDS of Hb genes in ARLEQUIN 3.5. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 ZHu et al.  Elevational divergence in major hemoglobin genes of great tits 3 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 pro- portional to haplotype frequencies. An elevational differentiation presents in a -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. scalars of 0.25 (mtDNA) and 1.00 (autosomal DNA), we initially Genetic diversity and gene flow in local populations performed IMa runs with flat priors for each parameter. On the We used a coalescent-based “isolation-with-migration” model im- basis of preliminary runs, we tuned upper bounds for priors that en- plemented in IMa (Hey and Nielsen 2007) to estimate inter- compassed the overall posterior distributions. We finally ran the population gene flow in the sampled P. minor populations. We analyses in 6 Markov-coupled chains with geometric heating focused on the gene flow strength in the same dataset (a -globin or (g1 ¼ 0.8; g2 ¼ 0.9), discarding at least 20,000 samples as “burn- reference loci) and asymmetric gene flow ratio scaled by migration in”. Each chain was continued until the effective sample size rates (m* ¼ 2 Nm) i.e., locus-specific effective population size rather (ESS) was >100. To ensure the convergence of parameter estimates, than mutation rates. We thus conducted analyses for the a -globin 2 replicated runs in each analysis were conducted by using different gene and joint reference loci (2 mtDNA and 4 autosomal introns, random seeds. 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). Genotypic profiles of hb genes and elevational clines of According to the biasing effects of substitution models, we fitted the hb alleles Hasegawa–Kishino–Yano (HKY) and Infinite Sites (IS) mutation To characterize possible genotypes differentiated between highland models for mtDNA and autosomal DNA regardless of calculated and lowland populations, we estimated genotypic frequencies in models (Strasburg and Rieseberg 2010). Final parameter estimations Hardy–Weinberg equilibrium (HWE) for haplotypes and genotypes of were scaled to effective population size multiplied by mutation rate the 2 Hb genes in ARLEQUIN 3.5 (Excoffier and Lischer 2010). To (h ¼ 4Nl) with estimated migration rates (m ¼ m/l) transformed to explore the threshold of elevational genetic variation, we conducted a 2Nm (M ¼ h  m/2) (Hey and Nielsen 2004). Given inheritance maximum likelihood estimation of Hb allele frequencies by assigning Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 4 Current Zoology, 2017, Vol. 0, No. 0 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 a -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 a -globin gene or neutral loci. (B-C) The genetic structures of major Hb genes and neutral reference loci (mtDNA and nucleotide introns) were illus- trated 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). individuals of P. minor populations into 2 parallel elevational tran- Results sects: T1 in the south and T2 in the north (Supplementary Figure S2). The population genetic structure along elevations We used the cline shapes to acquire the profile of highland replace- For both the great tit superspecies and P. minor, significant genetic ments and a modified Markov Chain Monte Carlo (MCMC) algo- divergence was found in Hb genes, especially for the a -globin gene rithm to determine support limits in ClineFit 2.0b (Porter et al. 1997). (Figures 1 and 2). The divergence of the a -globin gene partitioned The 4-parameter model for the cline center (c), cline width (w), min- the superspecies into 2 elevational clades: the monophyletic clade of imum allele frequency (P ) and maximum allele frequency (P )atthe L R high-altitude individuals from the phylogeographic clade A, B, C, ends of the cline shapes was estimated. We set a burn-in of 800 par- and E, and the other clade including the clade D and low-altitude in- ameter tries per annealing step and 1000 tries for the second run that dividuals from the clade A (Figure 1). Similarly, the network of the are combined in sampling; subsequently, 2000 replicates were saved a -globin gene in P. minor populations presented the monophyletic in the sampling parameter space, with 20 replicates between each highland clade including HA and HB, which was diverged from the save. To assess the coincidence of the cline center and the concordance lowland clade-A population (LA) (Figure 2). The b -globin gene, of the width, we combined the sampled parameter space of replicated however, did not exhibit any clear genetic structure. The single runs according to the maximum 2 LnL limits and subsequently gener- b21Ala/Thr substitution presented in the southwest mountains was ated support plots (not shown) and sigmoidal cline graphics in not fixed in frequency (Figure 2 and 4). Therefore, in further ana- Mathematica v10. The F for each non-synonymous substitution be- ST lyses of genetic variance, we mainly used the local population data- tween highland and lowland populations was calculated using locus- set to explore a potential genetic response to the altitudinal gradient by-locus AMOVA in ARLEQUIN 3.5 (Excoffier and Lischer 2010). in P. minor. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 ZHu et al.  Elevational divergence in major hemoglobin genes of great tits 5 Figure 3. Posterior probability of individual assignment to highland and lowland populations in STRUCTURE 2.3. Black, gray, and white represent individuals A A from population HB, HA and LA (Fig. 2). (A) 4 introns of reference loci, (B) 4 reference introns and a -globin gene, (C) 4 reference introns and b -globin gene, (D) 4 A A reference introns and 2 hemoglobin genes. The higher divergence occurs in the dataset including a -globin gene but lower with b -globin gene add-in. The genetic structure in differences of elevation, which could be populations compared with other loci (average 0.7266 0.21 SD). The largely due to differences in the a -globin gene, was further confirmed nucleotide diversity of both Hb genes and ND2 of clade A exhibited by STRUCTURE analysis on local populations (Figure 3). The different the lowest levels (p< 0.002, average 0.00526 0.0050 SD). Significant assignment to lowland or highland (>1,100 m a.s.l.) populations was negative Tajima’s D and Fu’s F values were present in different popu- A A most dramatic in average posterior probability of the a -globin (6.1% lations of b -globin and reference loci (Supplementary Table S3). In versus 88.2%,68.0% and 18.7% standard deviation [SD]) and less so a -globin gene of all populations, insignificant negative values were A A of joint a -and b -globin genes (8.6% versus 82.6%,610.9% and predominant, indicating lower than expected neutrality. Significant 22.0% SD). In contrast, the assignments showed much less divergence negative values predominated among nuclear loci of the LA popula- A A in response to elevation either for the single b -globin (43.6% and tion and in b -globin for both HA and LA populations, but for 55.3%) or for nuclear introns (55.7% and 44.8%). mtDNA in the HB population. In P. minor populations of clade A, sig- nificantly greater divergence was indicated by fixation indices (F )in ST the a -globin gene in contrast to other loci (F *> 0.5 versus< 0.1, ST Gene flow among elevationally differentiated Table 2). The significant inter-clade divergence of mtDNA was con- populations of P. Minor sistent with genetic structure. The pairwise F values of the a -globin ST The coalescent estimations of IMa showed unbalanced downslope gene showed significant differentiation among populations at different gene flow (from highland to lowland) in the joint reference loci elevations (Table 2). (assumed as neutral gene flow) among local P. minor populations All the Hb genotypes were profiled in HWE analyses (Table 1; Figure 2A). The highest gene flow was found in the down- (Supplementary Table S4). The genotype of each Hb gene exhibited slope (HA to LA) neutral loci and in 16.1, 7.3 and 18.5 times higher differences in frequency between highland and lowland populations. level than the upslope (LA to HA) and the gene flow between 2 clades The 6-allele genotype in the a -globin gene was completely fixed in (HB-LA and HB-HA) (Table 1). For the a -globin gene, in contrast, all individuals in the HB population (the same genotype in 46 downslope gene flow between clades (HB-HA) was 23.9 times higher chromosomes). Accordingly, both the observed genotypic heterozy- than upslope gene flow within the clade A (LA-HA). The lower gene gosity of the a -globin gene in HB was lower than that in the other flow of neutral loci is between 2 clades (HB versus HA/LA), while the 2 populations (H /H ¼ 0.753, 0.924 and 0.861 in HB, HA, and o e lower a -globin gene flow is between populations differed by eleva- LA). The b -globin genotype of most individuals was homozygous, tions (LA versus HB/HA). Relative to the restrained neutral gene flow, and a unique substitution of b21Ala/Thr occurred in HB popula- the much higher gene flow in a -globin from HB into HA is in accord- tions at a low frequency (23.5%). From the allele-specific analysis, ance with the elevational structure of local populations. a 57Gly/Ala (F ¼ 0.004) and b21Ala/Thr (F * ¼ 0.168, ST ST *P< 0.05) exhibited the lowest differentiation between highland Genetic divergence and genotypic profiles in major Hb and lowland populations (Figure 4). The a -globin genotypes of 2 genes of P. Minor distant alleles were fixed according to elevation (F * > 0.9): ST We identified nucleotide polymorphisms of all genes in the P. minor a35Ala/72Asn was predominant in highland populations (HA and A A populations (Supplementary Table S3). The haplotype diversity was HB). Moreover, a 49Asn/Ser and a 108Ala/Val also exhibited high the lowest in ND2 and b -globin gene (H < 0.5) of HA and LA levels of elevational differentiation (F * ¼ 0.845). d ST Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 6 Current Zoology, 2017, Vol. 0, No. 0 Table 1. Coalescent estimation of population genetic parameters among the P. minor populations using IMa Loci Populations h h m m m * m * m */m * H L H L H L L H (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 h ¼ 4Nl, N is the effective population size of each population, l is the mutation rate. m , the migration rate from lowland to highland (upslope); m , 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. m * and m *, were scaled m by effective population size by m  h/2 of upslope and downslope gene flow. The m */m * H L L H ratio is an index of the asymmetric gene flow. NA is an ineffective ratio of zero m values. ST Table 2. Pairwise F values for nucleotide haplotypes of the P. minor populations A A A A Pairwise a a -cds b b -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. Elevational variation in allele frequencies across complete elevational fixation (Supplementary Figure S3). And fre- quencies of the same allele in highland populations exhibited differ- parallel transects ent results in parallel transects that a 49/108 was flexible in T2 Specific a -globin allelic profiles in local populations of P. minor (Figure 5 and Supplementary Figure S3). In contrast, the a 22/35/57 were further determined by cline-fitting analyses along parallel ele- lacked differentiation along the altitudinal gradient for samples vational gradients (Figure 5 and Table 3). The highland alleles from all clades (Figure 1). The b 21Thr occurred at a lower fre- shifted along altitudinal gradients coincide at a threshold of 1,100 m quency from lowland to highland populations except for fixation at a.s.l. (6150 m SD) that is consistent with the genetic structure of the peak elevation (< 36% at lower than 3,000 m a.s.l.). superspecies sorted by elevation. The center fittings of cline shapes were similar across parallel transects at or above 1,000 m. The aver- age 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 Discussion consistent cline centers (1,060, 962–1,152 m a.s.l.) and widths (446, Elevational differentiation of genetic structure in 14–813 m) occurred in a 49/72/108 of the T1 transects, whereas A A the approximately coincident centers of a 22/35/108 were in T2. a -globin gene The average center of T2 was slightly higher (1,131.3 m a.s.l.) and We found an elevationally differentiated structure in the a -globin wider (707 m) than that of the T1 transect. In the more elevationally gene, which is distinctive from the mtDNA phylogeny and the un- fixed alleles (a49/72/108 of T1), higher differentiation (F ¼ 1.00, structured networks in nuclear loci. This topological mismatch be- ST P< 0.05) was apparent across elevational transects. Accordingly, tween functional gene and neutral genes might be induced by the a 49-72-108 across T1 exhibited the highest differentiation incomplete lineage sorting or introgressive hybridization between (F ¼ 1.00), and a 22 across T1 had the lowest differentiation sister lineages, as in sister Anatidae species (Natarajan et al. 2015). ST (F ¼ 0.36; Table 3). According to the genetic structure, the unique The elevationally differentiated genetic structure of the a -globin ST A A b -globin replacement without an elevation pattern could not be gene is also distinctive from the unstructured b -globin. This result used as the marker to reveal the elevation pattern in local popula- differs from most Andean birds that elevational differentiated tions (Figure 4 and Supplementary Figure S3). b -globin gene is contributed to the conspecific adaptive divergence The Hb allelic variations in the great tit complex samples com- (McCracken et al. 2009b; Bulgarella et al. 2012; Galen et al. 2015; A A prised the same 6 a -globin substitutions and 1 b -globin substitu- Natarajan et al. 2015; but see Wilson et al. 2012). The high- and tion as in the P. minor populations. In accordance with local low-altitude divergence in a -globin gene only occurred in cinna- A A populations, the allele frequencies of a 49/72/108 exhibited rela- mon teals as similar manner as b -globin gene divergence in other tively greater fixed elevational variation and a 72 exhibited waterfowl (McCracken et al. 2009b; Bulgarella et al. 2012; Wilson Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 ZHu et al.  Elevational divergence in major hemoglobin genes of great tits 7 Figure 4. Genotype frequencies in amino acid alleles of major Hb genes. F represents for differentiation of each allele between highland and lowland popula- ST A A tions, *P< 0.05. The striking heterozygous highland and lowland genotypes involve a 22Asp/Gly in LA population and b 21Ala/Thr in HB population, in contrast to low-level admixed in a35/57. The more homozygous a 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. et al. 2012). The adaptive divergence in Hb genes of Andean water- population expansion. The a -globin genetic variants may have birds consistently contributed to the dispersal of ancestral genotypes existed in the high-altitude eastern Himalayan population and dis- in lowland populations by colonization into the highlands. In con- persed to the lower habitats in central China earlier before the inter- clusion, the discrepancy in genetic structures between the 2 major clade divergence in neutral loci. The strong downslope neutral gene Hb genes might be resulted from adaptive evolution in the a -globin flow, in contrast, implies the historical rapid population expansion gene of great tits. was from the central highland into the southeast lowland population after inter-clade divergence. Thus, the restricted a -globin gene flow within the same clade may be a result of rapid population expansion Gene flow patterns among local populations of P. minor in the lowland population, offsetting the possible countervailing se- From coalescent gene flow estimation of gene flow in P. minor,we lection on high-altitude genetic variants after inter-clade divergence. found higher than expected a -globin gene flow between clades with restricted neutral gene flow. In contrast to background gene Elevational variation of alleles depends on flow revealed by neutral loci, the inter-clade downslope gene flow of a -globin gene is nearly absent between elevation-diverged popula- demographic history tions within the clade. Similar patterns exhibited in the restricted The genotypic profile of major Hb genes in HWE coincides with the A A b -globin gene flow in the adaptive elevational divergence with genetic structure of elevational segregation in the a -globin gene. highland-derived replacements in Andean waterbirds (McCracken However, the haplotypic H /H ratios of the a -globin gene in the o e et al. 2009b; Wilson et al. 2012; Natarajan et al. 2015). The a -glo- HB and LA populations deviated significantly from equilibrium. bin gene flow from subtropical highlands (HB) into temperate low- This result suggests subunit-specific allelic profiles. From elevational lands (HA) of P. minor is in a direction similar to the b -globin gene cline fitting, we identified the high/lowland fixed a49/72/108 geno- flow from tropical highlands into temperate lowland of south- types that coincide with possible standing variants a49Asn/72Asn/ American waterfowl (McCracken et al. 2009b; Munoz-Fuentes 108Ala, originating from an ancestral population in the eastern et al. 2013). The restricted a -globin gene flow in the same clade of Himalaya (Tietze and Borthakur 2012; Qu et al. 2015). From the P. minor may also have been homogenized by high gene flow of neu- ancestral states reconstruction of Hb genes in Paridae, the highland tral loci contingent on population history, as in the pattern shown a -globin haplotype was inferred to be the standing genetic variant for major Hb genes in Andean waterfowl (McCracken et al. 2009b). in P. minor (Zhu et al., unpublished data). In contrast, the lowland Our results suggest the role of a -globin gene flow in driving genetic fixed b21Thr is more likely a result of newly derived mutations in divergence along elevational gradients in population history of a highland eastern Himalaya populations and fixed after the rapid ex- widespread species. pansion. These results implied that standing genetic variation of the Given the demographic history of the great tit (Zhao et al. 2012; a -globin gene might have been fixed in clade B before the inter- Qu et al. 2015), we propose a plausible route for historical clade divergence. The heterozygous a49/72/108 genotypes in clade Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 8 Current Zoology, 2017, Vol. 0, No. 0 Figure 5. Elevational cline fitting of a -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 popula- tions. The coincident cline varied elevational shifts occurred in a49/72/108 in T2 transect and a22/35 in T1 transect. The genotypic frequency of a49Asn/Ser-72Asn/ His-108Ala/Val in T1 transect and a72Asn/His in both transects are fixed in high/low elevations outside the cline shift regions. Table 3. Maximum likelihood of elevational clines of aA-globin allele frequencies across parallel transects in P. minor Alleles Transect Centre/m Width/m P P F (*P < 0.05) L R ST a 22 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) a 35 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) a 49 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) a 72 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) a 108 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 lnL –lnL< 2 likelihood. max Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 ZHu et al.  Elevational divergence in major hemoglobin genes of great tits 9 A, in contrast, may be a mixed result of countervailing selection on Author contributions standing high-altitude variants and homogenizing neutral gene flow, X.Z., Y.G., and F.L. conceived the ideas; G.S. and Y.Q. collected samples; as the case in Andean waterfowl (McCracken et al. 2009b; Wilson Y.G. and Z.X. performed the experiments and analyzed the data; X.Z., Y.G., et al. 2012; Natarajan et al. 2015;). Therefore, we suggest that the G.S., G.D., and F.L. led the writing; X.Z., Y.G., G.S., Y.Q., D.G. and F.L. a -globin gene is likely under historical adaptive selection rather discussed the results and determined the conclusions. All authors agreed on than the b -globin gene in P. minor and the great tit complex. the main conclusions of the article. It is known that adaptive divergence in a wild species might not be determined by phenotypic changes of Hb , whereas the fixation of Supplementary material genetic variants of functional genes eventually depends on demo- graphic history (Olson-Manning et al. 2012; Raeymaekers et al. 2014; Supplementary material can be found at http://www.cz.oxfordjournals.org/. Ferchaud and Hansen 2016). The fixation of standing genetic vari- ants in highland populations (HB or HA) was inferred for most alleles References (8/10, except for a49/108 in T2) (Zhu et al., unpublished data) that might be homogenized by neutral gene flow. Similarly, much nar- Bandelt H-J, Forster P, Ro ¨ hl A, 1999. Median–joining networks for inferring rower shifts but unfixed highland a49Asn/108Ala alleles in HA and intraspecific phylogenies. Mol Biol Evol 16:37–48. Bulgarella M, Peters JL, Kopuchian C, Valqui T, Wilson RE et al., 2012. wider cline shifts (a22/35/72 of T2) may have been caused by neutral Multilocus coalescent analysis of haemoglobin differentiation between low– gene flow in clade A. In agreement with these results, the a22/35 were and high–altitude populations of crested ducks Lophonetta specularioides. absent from altitudinal clines in the full sampling, whereas only the Mol Ecol 21:350–368. completely diverged a72Asn/His suggested that an ancestral standing Cheviron ZA, Natarajan C, Projecto-Garcia J, Eddy DK, Jones J et al., 2014. polymorphism was fixed in all highland populations of the great tit Integrating evolutionary and functional tests of adaptive hypotheses: a case complex. Therefore, although the same elevational differentiation was study of altitudinal differentiation in hemoglobin function in an Andean detected in the genetic structure of the superspecies, specific allele pro- sparrow Zonotrichia capensis. Mol Biol Evol 31:2948–2962. files of in P. minor populations were affected by local homogenizing Colosimo PF, Hosemann KE, Balabhadra S, Villarreal G, Dickson M et al., gene flow. In a word, differed from the case of standing Hb variants 2005. Widespread parallel evolution in sticklebacks by repeated fixation of originated in low-altitude Andean waterfowl (Natarajan et al. 2015), ectodysplasin alleles. Science 307:1928–1933. Excoffier L, Lischer HEL, 2010. Arlequin suite ver 3.5: a new series of pro- the recurrently replaced alleles in deeply diverged lineages of Eastern grams to perform population genetics analyses under Linux and Windows. Himalaya and Central Asia can be interpreted as the fixation of stand- Mol Ecol Resour 10:564–567. ing high-altitude alleles in East-Himalaya. Ferchaud AL, Hansen MM, 2016. The impact of selection, gene flow and In summary, our study highlights the necessity of integrating demographic history on heterogeneous genomic divergence: three-spine gene flow and demographic history into exploring potential adap- sticklebacks in divergent environments. Mol Ecol 25:238–259. tive evolution in the genetic divergence of widespread species. This Fu YX, 1997. Statistical tests of neutrality of mutations against population preliminary investigation on genetic variants in a widespread passer- growth, hitchhiking and background selection. Genetics 147:915–925. ine superspecies provides an ideal model to further explore the adap- Fu YX, Li WH, 1993. Statistical tests of neutrality of mutations. Genetics 133: tive evolution of Hb genes in wild birds. 693–709. Galen SC, Natarajan C, Moriyama H, Weber RE, Fago A et al., 2015. Contribution of a mutational hot spot to hemoglobin adaptation in high- altitude Andean house wrens. Proc Natl Acad Sci USA 112: 13958–13963. Gill F, Donsker D, 2016. Ioc world bird list (v 6.4) [Cited 2016 May 21] Data archiving Available from: http://www.worldbirdnames.org/ All DNA sequences are accessible under GenBank accession num- Gosler A, Clement P, Christie DA, 2017. Great tit Parus major. In: del Hoyo J, bers (MF321902–MF322504), and sample information for all indi- Elliott A, Sargatal J, Christie DA, Juana d, editors. Hand Book of the Birds viduals is found in Table S1 of the Supplementary file. of the World Alive. Barcelona: Lynx Edicions. Hey J, Nielsen R, 2004. Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. Persimilis. Genetics 167: 747–760. Acknowledgments Hey J, Nielsen R, 2007. Integration within the Felsenstein equation for im- proved Markov chain Monte Carlo methods in population genetics. Proc We thank Shimiao Shao for providing comments on the early draft, Yalin Natl Acad Sci USA 104:2785–2790. Chen and Dezhi Zhang for revising the latest version of this manuscript. We Kumar A, Natarajan C, Moriyama H, Witt CC, Weber RE et al., 2017. also thank Per Alstro ¨ m for valuable comments on this work. We thank Stability-mediated epistasis restricts accessible mutational pathways in the Tianlong Cai for providing the elevation maps and Ruiying Zhang for provid- functional evolution of avian hemoglobin. Mol Biol Evol 34:1240–1251. ing sample information. Kvist L, Arbabi T, P€ ackert M, Orell M, Martens J, 2007. Population differenti- ation in the marginal populations of the great tit (Paridae: Parus major). Biol J Linn Soc 90:201–210. Kvist L, Martens J, Higuchi H, Nazarenko AA, Valchuk OP et al., 2003. Funding Evolution and genetic structure of the great tit Parus major complex. Proc R This work was funded by the Strategic Priority Research Programme of the Soc Lond B Biol Sci 270:1447–1454. Chinese Academy of Sciences (XDB13020300), State Key Programme of Librado P, Rozas J, 2009. Dnasp v5: a software for comprehensive analysis of Natural Science Foundation of China (NSFC 31330073), a grant from the DNA polymorphism data. Bioinformatics 25:1451–1452. Ministry of Science and Technology of China (2014FY210200) to F.L., and a Mairbaurl H, Weber RE, 2012. Oxygen transport by hemoglobin. Compr grant from the NSFC programme (J1210002) to Y.G. and X.J.Z. G.S. is sup- Physiol 2:1463–1489. ported by the NSFC fund (31471991) and a grant (Y229YX5105) from the McCracken KG, Barger CP, Bulgarella M, Johnson KP, Kuhner MK et al., Key Laboratory of the Zoological Systematics and Evolution of the Chinese 2009a. Signatures of high-altitude adaptation in the major hemoglobin of Academy of Sciences. five species of Andean dabbling ducks. Am Nat 174:631–650. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018 10 Current Zoology, 2017, Vol. 0, No. 0 McCracken KG, Bulgarella M, Johnson KP, Kuhner MK, Trucco J et al., Qu Y, Tian S, Han N, Zhao H, Gao B et al., 2015. Genetic responses to sea- 2009b. Gene flow in the face of countervailing selection: adaptation to high- sonal variation in altitudinal stress: whole-genome resequencing of great tit altitude hypoxia in the bA hemoglobin subunit of yellow-billed pintails in in eastern Himalayas. Sci Rep 5:14256. the Andes. Mol Biol Evol 26: 815–827. Raeymaekers JAM, Konijnendijk N, Larmuseau MHD, Hellemans B, De Monge C, Leon-Velarde F, 1991. Physiological adaptation to high altitude: Meester L et al., 2014. A gene with major phenotypic effects as a target for oxygen transport in mammals and birds. Physiol Rev 71: 1135–1172. selection vs. homogenizing gene flow. Mol Ecol 23:162–181. Munoz-Fuentes V, Cortazar-Chinarro M, Lozano-Jaramillo M, McCracken Schluter D, Clifford EA, Nemethy M, McKinnon JS, 2004. Parallel evolution KG, 2013. Stepwise colonization of the Andes by ruddy ducks and the evo- and inheritance of quantitative traits. Am Nat 163:809–822. lution of novel beta-globin variants. Mol Ecol 22: 1231–1249. Schluter D, Conte GL, 2009. Genetics and ecological speciation. Proc Natl Natarajan C, Projecto-Garcia J, Moriyama H, Weber RE, Munoz-Fuentes V Acad Sci USA 106 (1 Suppl): 9955–9962. et al., 2015. Convergent evolution of hemoglobin function in high-altitude Scott GR, Milsom WK, 2006. Flying high: a theoretical analysis of the factors Andean waterfowl involves limited parallelism at the molecular sequence limiting exercise performance in birds at altitude. Respir Physiol Neurobiol level. PLoS Genet 11:e1005681. 154:284–301. Olson-Manning CF, Wagner MR, Mitchell-Olds T, 2012. Adaptive evolution: Stephens M, Smith NJ, Donnelly P, 2001. A new statistical method for haplo- Evaluating empirical support for theoretical predictions. Nat Rev Genet type reconstruction from population data. Am J Hum Genet 68:978–989. 13:867–877. Storz JF, 2007. Hemoglobin function and physiological adaptation to hypoxia Opazo JC, Hoffmann FG, Natarajan C, Witt CC, Berenbrink M et al., 2015. in high-altitude mammals. J Mammal 88:8. Gene turnover in the avian globin gene families and evolutionary changes in Strasburg JL, Rieseberg LH, 2010. How robust are ‘isolation with migration’ hemoglobin isoform expression. Mol Biol Evol 32:871–887. analyses to violations of the IM model? A simulation study. Mol Biol Evol P€ ackert M, Martens J, Eck S, Nazarenko AA, Valchuk OP et al., 2005. The great 27:297–310. tit Parus major: a misclassified ring species. Biol J Linn Soc 86:153–174. Tietze DT, Borthakur U, 2012. Historical biogeography of tits (Aves: Paridae, Perutz MF, 1983. Species adaptation in a protein molecule. Mol Biol Evol Remizidae). Org Divers Evol (12):433–444. 1:1–28. Weber RE, 2007. High-ltitude adaptations in vertebrate hemoglobins. Respir Porter AH, Wenger R, Geiger H, Scholl A, Shapiro AM, 1997. The Pontiac Physiol Neurobiol 158: 132–142. daplidice - edusa hybrid zone in northwestern Italy. Evolution Wilson RE, Peters JL, McCracken KG, 2012. Genetic and phenotypic diver- 51:1561–1573. gence between low– and high–altitude populations of two recently diverged Pritchard JK, Stephens M, Donnelly P, 2000. Inference of population structure cinnamon teal subspecies. Evolution 67: 170–184. using multilocus genotype data. Genetics 155:945–959. Zhao N, Dai C, Wang W, Zhang R, Qu Y et al., 2012. Pleistocene climate Projecto-Garcia J, Natarajan C, Moriyama H, Weber RE, Fago A et al., 2013. changes shaped the divergence and demography of Asian populations of the Repeated elevational transitions in hemoglobin function during the evolution great tit Parus major: evidence from phylogeographic analysis and ecolo- of Andean hummingbirds. Proc Natl Acad Sci U S A (51):20669–20674. gical niche models. J Avian Biol 43: 297–310. Downloaded from https://academic.oup.com/cz/advance-article-abstract/doi/10.1093/cz/zox042/3929961 by guest on 13 July 2018

Journal

Current ZoologyOxford University Press

Published: Jul 6, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off