De Novo Assembly of a Draft Genome for the Parrotbill.

We first sequenced and assembled de novo a reference genome of a male Taiwanese vinous-throated parrotbill. The assembled genome, 1.06 Gb in length (genomic coverage = 110×, SI Appendix, Table S1), contains 6,508 scaffolds (2,034 scaffolds > 10 kb) with N50 size of 1.94 Mb. This draft genome includes a total of 72 Mb (6.8%) of repetitive sequences (SI Appendix, Table S2). We found 20,149 predicted genes, 9,581 of which could be functionally annotated. In total, we obtained 278 Mb of genes (26% of total genome) consisting of 253 Mb of introns and 25 Mb of coding regions.

Genetic Variation and Geographic Structure of Parrotbills in Taiwan.

We resequenced the genomes of 40 parrotbill individuals that were collected from four populations (10 individuals per population; Fig. 1A) with an average depth of 5.8× genome coverage (SI Appendix, Table S3). They formed two high- (>2,000 m) and low- (<100 m) altitudinal comparison pairs, one from the east side of the Central Mountain Range (CMR) in Taiwan and one from the west side. The altitudinal gradients between the two high- and low-altitude populations are steep (43 km and 55 km apart for the east and west population pairs, respectively). The two low-altitude localities were much warmer (over 10 °C higher annual mean temperature) and drier (50% lower annual precipitation) than the two high-altitude localities (SI Appendix, Table S3).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

(A) The vinous-throated parrotbill and four sites (red dots) on an east–west section of central Taiwan at which vinous-throated parrotbills were sampled and (B) the distribution of the F_ST and ΔF_ST in the east and west high-/low-altitude local population pairs of each 10-kb nonoverlapping genomic window that was aligned with the published genome of the zebra finch. Red dots on top of and within each panel represent candidate regions on the genome (n = 24); red horizontal lines indicate the top 1% of F_ST and ΔF_ST. EH, high-altitude population east of CMR; EL, low-altitude population east of CMR; WH, high-altitude population west of CMR; WL, low-altitude population west of CMR.

Comparing with the reference genome, 20,017,403 polymorphic sites (S) within the 40 resequenced genomes were found. The levels of intrapopulation genetic diversity are similar (π and mean individual heterozygosity, H¯, SI Appendix, Tables S3 and S5). However, the number of polymorphic sites (S), a demographically sensitive index, is about 12% greater in the low-altitude (mean S = 11,329,479) than in the high-altitude populations (mean S = 10,132,251). This is consistent with the conjecture that populations at or near glacial refugia retain more stable and larger effective population sizes (N_e) (hence more genetic variants) than populations at the expanding edges during postglacial expansion (20, 21).

The population structure was revealed in an unrooted neighbor-joining population tree (SI Appendix, Fig. S3) based on genome-wide F_ST (fixation index, a genetic measurement of population differential) values among all four populations (SI Appendix, Table S6), which indicates that pairs of populations on the same side of the CMR are clustered together. Thus, the two population pairs can serve as independent samples for illustrating parallel adaptation along altitudinal gradients (22).

Demographic History of the Vinous-Throated Parrotbill in Taiwan.

Bottleneck and other events could leave genetic signatures resembling those of selection (23). We used the diploid sequence of the reference genome and a high-coverage genome sequence (24.5×) of an additional male Taiwanese parrotbill to test for recent severe bottleneck events. Results of pairwise sequentially Markovian coalescent (PSMC) analysis (24) indicate that their N_e varied through time, although without evidence of a severe population bottleneck (SI Appendix, Fig. S4). N_e was ∼30,000–100,000 from the middle of the last glacial period (∼50–70 thousand years ago) to the present.

Genomic Regions Associated with Divergent Selection Between Different Altitudes.

