Parallel genomic responses to historical climate change and high elevation in East Asian songbirds

Extreme environments present profound physiological stress. The adaptation of closely related species to these environments is likely to invoke congruent genetic responses resulting in similar physiological and/or morphological adaptations, a process termed “parallel evolution” (1). Existing evidence shows that parallel evolution is more common at the phenotypic level than at the genic level (2, 3). Despite some clear examples of parallelism at the same amino acid substitution (4, 5), parallel phenotypes are often generated by nonparallel genetic changes (6, 7), especially for polygenic traits (8). Given the universality of gene-level nonparallelism, recent studies demonstrate that molecular evolution is more likely to be parallel at levels above the gene [i.e., genomic parallelism (9, 10)]. However, similar to genic parallelism, parallel genomic evolution also has inconsistent outcomes, ranging from the absence of any parallelism (11) and similarity in biological functions or pathways (9) to genome-wide abundant parallelism (12). Theoretical studies suggest that these inconsistent adaptive solutions may be influenced by differences in standing genetic variation (i.e., heterozygosity) (13, 14). A loss of heterozygosity within small populations might reduce opportunities for parallel evolution across replicated species by decreasing the probability that the same alleles are available in response to shared selection pressures (1517).

Levels of heterozygosity are influenced by historical demography. The genetic legacy of shifts in species’ ranges and population size due to Quaternary glacial fluctuations is well documented (18), and many terrestrial vertebrates from glaciated regions experienced historical population bottlenecks as a result of glaciation (19, 20), including birds (21). However, organisms currently inhabiting previously glaciated regions may have different colonization histories and therefore different demographic histories, resulting in different levels of heterozygosity among cooccurring extant species. Species complexes consisting of species with variable demographic histories inhabiting extreme environments provide an opportunity to test the hypothesis that levels of heterozygosity mediate the extent to which adaptive genomic responses occur in parallel. We put this hypothesis to a test through comparing genomes of closely related, codistributed tit species across East Asia. As western populations of widespread species may have relied on higher standing variation in comparison with western endemics, we test whether the former show also a higher degree of genic parallelisms in adaptation to high altitude (HA).

East Asia is a region of extreme environmental variation. Most of the east is at or near sea level, but western East Asia is dominated by the Qinghai–Tibet Plateau (QTP), with an average elevation of 4,500 m above sea level (m a.s.l.) (Fig. 1A). During the Pleistocene, this western region was covered by glaciers, whereas the eastern lowlands were largely ice free (22, 23). Songbirds found in the western highlands face several extreme environmental conditions, including low availability of oxygen and colder annual temperatures that exert strong selective pressures (7, 24). This spatial environmental heterogeneity over geologic time provides a rare opportunity to investigate the similarity of genomic responses to the combination of historical climate change and high-elevation extreme environments.

The tits and chickadees (Passeriformes: Paridae) are an avian lineage with limited dispersal ability that originated at middle elevations (∼2,200 m a.s.l.) (SI Appendix, Fig. S1) in the Sino-Himalayan region (25), subsequently expanding into both the eastern lowlands and western highlands (up to 5,500 m a.s.l.) throughout East Asia. Phylogeographic studies (based on mitochondrial DNA) of wide-ranging East Asian tit species that span east to west in East Asia show a consistent pattern of lack of gene flow and east–west population differentiation that dates to the Early- to Mid-Pleistocene or earlier, likely driven by both ecological differences between west and east as well as physical barriers to gene flow (2628). Moreover, western high-elevation parids show evidence for adaptation to life at extreme elevations (7, 29). The QTP-native Pseudopodoces humilis and the widespread Parus major evolved distinct adaptive responses for increased energy metabolism to cope with high elevation (30, 31). Likewise, although multiple high-elevation parids show parallel evolution for increased blood–oxygen affinity, these are due to unique divergent substitutions in hemoglobin genes (7). However, a comparative transcriptomic analysis between two high-elevation parids (Lophophanes dichrous and Periparus rubidiventris) detected similarity in genes associated with high-elevation adaptations at the gene expression level (32). Here, we analyzed genome data from all 19 extant East Asian tit species to identify patterns of genetic diversity in response to Pleistocene glaciations and to test for the degree of parallelism in closely related species’ responses to extreme high-elevation environments, as well as the effects of genetic diversity on adaptive parallelism.

Results

Sequence Analysis and Time Trees.

We sequenced whole genomes from 19 parid species (n = 87) and two species from the closely related Remizidae (n = 4; Fig. 1A) to a mean genome-wide depth of coverage of 18.7× (range:10.2 to 47.8×) per individual. On average, 95% (range: 71 to 98%) of the reads aligned to the chromosome-assembled P. major reference genome (33) (SI Appendix, Fig. S2) and resulted in a dataset of 43,458,033 single-nucleotide polymorphisms (SNPs; SI Appendix, Table S1).

Phylogenies constructed from this SNP dataset consistently recovered a single identical topology with fully supported nodes for all 19 parids species, irrespective of tree-building methods (both genome-wide concatenation and coalescent-based species tree approaches) or whether we used autosomal SNPs or only Z-linked SNPs (SI Appendix, Fig. S3). We converted this phylogeny to a time tree using coding genes from the de novo assemblies as follows: For each species (or population, for widespread pecies that span both western and eastern East Asia), we generated a “rough” de novo genome assembly from which we recovered orthologous protein-coding genes. After translating nucleotide sequences to amino acid sequences, we identified a set of 1,885 one-to-one orthologs among the 24 assemblies (SI Appendix, Table S2). Of these orthologs, we found that 246 had evolved in a clock-like manner. We then concatenated these 246 genes, constrained the above topology, and enforced a rate-calibrated molecular clock. The resulting time tree showed that Paridae diverged from Remizidae ∼26 Mya, with subsequent splits between extant eastern and western taxa spanning a range from ∼13 to 1.3 Mya (Fig. 1B). Western and eastern taxa did not form independent clades; instead, they were distributed across the tree. Likewise, our principal component analysis and multidimensional scaling analysis rejected the hypothesis that western and eastern genomes formed distinct clusters (SI Appendix, Fig. S4 A and B). Nodes connecting eastern and western species dated to earlier than 3 Mya (i.e., prior to the Pleistocene), except for intraspecific nodes within three widespread species (P. major, Poecile montanus, and Periparus ater) that exhibited shallower divergence times, as expected for within-species comparisons. Here, estimated intraspecific divergence times between eastern and western populations ranged from 1.8 to 1.3 Mya (Fig. 1B). These results more or less agree with previous phylogenetic dating efforts for Paridae (25, 34, 35) and confirm that the most recent intraspecific splits between western and eastern populations date to the Early- to Mid-Pleistocene.

Patterns of Genetic Diversity.

East Asian parids had relatively low levels of genome-wide nucleotide diversity (π) (average: 0.74 × 10−3; SI Appendix, Table S1). Widespread P. major had the highest π (1.77 × 10−3), whereas the Taiwan island–endemic Machlolophus holsti had the lowest π (0.07 × 10−3), consistent with expectations from prolonged insularization. We found an extremely high correlation between π and heterozygosity (H) (R2 = 0.97, P = 6.76 × 10−16; SI Appendix, Fig. S5A), so we used the latter to analyze landscape-level patterns of genomic diversity, as H is less likely to be influenced by sample size and uneven geographic sampling. Kriging-smoothed estimates of phylogroup-wide H were significantly lower in all western birds than in all eastern birds (Wilcoxon rank-sum = 634, P = 0.008; Fig. 2A), indicating congruent and species-independent patterns of genomic diversity across all East Asian parids. Specifically, species endemic to the western highlands had consistently lower H than species found exclusively in the eastern lowlands (Wilcoxon rank-sum = 336.5, P = 0.009; Fig. 2B). However, the three widespread species (P. major, Pe. ater, and Po. montanus) that ranged from the eastern lowlands to the western highlands had a significantly higher H than western endemic species (Wilcoxon rank-sum = 565, P = 0.001) but similar to that observed in eastern endemic species (Wilcoxon rank-sum = 343.5, P = 0.12), although their western populations had a significantly lower H than conspecific eastern populations (Fig. 2B). The same trend of reduced H among western endemic species was observed in both coding and noncoding regions, and in fact, among western endemics H was significantly lower in coding regions relative to the remaining genome (i.e., noncoding regions; Wilcoxon rank-sum = 7,225, P < 2.2 × 10−16) (SI Appendix, Fig. S5B).

Fig. 2.
Fig. 2.

Landscape patterns of lineage-wide genetic diversity and historical demography. (A) Spatial distribution of genomic diversity based on the heterozygosity (H) of each individual. Western parids have significantly lower H than eastern parids (**P < 0.01). The white contour lines indicate the extent of glaciers during the LGM (0.023 to 0.018 Mya). The map of glaciers is from www.deviantart.com/atlas-v7x. (B) Comparisons of H among western endemic, eastern endemic, and widespread species. The H of widespread species is significantly higher than that of western endemics but similar to that of eastern endemics. Within widespread species, western populations have significantly lower H than all eastern populations; *P < 0.05; NS, nonsignificant. (C) Historical demography of western and eastern endemic parids inferred by PSMC. All nine western endemic species (red curves) experiences a contraction of effective population size (i.e., Ne contraction) during the LGP (0.11 to 0.012 Mya): six (d, f, g, h, k, l) declines from a high point at or near the beginning of the LGP, while in three (b, c, j), the decline starts at or near the LGM. Five (2, 3, 5, 9, 10) of six eastern endemic species (black curves) experiences increase in Ne during the LGP, the sixth (7) experiences a decline. (D) Historical demography of the three widespread species. The eastern population of each species experiences demographic expansion during the LGP, yet two western populations (Pe. ater and Po. montanus) also experience expansion (in contrast, western P. major experienced a contraction). The lowercase letters and numbers above the lines correspond to the different species (Fig. 1B) whose localities are labeled on the map in A. A total of 100 bootstraps for the demography of each species are shown in SI Appendix, Fig. S6. The legend shows changes in temperature during the past 120,000 y and indicates the intervals of the LGP and LGM glacial periods. ΔTs is the global mean surface temperature anomaly from Zachos et al. (88) and Hansen et al. (89) plotted by year.

Patterns of Historical Demography.

Pairwise sequentially Markovian coalescent (PSMC) estimates of effective population size (Ne) over time showed markedly different patterns in the demographic history of western and eastern species during the Pleistocene. All nine western endemic species that occurred only in western East Asia showed a significant decline in Ne during this period. For six of these, the decline began at or near the beginning of the Last Glacial Period (LGP, 0.11 Mya); for the other three species, demographic contraction began at the Last Glacial Maximum (LGM, 0.023 to 0.018 Mya) (Fig. 2C and SI Appendix, Fig. S6A and Table S3). In contrast, five of six eastern endemic species showed significant demographic expansion during the Pleistocene. Here, three showed expansions starting at or near the LGM, while for the remaining two expansion started earlier during the LGP (Fig. 2C and SI Appendix, Fig. S6B). In all three widespread species (P. major, Po. montanus, and Pe. ater), eastern populations showed significant population expansion during the LGP, while two of the three western populations showed significant expansion during the LGP (Fig. 2D and SI Appendix, Fig. S6C), rather than the contraction seen in western endemic species (Fig. 2C). This finding, along with the lower H observed in western relative to eastern populations of each widespread species (Fig. 2B), suggests that all three species underwent population expansion from lowland areas into the western region during the Pleistocene, and not the reverse, which would also agree with the split dates obtained in our time tree (Fig. 1B). Finally, the Taiwan-endemic species (Ma. holsti) showed significant and constant population expansion throughout the LGP (SI Appendix, Fig. S6D), which was probably consistent with population recovery after a founder event (36).