For the east population pair, we found 2,812 outlier regions (0.715 ≥ F_ST ≥ 0.271; the top 1% of F_ST values; mean genome-wide F_ST = 0.108, 95% CI = 0.108–0.109); for the west population pair, we found 3,051 outlier pairs (0.536 ≥ F_ST ≥ 0.140; mean genome-wide F_ST = 0.052, 95% CI = 0.051–0.052) (Fig. 1B). Of these regions, 152 are shared between the east and west populations. It is important to consider that both ecologically divergent selection and purifying selection would reduce genetic diversity on functionally conserved regions and lead to high F_ST in the linked genomic region. Whereas divergent selection should only affect populations in different environments, by contrast, purifying selection would affect all populations (25). For the Taiwanese parrotbills, F_ST of 10-kb windows was highly correlated between all population pairs (SI Appendix, Fig. S5), implying that the high-F_ST regions across genome is largely the result of purifying selection. Purifying selection on linked genomic regions is further supported by the negative correlation between F_ST and π values (SI Appendix, Fig. S6), positive correlation between F_ST values and the level of linkage disequilibrium, r² (SI Appendix, Fig. S7), and negative correlation between π and r² values (SI Appendix, Fig. S8). Therefore, we employed ΔF_ST analysis to reduce the effect of purifying selection on genetic differentiation of genomic regions.

Only 24 10-kb nonoverlapping autosomal regions in the top 1% for both F_ST and ΔF_ST values (F_ST ≥ 0.464 and 0.344 in the east and west population pairs, respectively; and ΔF_ST ≥ 0.164 and 0.101 in the east and west population pairs, respectively) were considered candidate regions for adaptation to altitudes along both sides of the CMR (Fig. 1B). These regions have high levels of genetic differentiation between different altitudes and are less likely to be affected by purifying selection. Only six of the candidate regions are concatenated into two longer divergent genomic islands that were 20 and 40 kb in length.

Seventeen candidate genes were identified (SI Appendix, Table S8). Among them, four are involved in oxygen cascade: VAV3 and COL15A1 are related to angiogenesis; IGF2 to respiratory system phenotype; and TPPP to hemoglobin content. Three candidate genes are related to aspects of metabolism that could be involved in thermoregulation: SUPT7L (26) and HBP1 are related to lipid metabolism; and OLA1 to ATP hydrolysis. High‐altitude mammals have physiological adaptations in their circulatory, hematologic, respiratory, and thermoregulatory systems as well as in their metabolism in comparison with their sea-level counterparts (27, 28). Therefore, the genes we identified may be associated with adaptation between high and low altitudes. Among them, VAV3 (29), TPPP (30), IGF2 (31), and COL1A1 (32) have been reported to be associated with altitudinal adaption or hypoxia syndrome in other organisms. We also found that HBP1 (33) and OLA1 (34) are associated with temperature acclimation and heat shock phenotypes, respectively, in other vertebrates. The phenogenotypic associations of these genes in other systems not only validate our inferences, but also support the view that multiple quantitative traits, including oxygen cascade and thermal regulation, were involved in the adaptation to different altitudes.

Although adaptation could have been achieved through regulatory evolution of individual genes (e.g., those mentioned above), functionally related genes could also have changed expression simultaneously in the adaptive process. To examine this possibility, we performed an enrichment analysis on gene ontology (GO) terms. Among the 23 GO terms in which the 17 candidate genes were enriched [false discovering rate (FDR) < 0.05; P < 0.0019], many were found to be associated with glycogen metabolism (7 GO terms) and epigenetic control of gene expression (11 GO terms), especially histone acetylation (8 GO terms) for the latter group (SI Appendix, Table S9). This observation is consistent with previous studies showing reduced activity of glycogen metabolism accompanied by increased metabolism of fatty acids in organisms moving to a higher altitude (35, 36). This result also implies that epigenetic control, especially histone modification involving acetylation, played a critical role in the adaptation of vinous-throated parrotbills to different altitudes and their ability to cope with rapid environmental changes (37). Consistently with that possibility, the involvement of chromatin modifications in high-altitude adaptation has been observed in several mammals (38, 39).

Predominant Role of Noncoding Genomic Regions in Adaptation to Different Altitudes.

Within the 24 candidate regions, 35 (average F_ST = 0.767; range: 0.631–1.000) and 32 (average F_ST = 0.628; range: 0.461–0.832) SNPs with the highest F_ST values in the east and west altitudinal pairs, respectively, were considered as candidate SNPs (SI Appendix, Tables S10–S13). The F_ST values of these candidate SNPs are significantly higher than those for all SNPs in both altitudinal pairs (Welch two-sample t test: t = 19.795 in the east pair and 21.592 in the west pair, P < 1 × 10⁻¹⁵). Furthermore, we found that most of the candidate SNPs are fixed or nearly fixed (frequency of major allele ≥ 0.9) in one of the four populations (34 and 27 of the candidate SNPs for the eastern and western populations, respectively; SI Appendix, Tables S10–S13), suggesting that they may have had major effects on the associated phenotypes. The higher polymorphism levels of the other candidate SNPs suggest that they are probably under soft sweep or polygenic selection and thus unable to drive these advantageous genetic variants to fixation (40). Noticeably, among these candidate SNPs, only 10 (29.4%) and 4 (14.8%) of them were fixed or nearly fixed in both high and low altitudes in the east and west populations, respectively. This implies that the intensity of divergent selection could be asymmetric at different altitudes for these candidate SNPs.