Genome-Wide Patterns of Elevation-Associated Differentiation.

At our sampling locations ranging from sea level to 4,550 m a.s.l., available oxygen (partial pressure of oxygen: PO2) declines sharply with increased elevation from 21 to 12% (R2 = 0.99, P < 2.2 × 10−16; SI Appendix, Fig. S8A), representing an almost 50% reduction in inspired O2. The annual mean temperature also declined with elevation (R2 = 0.45, P = 9.85 × 10−13; SI Appendix, Fig. S8A), but the relationship was more variable. For small-bodied endotherms such as tits, the physiological stress of hypoxia and cold tolerance at extreme elevation are inexorably linked because these animals cannot reduce metabolic activity to mitigate against reduced oxygen availability given their thermogenic requirements (7, 37, 38).

Following prior work (39), we defined extremely high elevations as above 3,000 m a.s.l. Therefore, we classified 37 out of 45 western birds sampled at or above 3,000 m a.s.l. as HA parids and classified 38 out of 42 eastern birds sampled at or below 1,000 m a.s.l. as low-altitude (LA) parids (SI Appendix, Fig. S8A and Dataset S1). All HA parids are from western populations, and all LA parids are from eastern populations; the three widespread species have both HA and LA populations. We compared allele frequency (AF) distributions between all HA and all LA individuals as a first-step filter to detect shared alleles that might be responsible for parallel adaptation to high elevation. We found that none of the more than 43 million SNPs in our dataset (see Sequence Analysis and Time Trees) were differentially fixed (SI Appendix, Fig. S8B), and more generally, we observed low levels of the fixation index (FST) between HA and LA parids: only 5,589 (0.013%) SNPs had FST values above 0.6 (SI Appendix, Fig. S8C). Low-fixation rates were also observed between HA endemics and their LA relatives as well as between HA populations and LA populations of widespread species (SI Appendix, Fig. S8C). However, in each one-to-one HA–LA comparison, we consistently found a high proportion of fixed SNPs (average: 46.49%, range: 18.19 to 70.12%) (SI Appendix, Fig. S8D) due to either abundant adaptive substitutions unique to each species or the potential effect of small sample sizes on FST calculation.

To further evaluate whether high-elevation adaptation among East Asian parids was the result of shared genomic processes, we calculated FST for each comparison in 50-kb sliding windows, which allowed us to identify highly differentiated genomic regions (i.e., the genomic windows in the top 1% percentile of FST values) shared among HA–LA comparisons (Fig. 3A and SI Appendix, Fig. S9A). Among these HA–LA pairs, we found a low proportion (average: 11.1%, range: 0.5 to 31%) of outlier windows shared by any two pairs (SI Appendix, Fig. S9B). All four HA endemic taxa had only two windows in common, while the three HA–LA comparisons involving widespread taxa had none (SI Appendix, Fig. S9C). While some HA taxa had a single obvious LA comparative taxon/population, in other cases, two or three LA taxa were equally closely related to a given HA taxon (Fig. 1B and SI Appendix, Fig. S1). Therefore, to avoid any bias introduced by multiple comparisons with the same LA taxa, we chose one HA–LA paired comparison at random from each HA endemic, plus three HA–LA population comparisons of widespread species, a total of seven HA–LA comparisons, for the following analyses (Fig. 3A). At the conclusion of the study, we repeated the analyses with all possible relevant HA–LA comparisons and found consistent results. Finally, we observed that within a single HA-endemic species and its multiple comparisons, the number of shared windows (average: 31%, range: 12 to 48%; SI Appendix, Fig. S9D) was significantly higher than that across HA–LA comparisons (Wilcoxon rank-sum = 116, P = 0.001).

Fig. 3.
Fig. 3.

FST outlier windows are enriched with HA-related genes. (A) Window-based comparisons of each HA taxon with its LA relative. Red dashed lines indicate cutoff values of the ZFST distribution in autosomes and the Z chromosome. A total of 202 highly divergent windows are identified for each HA species. Illustrations of parids are from ref. 73. (B) Number of potentially HA-related genes found in FST outlier windows by GO term. The 17 HA-related GO terms are grouped by the four physiological stages of the O2 transport cascade and/or thermogenesis plus a fifth “hypoxia” category. The unfilled bar represents the log-transformed P value for the enrichment of GO terms, and the color-filled bar represents the observed number of genes associated with that GO term found in outlier FST windows. Enrichment is a measure of the overrepresentation of particular genes by GO terms. Phu, Ps. humilis; Psu, Po. superciliosus; Ldi, Lo. dichrous; Pru, Pe. rubidiventris; Pat, Pe. ater; Pma, P. major; Pmo, Po. montanus.

To investigate the role of selection in causing high FST in our HA–LA comparisons, we examined whether outlier windows were in regions of exceptionally high absolute divergence (dXY “islands”) or exceptionally low absolute divergence (dXY “valleys”). Among the four HA–LA comparisons involving HA endemics, all sets of FST outlier windows had significantly higher dXY values than nonoutlier windows (SI Appendix, Fig. S9E) and had significantly more overlap with dXY “islands” than with dXY “valleys” [two-sample t test, t(3.66) = –3.84, P = 0.011; SI Appendix, Fig. S9F]. In contrast, among the three widespread taxa comparisons, significantly lower dXY was observed in two of the three FST outlier sets compared with nonoutliers (SI Appendix, Fig. S9E), and all three sets showed more overlap with dXY “valleys” than with dXY “islands” but without significance [two-sample t test, t(4) = 1.73, P = 0.08; SI Appendix, Fig. S9F]. This result appears to be consistent with the pattern of more recent adaptive divergence in widespread taxa and more long-term adaptive divergence in HA endemics because recent linked selection may result in dXY “valleys,” whereas long-term extreme selection that can fix a large number of substitutions may cause dXY “islands” (40).

Polymorphic Genes in Highly Differentiated Windows Are Disproportionally Assigned to Gene Ontology Terms Associated with High-Elevation Adaptation.

Elevated FST may arise from both selective and stochastic processes. To address to what extent differentiation may be the consequence of high-elevation adaptation, we tested for enrichment of HA-related gene functions using Gene Ontology (GO) enrichment analyses. Indeed, HA-related GO terms (averageobs: 31.14, range: 24 to 40) were significantly enriched in outlier windows (SI Appendix, Fig. S10A), compared with 100 randomly selected sets of genes (averageexp: 2.64, range: 0 to 9) (Wilcoxon rank-sum = 700, P = 4.06 × 10−6; SI Appendix, Fig. S10B). To identify all genes possibly associated with adaptation to high elevation from the outlier windows, we screened all coding sequences against a database of 17 parent-level GO terms relevant to the oxygen transport cascade and/or thermogenesis and therefore could be related to high-elevation adaptation (Materials and Methods) (Fig. 3B). Including dependent (i.e., “children”) terms in this analysis, we obtained a list of 973 HA-related GO terms that we used for all subsequent GO analyses (Dataset S2). This set included a total of 230 HA-related genes from outlier windows of all seven comparisons. These genes formed a network connecting many significantly enriched pathways, some of which were involved in physiological processes that may be associated with high-elevation adaptation, such as the following: “HIF-1 signaling pathway,” “vascular smooth muscle contraction,” “fatty acid biosynthesis,” “phosphatidylinositol signaling system,” “glucagon signaling pathway,” and “cAMP signaling pathway” (Fig. 4A and Dataset S3). Across the seven comparisons, a total of 212 pathways were found showing high overlap; 19 occurred in all comparisons, 52 occurred in four comparisons of HA endemics, and 30 occurred in three comparisons of widespread species (Fig. 4B).

Fig. 4.
Fig. 4.

East Asian parids show parallel adaptive responses to high elevation at the level of biological function. (A) A Kyoto Encyclopedia of Genes and Genomes (KEGG) network for 230 enriched HA-related genes from all outlier windows of the seven comparisons. These pathways, especially the “HIF-1 signaling pathway,” “vascular smooth muscle contraction,” “fatty acid biosynthesis,” and “phosphatidylinositol signaling system,” may be associated with high-elevation adaptation. The size and color of circles indicate P values and different groups, while the line linking any two pathways shows the connectivity scored by a kappa score (see details in Dataset S3). (B) High similarity in KEGG pathways observed across comparisons. Nineteen pathways cooccur in all seven comparisons, 52 pathways cooccur in the four comparisons of HA endemics, and 30 pathways cooccur in the three comparisons of widespread species. (C) High FST outlier windows have a lower number of total genes compared with random windows. (D) However, outlier windows have significantly more HA-related genes than random windows. (E) Additionally, outlier windows have significantly more enriched HA-related GO terms than random windows. Cell color in CE indicates the scaled z-score (by SD) within each HA–LA comparison. (F) FST outlier windows have similar overlap of enriched HA-related GO terms between HA endemics and widespread taxa, which is significantly higher than random windows. Cell color indicates the log-transformed P value generated by the hypergeometric test. (G) Likewise, enriched HA-related GO terms have significantly higher semantic similarity in FST outlier windows than in random windows, but no difference is observed between HA endemics and widespread taxa. Cell color indicates the similarity degree. In both F and G, outlier windows and random windows (solid boxes) as well as HA endemics and widespread taxa (dashed boxes) are compared; ***P < 0.001; NS, nonsignificant. Phu, Ps. humilis; Psu, Po. superciliosus; Ldi, Lo. dichrous; Pru, Pe. rubidiventris; Pat, Pe. ater; Pma, P. major; Pmo, Po. montanus.

We then asked whether outlier windows had a higher frequency of HA-related genes than randomly sampled genomic regions. Indeed, FST outlier windows had significantly more HA-related genes (averageobs: 51.9, range: 40 to 66) than random windows (averageexp: 17, range: 9 to 26) [two-sample t test, t(6) = 11.58, P = 7.20 × 10−8; Fig. 4D], although overall, these outlier windows had, on average, fewer genes (averageobs: 230.7, range: 197 to 291) than random windows (averageexp: 276.6, range: 230 to 324) [two-sample t test, t(6) = –3.37, P = 0.006; Fig. 4C]. Similarly, outlier windows harbored significantly more HA-related genes than random windows when we controlled for an equal number of genes (Wilcoxon rank-sum = 0, P = 9.93 × 10−6; SI Appendix, Fig. S10C). Furthermore, there were significantly more GO terms enriched by HA-related genes in FST outlier windows (averageobs: 335.9, range: 218 to 583) than in random windows (averageexp: 78.7, range: 0 to 568) [two-sample t test, t(6) = 5.03, P = 0.0003; Fig. 4E]. These results are consistent with the notion that the high differentiation observed in FST outlier genomic regions is at least partially the result of adaptation to extreme elevation in East Asian parids.

Partial Overlap of HA-Related Genes.