Intriguingly, none of these candidate SNPs is located within coding regions: eight are located in intronic regions of the HBP1 and AAK1 genes, and the other 59 are located in intergenic regions. The disproportionate number of candidate SNPs located in intergenic regions (Fisher’s exact test: P = 0.008) supports the predominant role of regulatory regions in adaptation (41, 42). Alleles in regulatory regions are often codominant in their effect, have mild pleiotropic consequences, and modify only the expression of genes rather than their functions (43).

Standing Genetic Variation as the Predominant Source of Adaptation.

The prerequisite for inferring standing variation using the comparative approach is recent coancestry between sampled taxa. We found that the mainland and Taiwan parrotbill populations diverged early in the last glacial period with estimated N_es of 482,065 for Taiwan and 1,527,174 for the mainland (SI Appendix, Table S14). Their recent divergence is also supported by the low genome-wide net genetic distance (D_a) between them, 0.0007 (D_xy = 0.0051 between Taiwan and mainland populations; πs of Taiwan and mainland populations are 0.0040 and 0.0048, respectively). This implies a divergence time (t = D_a/2μt; t, divergence time; μ, mutations per site per year) of about 152,000 y, meaning that genetic variants shared between the extant Taiwan and mainland parrotbill populations are unlikely to have arisen independently. However, the F_ST values between the mainland population and the low-altitude populations in Taiwan are lower (mean = 0.122) than those between the mainland and the high-altitude ones (mean = 0.133). This supports the view that the high-altitude populations are the leading edge of a postglacial population that possesses lower genetic diversity and a smaller N_e than the low-altitude refuge population.

Assuming genetic panmixia across the vinous-throated parrotbill’s mainland range (44), we further resequenced the genomes of 40 parrotbill individuals (SI Appendix, Table S4) to represent the species’ genetic diversity in the mainland. We found ∼1.3 times as many polymorphic sites in the mainland (S = 27,879,100) as in the Taiwan population (S = 20,017,403). About 66.2% of polymorphic sites in Taiwan (S = 13,259,879 sites) were also polymorphic in the mainland. These shared polymorphic sites likely represent common standing SNPs (minimum allele frequency, MAF, of SNPs ≥ 1.25%) in both the mainland and Taiwan populations, while the private genetic variants of each population are more likely to be rare standing SNPs (MAF < 1.25%) that our sampling scheme failed to detect in the other population, together with some new mutations that occurred since the populations split. About 94% of the candidate SNPs in the Taiwan populations are also polymorphic (or shared) in the mainland populations (Fig. 2A), significantly more than the 66.2% of all Taiwan SNPs (Fisher’s exact test: P = 1.6 × 10⁻⁵ and P = 0.002 for the east and west altitudinal pairs, respectively; Fig. 2A). The genomic regions around functional elements including genes should experience different evolutionary trajectories from those around nonfunctional regions because of the distinct densities of selection targets within them (45). Therefore, we compared the proportion of shared SNPs in a range 1 Mb upstream and downstream of all genes in the genome with the proportion of the candidate SNPs themselves that were shared. The proportion of candidate SNPs that was shared was still significantly higher (Fisher’s exact test: P < 0.001 and P = 0.004; Fig. 2A), confirming that shared standing genetic variation was the predominant genetic source of adaptation.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

(A) The proportion of polymorphic SNPs shared with the mainland population is significantly higher for candidate SNPs inferred in both high-and low-altitude population pairs than that of the entire genome and SNPs within 1 Mb downstream and upstream regions of all genes in the genome (within 1 mb) of the Taiwan population (East, east population pair; West, west population pair; Fisher’s exact test, candidate SNPs: P = 1.6 × 10⁻⁵ and P = 0.002 for the east and west high-/low-altitude population pairs, respectively; all of the SNPs within 1 Mb regions: P < 0.00001 for both east and west high-/low-altitude population pairs). (B) The minimum allele frequency (MAF) of candidate SNPs in the mainland population is significantly higher than that of noncandidate SNPs (Welch two-sample t test: t = −10.524, df = 63, P = 1.644 × 10⁻¹⁵). (C) The MAF of shared variants is significantly higher than private variants (Welch two-sample t test, t = −2,909.8, df = 19,750,000, P < 2.2e-16). PV, private variants; SV, shared variants.