The number of HA-related genes found in outlier windows varied little among seven pairs (average: 51.9, range: 40 to 66) (χ2 = 7.27, df = 6, P = 0.30), but we recovered only partial overlap among HA-related genes across the seven pairs. On average, any two HA–LA parid pairs shared 10.9 HA-related genes (range: 0 to 22; SI Appendix, Fig. S11A), with only one gene (PLB1) found in six of the seven sets of FST outlier windows (SI Appendix, Fig. S11B), and no gene was found in all seven pairwise comparisons (SI Appendix, Fig. S11C). Pairs of HA endemic species shared slightly more HA-related genes than did pairs of widespread species (nendemic = 14.8 versus nwidespread = 6.7), but the difference was not significant [two sample t test, t(7) = 1.79, P = 0.12], according with the similarity in the number (nendemic = 52.5 versus nwidespread = 51) of HA-related genes between them [two-sample t test, t(5) = 0.54, P = 0.62]. Moreover, we found that four HA-related genes (PLB1, MORC2, RBKS, and GLP1R) occurred in all outlier window sets of the four HA endemics, but only one HA-related gene (NDST2) occurred in outlier windows for all three widespread species (SI Appendix, Fig. S11D). PLB1 (41) and MORC2 (42) are associated with lipid metabolic processes; RBKS (43) and NDST2 (44) take part in carbohydrate metabolism, whereas GLP1R regulates blood pressure (45). These results may suggest only limited parallelism among particular genes potentially involved in high-elevation adaptation.

Substantial Overlap of HA-Related GO Terms.

In outlier windows, HA-related GO terms had significantly higher pairwise overlap [two-sample t test: t(40) = 5.95, P < 5.55 × 10−7] and semantic similarity [two-sample t test: t(40) = 5.15, P < 7.41 × 10−6] than HA-unrelated GO terms. Likewise, the number (averageobs: 122, range: 51 to 212) of HA-related GO terms cooccurring between any two HA–LA comparisons was significantly higher than expected by chance (averageexp: 7, range: 0 to 182) (Wilcoxon rank-sum = 806,722, P < 2.2 × 10−16; Fig. 4F), and the set of HA-related GO terms cooccurring between any two HA–LA comparisons also had significantly greater semantic similarity (averageobs: 0.7, range: 0.58 to 0.77) than expected by chance (averageexp: 0.39, range: 0 to 0.81) (Wilcoxon rank-sum = 636,880, P < 2.2 × 10−16; Fig. 4G). In contrast, enriched HA-unrelated GO terms were significantly less likely to cooccur in outlier windows (averageobs: 54, range: 32 to 77) than expected by chance (averageexp: 110, range: 25 to 288) (Wilcoxon rank-sum = 62,648, P < 2.2 × 10−16; SI Appendix, Fig. S12A), and between pairs of HA–LA comparisons, sets of cooccurring GO terms unrelated to HA also had significantly less semantic similarity in outlier windows (averageobs: 0.55, range: 0.25 to 0.76) than did sets from random windows (averageexp: 0.65, range: 0.35 to 0.85) (Wilcoxon rank-sum = 199,654, P = 1.23 × 10−8; SI Appendix, Fig. S12B). Among pairs of both HA endemics and widespread taxa, we consistently observed high but not significantly different overlap in HA-related GO terms [two-sample t test: t(7) = 1.10, P = 0.31] and low, but not significantly different, overlap in HA-unrelated GO terms [two-sample t test: t(7) = –0.44, P = 0.68]. Likewise, high semantic similarity in HA-related GO terms [two-sample t test: t(7) = 1.31, P = 0.23] and low semantic similarity in HA-unrelated GO terms [two-sample t test: t(7) = 0.55, P = 0.60] were found equally in both sets of pairs. Thus, while individual genes showed little overlap among pairwise comparisons, GO terms related to extreme elevation adaptation showed substantial overlap in these same comparisons.

Biological Functions and Processes Associated with Genomic Outlier Windows.

The 17 GO terms (and their dependents) varied in the frequency in which they were enriched among our seven HA–LA comparisons (Fig. 3B and SI Appendix, Table S4). Terms associated with “respiration,” “circulation,” and “oxygen diffusion” in tissues were less often enriched (on average, between 1 and 2.8 out of 7 comparisons), although we noted that within the circulation terms, two terms more directly associated with circulatory system development (i.e., “vascularization”) were enriched in six of seven comparisons. The three most frequently enriched GO terms were all associated with utilization of oxygen in tissues and cells (GO:0046034, GO:0006629, and GO:0005975; on average five out of seven comparisons), as was the hypoxia-specific term (GO:0036293, six of seven comparisons). Importantly, GO terms for both “lipid metabolic processes” (GO:0006629) and “carbohydrate metabolic processes” (GO:0005975) were enriched in all seven HA–LA comparisons; these processes serve dual functions related to oxygen utilization and thermogenesis.

Evidence for Selection in HA-Related Genes.

We found evidence for positive selection on amino acid changes in just one HA–LA comparison and only in 1 gene out of 230 HA-related genes (Dataset S4). This gene, MORC2, is involved in lipid metabolism (42), a process that has been shown previously to be important for high-elevation adaptation in several vertebrates (30, 46). We note that the codon-based selection tests used here (47) have less power and are not able to identify adaptive evolution in regulatory and other noncoding regions unlike selective sweep tests, which require population-level sampling beyond the scope of this project. Nonetheless, our observation of overrepresentation of genes with functions related to high-elevation adaptation in outlier windows, but limited evidence for positively selected amino acid substitutions in those genes appears consistent with hypotheses that emphasize selection on gene expression and regulatory processes rather than on novel amino acids in explaining vertebrate adaptation to high elevation (32, 46, 4850). However, additional studies at population levels are needed to disentangle which genomic regions selection truly acts on.

Discussion

Shared History Results in Parallel Patterns of Genome-Wide Diversity.

The genomes of East Asian parids have been shaped by both history and environment. The demographic legacy of Pleistocene climate change had a parallel effect on genome-wide patterns of diversity. All nine western endemic species experienced a significant reduction in effective population size during the LGP, potentially consistent with glacial expansion in western East Asia. Glaciation of the western highlands would have made most of the region uninhabitable, with endemic parids restricted to isolated refugia relative to their current distributions (51), as demonstrated by declining coalescent-based population size estimates from the Pleistocene toward the present for all nine species of western-restricted East Asian parids (Fig. 2C). In contrast, eastern parid taxa, which are exclusively lowland, do not show a similar reduction in genome-wide diversity. The eastern lowlands of East Asia were unglaciated during the LGP (Fig. 2A), resulting in relatively stable environments throughout the Quaternary (22, 23). We found no evidence for Pleistocene population size contraction in eastern parids, in agreement with previous studies (26, 27, 52, 53). In fact, most eastern populations showed evidence of growth since the LGM, which is different from previously hypothesized post-LGM expansion scenarios based solely on mitochondrial DNA (mtDNA) population variation (27, 51) and is consistent with a scenario of population growth coincident with poleward expansion of a suitable habitat for lowland eastern parids. Previous studies (31, 54) have hypothesized that reduced genetic diversity in high-elevation animal populations may be a direct consequence of adaptation to life at HA and/or the founder events that permitted these populations to become established in the first place. For East Asian parids, the extremely low levels of fixed SNPs between all HA and all LA birds as well as the commonality of low genomic diversity in both HA and non-HA western parids appear to reject this alternative hypothesis. Instead, landscape-wide differences in genomic diversity are more likely the result of different demographic histories between highland and lowland tits, driven probably by Pleistocene glacial patterns.

Shared Adaptive Pressures Result in Parallel Evolution at Functional, Not Genic, Levels.

For small-bodied and active endotherms, such as East Asian parids, the extreme elevation found across much of the QTP and the Himalayas poses considerable physiological challenges to survival. Across seven HA–LA comparisons of East Asian parids, we found that the most divergent genomic regions were significantly overrepresented by genes with biological processes and functions related to hypoxia and thermoregulation (Fig. 4D). Because the chromosomal locations of these outlier regions show little overlap among HA–LA comparisons (Fig. 3A) and given the high level of chromosomal synteny known generally in birds (55) and more specifically in tits (30, 33), our results appear to support the notion that little overlap exists among the particular genes involved in high-elevation adaptation. Indeed, we identified a total of 230 genes across all East Asian tits that may be involved in high-elevation adaptation, but on average, only 10.9 of these genes cooccurred in any two HA–LA comparisons, and just one gene was found in outlier windows shared among six of our seven comparisons. This result probably mirrors hemoglobin adaptation in HA birds: While occasional cases of parallel evolution at particular hemoglobin amino acid substitutions have been observed (56), more generally, derived hemoglobin proteins in HA birds have parallel functions that are the result of distinct molecular evolutionary solutions (7, 50, 57). Such genetic nonparallelism has been regarded to be probably caused by the low levels of heterozygosity of replicated species (1517). Our comparisons do not show increased genic parallelism across the three widespread parids with higher heterozygosity compared with the four HA endemics. We did, however, find some evidence that FST outliers may be more frequently associated with selection in HA endemics compared with widespread taxa, potentially due to more long-term adaptive divergence in endemics. Indeed, the pattern of limited overlap of HA-related genes among pairs of both HA endemics and widespread taxa suggests that other aspects of genomic variation or limited number of genes identified by small sampling may lead to little genic parallelism across parids. Likewise, the high similarity of enriched HA-related GO terms between HA endemics and widespread taxa, which is significantly higher than expected by chance, is likely consistent with the notion of genomic rather than genic parallelism (10, 12, 58).

Rapid Metabolic Adaptation to High Elevation.

In East Asian parids, circulatory and metabolic pathways were enriched in nearly all comparisons, whereas respiratory and tissue diffusion pathways were seldom enriched. Importantly, we found that both carbohydrate and lipid metabolic processes were enriched in all seven HA–LA comparisons; these processes consistently show the importance of avian adaptation to high elevations (10, 12, 30, 31, 58). Recently, it has been argued that the evolution of carbohydrate metabolism precedes adaptation for increased lipid metabolism among high-elevation vertebrates (10, 31, 59). In contrast, we find similar numbers of enriched carbohydrate-metabolism HA-related genes and lipid-metabolism HA-related genes in both HA endemics, which presumably have inhabited the QTP for millions of years, and HA populations of widespread parids, which separated from eastern populations only during the early to Mid-Pleistocene (Fig. 3B). This suggests either that adaptive evolution of lipid and carbohydrate metabolism is likely coincident rather than serial in East Asian parids and/or that the pace of high-elevation metabolic adaptation is relatively rapid, as has been shown recently for another East Asian songbird (58) and Andean waterfowl (10). The absence of a difference in the number of HA-related genes between HA endemics and widespread taxa is consistent with a hypothesis of rapid high-elevation adaptation.

In conclusion, East Asian parids show parallel diversity and demographic patterns in response to Pleistocene climate change and parallel genomic adaptation to extreme HA environments. However, the particular genes for high-elevation adaptation appear unique to each HA species. Additionally, we find no evidence for higher genic parallelism in widespread taxa compared with western endemics despite significantly lower heterozygosity in both coding and noncoding genetic regions in endemics, contrary to the prediction that parallelism increases with standing variation. However, the taxonomic and geographic scope of our study led to small sample sizes, which limited our power to detect selection, which may thus affect our interpretation of parallel adaptation, especially at the levels of gene or amino acid substitution. Future studies require population-level sampling for further investigation of parallel evolution in high-elevation birds and should focus on noncoding regions (32), recognizing that genetic adaptation to high elevation is likely polygenic, with many loci exerting relatively minor effects (58, 60). Researchers should also consider at what scale (SNP, gene, or genome) adaptation is likely to occur when addressing the question of parallelism or convergence in evolutionary processes.

Materials and Methods

Sampling and Sequencing.

We sampled 91 individuals of parids and penduline tits. Field sampling conformed to the National Wildlife Conservation Law with permission from the local forestry administration. We extracted total genomic DNA from muscle and blood samples using the Tissue/Cell Genomic DNA Extraction Kit (Aidlab Biotechnologies). We prepared genomic libraries with a mean insert size of 350 base pairs (bp) according to the manufacturer’s protocols, which were sequenced on an Illumina HiSeq X Ten machine to obtain 150-bp paired end reads at Berry Genomics. Prior to downstream analyses, we cleaned the reads by removing adapter sequences, as well as reads with ≥3% unidentified nucleotides and reads with over 50% of bases having Phred quality scores <3. Detailed information, including sample identification number, collection location, sequencing depth of coverage, National Center for Biotechnology Information (NCBI) Sequence Read Archive accession numbers, and more is available in Dataset S1.