In the mainland parrotbill population, the MAF of candidate SNPs (average MAF = 0.300, σ² = 0.021) was significantly higher than that of other SNPs (average MAF = 0.107, σ² = 0.016; Welch two-sample t test: t = −10.524, df = 63, P = 1.644 × 10⁻¹⁵; Fig. 2B). Assuming that the mainland population with its large long-term N_e (SI Appendix, Table S14) is more likely to retain the ancestral allelic frequencies, our results support the view that standing variants with high initial frequencies facilitate a swift response to selection posed by changing environments (46, 47). Because the N_e of parrotbills is relatively large in Taiwan, genetic drift should only have a limited effect on the dynamics of allelic frequency in the population. Because the average MAF of private SNPs in Taiwan (average MAF = 0.038, σ² = 0.002) is significantly lower than that of shared ones (average MAF = 0.158, σ² = 0.019; Welch two-sample t test, t = −2909.8, df = 19,750,000, P value < 2.2e-16; Fig. 2C), these private SNPs are likely to have preexisted in the ancestral population in low frequency. The low proportion of private SNPs in candidate SNPs implies that rare genetic variants in the ancestral population are less likely to enable response to rapid environmental changes as predicted by the theory (4).

The significant role of shared standing genetic variation in adaptation revealed in this study (Fig. 2) has important implications for our understanding of adaptation and biological conservation. Our results suggest that to cope with environmental changes quickly, a population may benefit significantly from the short waiting time needed to exploit standing genetic variation, especially for variants with high initial frequencies. Our results are consistent with the accumulating evidence that species with high levels of genetic diversity can persist in changing environments (48) and that they are more likely to colonize novel environments (49). Therefore, a species’ evolutionary potential can be constrained by a lack of genetic variation (50). This underscores the need to preserve the genetic diversity of species in the face of accelerating environmental changes. Conversely, species with low genetic variation could fail to cope with accelerating anthropogenic environmental changes, which warrants more attention to their conservation status.

Both High and Low Altitudes Could Have Been Harsh for the Parrotbill Since the End of the Last Glacial Period.

During the course of postglacial range expansion, parrotbills at the leading expansion edge are thought to have encountered different environmental challenges from those at the trailing edge, and high altitudes are conventionally considered to have harsher environments than low altitudes (11). We would thus expect high-altitude parrotbills to be under more stringent recent selection. To test the “harsh high altitudes” hypothesis, we examined the polymorphic patterns of the candidate SNPs at different altitudes and found similar numbers of fixed or nearly fixed candidate SNPs on both sides of the island (Fisher’s exact test P = 0.809 and 0.803 for the east and west pairs, Fig. 3). Assuming that fixation probability depends on the strength of directional selection in a given environment, this suggests that the strength of selection pressures at different altitudes was comparable. Evidence of strong selection pressure in the low-altitude populations also arose from the results of the Tajima’s D test that we used to detect recent selective sweep [<0.1 N_e generations (51, 52)], in 10-kb sliding regions across the genome. Recent selective sweeps that occurred since the last glacial period as indicated by the lowest 1% of Tajima’s D value in each local population were detected for 6 (Tajima’s D ranges from: −1.3569 to −1.6825) out of 24 candidate regions, with similar frequencies in population pairs on both sides of the CMR (none and four for the east high- and low-altitude populations, respectively; one and two for the west high- and low-altitude populations, respectively; Fisher’s exact test: P = 0.112 and P = 1.000 for the east and west population pairs, respectively).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

The proportion of candidate SNPs that are fixed or nearly fixed (90% ≤ the major allele frequency < 100%) in the low- and high-altitude populations.

Our results suggest that during postglacial expansion, the population remaining at low altitudes (the trailing edge of the expansion) could also have experienced strong selection pressure. These low-altitude selection pressures may have been caused by increased temperatures, whereas those at the high-altitude leading edge of the expansion arose from the lower oxygen partial pressure. This finding calls for more research on adaptation at low altitudes or latitudes where populations at the trailing edge of postglacial expansion are located and are considered as long-term reservoirs of species’ genetic diversity and cradles of speciation (53).