Read Aligning and Variant Calling.

We aligned reads to the P. major reference genome (Assembly: GCA_001522545.2) (33) using BWA version 0.7.15 (61) with default parameters. Alignments were reordered, sorted, and marked for PCR duplicates using PICARD version 2.8.3 (broadinstitute.github.io/picard). We additionally performed local realignment for insertion-deletion (InDel) polymorphisms using the Genome Analysis Toolkit (GATK) version 3.7 (62). Variants were first called across all samples independently using both HaplotypeCaller in GATK and mpileup in Samtools version 1.3.1 (63). Sites identified as variants in either tool were subjected to a second round of variant calling in HaplotypeCaller to confirm genotypes.

To further improve variant quality, we conducted hard filtration for raw variant calls. We retained only variant sites for downstream analyses when they met the following criteria: quality > 1,000, 900 < overall depth (DP) < 2,700, quality by depth (QD) > 10, mapping quality (MQ) > 55, Fisher strand bias (FS) < 10, mapping quality rank sum (MQRankSum) > –5, read position rank sum (ReadPosRankSum) > –10, and symmetric odds ratio (SOR) < 2. Cutoffs for these parameters were chosen based on their density distributions. We also removed short InDels and retained only SNPs that had <10 individuals with missing data.

Phylogenetic Construction and Dating.

We reconstructed maximum-likelihood (ML) phylogenetic trees separately for autosomal and Z-chromosome SNPs (64), and we also generated one coalescent-based species tree based on 50 subsamples from the genome-wide SNP set (∼50,000 SNPs per replicate). For both the autosome and Z-chromosome datasets, we concatenated all SNPs and constructed ML phylogenies using RAxML version 8.2.10 (65) with 100 bootstrap replicates under the ASC_GTRGAMMA model. We also constructed a species tree using a coalescent-based method in ASTRAL-II (66) based on 50 random subsets of ∼50,000 SNPs. All tree-building approaches generated a single topology (Results), so we dated nodes on this tree as follows. Following Jarvis et al. (67) (see SI Appendix, Methods for details), we generated a dataset of clock-like evolving protein coding genes (Results) and enforced the above topology, which was loaded into the MCMCtree program in PAML version 4.9 (68) and generated both uncorrelated and autocorrelated relaxed clocks. We estimated posterior distributions of divergence time by Markov chain Monte Carlo sampling, with samples drawn every 2,000 steps over a total of 108 steps after a discarded burn-in of 2 × 107 steps.

Genetic Diversity Calculation.

We calculated nucleotide diversity (π) and heterozygosity (H) for each species or population in VCFtools version 0.1.13 (69). We used individual heterozygosity across all species and populations to simulate the spatial distribution of lineage-wide genetic diversity using kriging routines (70) in ArcGIS version 10.0 (ESRI) based on the heterozygosity at each locality. We then tested differences in the heterozygosity of genome-wide (both noncoding and coding) regions across western endemics, widespread species, and eastern endemics.

Demographic History Inference.

For each species and geographic population, we selected the individual with the highest sequencing depth to infer population size over time using PSMC (71). We generated consensus sequences that were filtered by removing sites with mapping quality below 25, low-inferred consensus quality (below 20), and low-read depth (less than 6, one-third of the mean coverage). According to suggestions by Nadachowska-Brzyska et al. (21), PSMC parameters (−p and −t options) were set as “−N25 −t5 −r5 −p4 + 30 × 2 + 4 + 6 + 10” with 100 bootstrap replicates. PSMC results were scaled using species-specific mutation rates and generation times (SI Appendix, Table S5). The synonymous substitution rate per synonymous site was calculated using the codeml program of PAML (72) based on the time tree (Fig. 1B) and was used to calculate the species-specific mutation rate (μ) for each species. We used the age of sexual maturity (73) multiplied by two as a proxy for the generation time (g) (74).

We scaled the relative population size to the effective size (Ne) using the approximated mutation rate (μ) based on the equation Nk = (θ0 × λk) / 4μs provided by H. Li and R. Durbin (github.com/lh3/psmc). We used generation time to convert coalescent time to calendar time using the following equation [Tk = (θ0 × tk × g) / 2μs]. θ0, tk, and λk are estimated by the PSMC model, and s is the bin size, which we set at the default value (s = 100). We used this equation to generate Ne estimates three times in the past 0.012 Mya (end of the LGP), 0.023 Mya (beginning of the LGM), and 0.11 Mya (beginning of the LGP). We used Wilcoxon rank-sum tests to test for significant differences between Ne estimates at the beginning of the LGP or LGM and at the end of the LGP (SI Appendix, Table S3).

Linkage Disequilibrium Analysis.

To select a suitable window size for genome-wide scanning, we estimated the levels of linkage disequilibrium (LD) within species for all species with a sample size of at least five. We phased haplotypes for each species using BEAGLE version 4.1 (75). Levels of LD (r2) were then calculated using VCFtools between pairs of loci in every 50-kb window of the genome. Haplotype r2 values were used to measure LD decay by averaging the r2 values in each bin of 500 bp over 20 kb. LD decreased rapidly within 5 kb and then slowly declined until 15 kb, suggesting that a window size of greater than 15 kb would avoid autocorrelation between adjacent windows (SI Appendix, Fig. S7).

Genome-Wide Scanning and Annotation of SNPs.

We scanned the whole genome to identify adaptive genetic variation at different levels using several FST analyses, including both SNP- and region-based FST calculations. We calculated SNP-based FST between all HA and all LA individuals, between HA endemics and their LA relatives, between HA and LA populations of widespread taxa, and between each HA species and its LA relative using VCFtools. We also calculated AF for all HA individuals and all LA individuals. Additionally, for each of the seven HA–LA pairs, we calculated FST in 50-kb sliding windows, a window size >3× greater than indicated by the maximum observed LD (SI Appendix, Fig. S7). We Z-transformed windows of autosomes and Z-chromosomes separately. We defined outlier windows as those with ZFST values exceeding the 99th percentile of the distribution. Therefore, each pair had 186 autosomal outlier windows and 16 Z-chromosome outlier windows. Meanwhile, we estimated the divergence patterns for these outlier windows using the dXY statistic calculated by an egglib_sliding_windows.py Python script (github.com/johnomics) and defined windows in the top 1% percentile of dXY values as dXY “islands” and the bottom 1% windows as dXY “valleys.” For comparison purposes, we generated 20 replicates of randomly selected sets of 186 autosomal and 16 Z-chromosome windows for each pair. For both the outlier and random replicate window sets, we annotated all SNPs using SnpEff version 4.3 (76) based on the location of SNPs in the P. major genome. SnpEff determines whether SNPs are found in coding or noncoding regions, and for coding-region SNPs, it determines whether the substitution is synonymous or nonsynonymous.

Functional Enrichment and Similarity Analyses.

Because they are metabolically active and homeothermic, most birds living at HA cannot compensate for hypoxia and low temperatures by behavioral change and must adapt physiologically to reduced oxygen availability and cold stress (7, 37). Physiological adaptations are most likely to occur at the four stages of the oxygen transport cascade (breathing/pulmonary oxygen diffusion, circulation oxygen delivery, tissue oxygen diffusion, and tissue oxygen utilization), in which oxygen is transferred via either diffusion or convection between atmospheric air and oxidative phosphorylation (37) or at ATP/carbohydrate/lipid metabolic processes (note: these three metabolic processes are also part of the oxygen transport cascade: tissue utilization) that are relevant to thermogenesis. Prior to all analyses, we identified GO molecular functions and biological processes associated with these four stages (48, 7779). This resulted in a list of 16 parental GO terms, to which we added an additional parental term directly associated with hypoxia (GO: 0036293; Fig. 3B). We used the QuickGO database (www.ebi.ac.uk/QuickGO/) (80) to identify direct descendants of these 17 parental GO terms (child terms) and filtered out only the GO terms found in the union of human, mouse, and chicken gene ontologies, which resulted in a final set of 973 GO terms (Dataset S2). We used topGO version 2.32.0 package (81) to annotate and enrich the SnpEff-identified variable coding genes by GO terms; for comparison, we also generated 100 random replicate gene sets (with the same count of genes as in the outlier set [e.g., n = 230]) from nonoutlier windows. From this, we separated all variable genes into two groups: those assignable to one of 973 potential high-elevation GO terms (HA-related) and those not assignable to the 973 GO terms, which we defined as “HA-unrelated.” The HA-related genes were used to construct a KEGG network in the CYTOSCAPE plugin ClueGO version 2.3.5 with kappa statistics (82). The network specificity was set to “medium,” and network connectivity was set to a kappa score of ≥0.4. Finally, we used topGO to calculate whether certain GO terms were “enriched” (i.e., were overrepresented among the list of outlier window genes more than could be expected by chance). We used a cutoff P value of 0.05 corrected for the false discovery rate using the Benjamini and Hochberg method (83).

For the set of seven HA–LA paired taxa, we compared 1) the number of coding genes; 2) the number of coding genes associated with the 973 HA-related GO terms; and 3) the number of enriched HA-related GO terms found in outlier windows with the average found in the 20 random replicates using two-sample t tests. To test whether the same HA-related enriched GO terms occurred in each of the seven outlier comparisons more often than expected by chance, we conducted pairwise hypergeometric tests using the dhyper function in R; these results were compared with the outcome of pairwise hypergeometric tests between the 20 random window replicates. Meanwhile, we measured the semantic similarity between HA-related enriched GO term sets using the GOSemSim version 3.11 package (84). The semantic comparison provides a quantitative way to compute similarities between GO terms, including functional similarity, interaction, and affiliation. Additionally, to compare these results with null expectations, we repeated the hypergeometric and semantic similarity analyses but with “HA-unrelated” enriched GO terms. Finally, we estimated differences in hypergeometric test values and semantic similarity between HA endemics and widespread taxa for both HA-related and HA-unrelated GO terms.

McDonald–Kreitman Test.

We constructed consensus sequences for each HA-related gene from HA–LA comparisons using vcf-consensus (69) based on our SNP data and aligned each HA-related gene using MUSCLE version 3.8.31 (85). Based on these alignments of HA-related genes, we conducted the McDonald–Kreitman test (47) in PopGenome version 2.7.2 (86) to detect signals of positive selection by comparing the ratio of fixed nonsynonymous sites (Dn) and fixed synonymous sites (Ds) with the ratio of polymorphic nonsynonymous sites (Pn) and polymorphic synonymous sites (Ps) using Fisher’s exact test. Under positive selection, the ratio of Dn to Ds is significantly larger than the ratio of Pn to Ps.

Ecological Analyses.

The atmospheric partial pressure of oxygen (PO2) was calculated by a simple formula based on altitude: PO2 = 21% × 101.325 × (1 – elevation / 44,300)^5.256, assuming a constant atmospheric pressure (101.325 kPa) and temperature (15 °C) at sea level. To estimate the range size of parid species, we calculated the geometric area of distribution maps downloaded from BirdLife International (datazone.birdlife.org) in ArcGIS. We estimated the annual temperature for each sampling locality based on 19 bioclimatic variables with a 30-s resolution from WorldClim (87) using the “Zonal” tool in ArcGIS.

Read more here: Source link