Genomic architecture underlying morphological and physiological adaptation to high elevation in a songbird

Organisms often acquire physiological and morphological modifications to conquer ecological challenges when colonizing new environments which lead to their adaptive evolution. However, deciphering the genomic mechanism of ecological adaptation is difficult because ecological environments are often too complex for straightforward interpretation. Thus, we examined the adaptation of a widespread songbird—the rufous‐capped babbler (Cyanoderma ruficeps)—to a relatively simple system: distinct environments across elevational gradients on the mountainous island of Taiwan. We focused on the genomic sequences of 43 birds from five populations to show that the Taiwan group split from its sister group in mainland China around 1–2 million years ago (Ma) and colonized the montane habitats of Taiwan at least twice around 0.03–0.22 Ma. The montane and lowland Taiwan populations diverged with gene flow between them, suggesting strong selection associated with different elevations. We found that the montane babblers had smaller beaks than the lowland ones, consistent with Allen's rule, and identified candidate genes—COL9A1 and SOX11—underlying the beak size changes. We also found that altitudinally divergent mutations were mostly located in noncoding regions and tended to accumulate in chromosomal inversions and autosomes. The altitudinally divergent mutations might regulate genes related to haematopoietic, metabolic, immune, auditory and vision functions, as well as cerebrum morphology and plumage development. The results reveal the genomic bases of morphological and physiological adaptation in this species to the low temperature, hypoxia and high UV light environment at high elevation. These findings improve our understanding of how ecological adaptation drives population divergence from the perspective of genomic architecture.


| INTRODUC TI ON
When species adapt to new habitats and accumulate genetic changes that contribute to their local adaptation, they may show genomic incompatibility between diverging populations that may then evolve into new species (Nosil & Feder, 2012;Schluter, 2001). However, adaptation-driven divergence (Nosil, 2012;Schluter, 2009), which seems conceptually straightforward, can be difficult to perceive and interpret in nature because ecological contexts are often intricate. For species that are distributed along an elevational gradient, populations have to physiologically and/or morphologically adapt to distinct habitats across a continuum of environmental change.
Analysing species along this kind of gradient thus provides a simple system, ideal for examining the mechanism of ecological adaptation.
Animals at higher elevations may have smaller body extremities, such as bird beaks, for reducing heat loss (Symonds & Tattersall, 2010).
Temperature, together with other environmental factors such as oxygen level and high ultraviolet (UV) light that change with elevation, could have also caused physiological adaptations in respiratory, cardiovascular, immune and metabolic systems (Mishra & Ganju, 2010;Scott et al., 2008;Siebenmann et al., 2017). Although these phenotypic changes imply natural selection for beneficial alleles in the corresponding genes, it is difficult to tell how the evolutionary processes actually unfolded in wild populations.
Genomic sequencing provides a powerful tool not only for reconstructing divergence history, but also for testing hypotheses regarding the origin and association of mutations that contribute to physiological and morphological adaptation to new environments (Seehausen et al., 2014), including high elevations (Hu et al., 2022;Lai et al., 2019;Martini et al., 2021;Qu et al., 2013;Qu et al., 2021;Storz, 2021). Studies have shown that the arrangement of mutations in the genome is important to ecological adaptation because structural genomic variations often affect multiple genes together (Malinsky et al., 2015;Zong et al., 2021). Outstanding questions regarding the role of genome architecture in ecological adaptation include whether mutations in coding or noncoding regions tend to be subject to divergent selection for adaptation (Hu et al., 2022;Lai et al., 2019;Malinsky et al., 2015;Nosil & Feder, 2012;Wolf & Ellegren, 2017), whether chromosomal inversions facilitate the accumulation of adaptive genetic differences between diverging lineages (Kirkpatrick & Barton, 2006;Poelstra et al., 2014), and whether genes related to ecological difference tend to be located in sex chromosomes or associated with those under sexual selection that hence facilitate the evolution of reproductive isolation (Bourgeois et al., 2020;Edwards et al., 2015;Malinsky et al., 2015;Qvarnström & Bailey, 2009). These questions remain poorly understood, especially in nonmodel species, and we aim to examine them by studying altitudinal adaptation in species with wide elevational ranges.
Early studies on the genomics of ecological divergence were limited to a few well-studied species, such as threespine sticklebacks (Gasterosteus aculeatus, Jones et al., 2012), lyrate rockcress (Arabidopsis lyrata, Turner et al., 2010), stick insects (Timema cristinae, Gompert et al., 2014) and cichlid fish (family Cichlidae, Malinsky et al., 2015;McGee et al., 2016). Although there have recently been more studies exploring genomic differences between sister or closely related species of wild birds (Ellegren et al., 2012;Martini et al., 2021;Poelstra et al., 2014;Qu et al., 2021;Vijay et al., 2016), these well-differentiated species are at the extreme end of the speciation continuum, and thus provide little information regarding the process of population divergence towards speciation (Marques et al., 2016). Therefore, there is currently a demand for research on how genomic differentiation contributes to ecological niche shifts among intraspecific populations for more nonmodel species (e.g., Lai et al., 2019;Qu et al., 2020). Such research will provide profound insights into how environmental heterogeneity enriches species diversity in nature.
Taiwan is a subtropical island characterized by numerous high mountains, and most birds in Taiwan are either confined to one or a few bioclimatic zones or perform seasonal elevational migration. In contrast, the rufous-capped babbler (Cyanoderma ruficeps) has a wide elevational range, spanning several bioclimatic zones (from sea level to nearly 3,000 m above sea level [a.s.l.]) in Taiwan. Interestingly, the plumage of rufous-capped babblers shows different levels of UVreflectance and brightness, probably due to different light environments across elevations, suggesting altitudinal adaptation in this bird (Fang et al., 2022). The plumage divergence could also affect the visual system, since plumage traits are critical to intraspecific communication such as mate choice (Maia et al., 2016;Poelstra et al., 2014). If and how the phenotypic changes are reflected in genomic differentiation is unclear.
Studies have proposed several phenotypic changes that may facilitate animal adaptation to shifting environmental regimes along the elevational gradient (Mishra & Ganju, 2010;Scott et al., 2008;Siebenmann et al., 2017;Symonds & Tattersall, 2010). For instance, Allen's rule states that species that live in cold environments tend to have smaller body extremities relative to their overall body sizes than those living in warm environments (Allen, 1877). Because beaks could serve as thermoregulators in birds, positive relationships between beak size and environment temperature in several avian species have been observed (Fan et al., 2019;Laiolo & Rolando, 2001;Symonds & Tattersall, 2010;Tattersall et al., 2009). However, the molecular mechanism of beak size divergence in response to temperature change is still poorly understood and warrants further research. Examining the change in beak size along an elevation gradient in the rufous-capped babbler with an available reference genome  is a promising attempt to further understand this molecular mechanism. Altitudinal adaptation may also involve physiological changes such as those in respiratory and cardiovascular systems to survive in the hypoxic conditions at high elevation. In addition, UV light intensity and low temperature in montane areas could induce changes in immune and metabolic systems (Mishra & Ganju, 2010;Scott et al., 2008;Siebenmann et al., 2017). Because of these features, the rufous-capped babbler is thus a valuable system for examining the molecular mechanism and structural genomic variations underlying the altitudinal adaptation of animals during their initial stages of ecological divergence.
In this study, we resequenced genomes of rufous-capped babblers from Taiwan and mainland China to investigate the divergence history and the genomic signature of altitudinal adaptation in the Taiwanese populations. In particular, (i) we analysed the beak morphology and genomic data of rufous-capped babblers collected across elevations to test Allen's rule and explore the underlying molecular mechanism. We also (ii) investigated the genomic bases of other phenotypic changes associated with altitudinal adaptation, and (iii) examined whether genetic variation in coding or noncoding regions dominated altitudinal adaptation, (iv) whether chromosomal inversions were associated with elevational divergence and (v) whether altitudinally divergent genes accumulated in sex chromosomes and/or were associated with traits related to mate choice.
We aimed to decipher adaptive changes in response to environmental challenges from a perspective of genomic architecture. Studying how altitudinal adaptation affects population divergence is critical for understanding the evolution of species in montane areas and predicting their fates in light of ongoing global warming.

| Elevational differences in beak morphology
We examined whether rufous-capped babblers exhibit elevational variation in their beak morphology consistent with Allen's rule.
According to the rule, we expected thinner/shorter beaks and smaller beak surface areas relative to body sizes in montane birds than in their lowland counterparts. To test these predictions, in 2013-2019 we captured 116 rufous-capped babblers from a lowland (L1 in Figure 1a: 121.65°E, 24.15°N, 28 m a.s.l.; N = 32), midelevation (121.45°E, 24.21°N, 1,223 m a.s.l.; N = 61) and montane locality (H1 in Figure 1a: 121.30°E, 24.17°N, 2668 m a.s.l.; N = 23) along an elevation gradient. We used a digital caliper to take wing length, beak length, beak width and beak height to the nearest 0.01 mm from live birds, with which we then estimated beak surface area using the formula, ((width + height)/4) × length × 3.14159. To identify the gender of rufous-capped babblers-a sexually monomorphic species-we collected blood samples during bird banding to conduct molecular sexual typing in the laboratory (Method S1). Individual birds were released where they were caught after the above treatments.
We scaled aforementioned beak variables by an individual bird's wing length, which is the recommended measure for quantifying intraspecific body size variation of passerine birds (Gosler et al., 1998).
We then used two-way ANOVA to evaluate whether birds of different sexes showed different elevational trends in their beak morphologies. After normality evaluations with both quantile-quantile plots and Shapiro-Wilk tests ( Figure S1), we applied a rank transformation to each wing-scaled beak variable to fulfil the normality assumption of ANOVA as much as possible. We also conducted Levene's tests by which we confirmed the homogeneity of variance between groups (p > .05). We found significant intersexual trend differences only for beak height (p < .05). For this variable, we further used one-way ANOVA to evaluate elevational variations in males and females separately. Tukey post hoc tests were then performed once significant results from ANOVA were obtained. We conducted Tukey post hoc tests with the two sexes pooled for beak length, beak width and beak surface area, and separately for each sex for beak height.

| Whole genome sequencing
In addition to L1 and H1 described above, we sampled rufouscapped babblers from an additional lowland locality (L2 in Figure 1a

| Sample mapping and base quality score recalibration
The trimmed reads were mapped against a rufous-capped babbler reference genome  using the MEM model of bwa version 0.7.17 (Li & Durbin, 2009). To improve alignments around indels, we performed local realignments using gatk 3.6.0 (McKenna et al., 2010) with default settings. We then marked duplicated reads using the MarkDuplicates function of picard version 2.18.14 (http://broad insti tute.github.io/picar d/).

IndelRealigner
We performed multistep base quality score recalibration (BQSR) before final single-nucleotide polymorphism (SNP) calling. Due to the lack of reference SNP data sets for the rufouscapped babbler, we performed BQSR twice as follows to improve quality score recalibration: (i) variants were first identified using gatk 4.0.6 HaplotypeCaller (Poplin et al., 2018) and bcftools version 1.9 mpileup (Danecek & McCarthy, 2017), respectively; (iii) the concordantly identified variants were then filtered by gatk 4.0.6 VariantFiltration with the following criteria: we removed sites that had mapping quality scores of <40, variants calling confidence of <2 or signatures of sequencing strand preference (FisherStrand >60 and StrandOddsRatio >3.0); and (iii) based on the remaining variants' nucleotide base quality scores, we adjusted whole-alignment nucleotide base quality scores using gatk 4.0.6 BaseRecalibrator.

| SNP calling
We applied bcftools mpileup with default settings to the recalibrated alignment for final SNP calling. The obtained SNPs were further filtered with the following criteria: (i) multi-allelic variants were excluded using the SelectVariants function of gatk; (ii) SNPs were treated as missing data if the sequencing coverage was less than one-third or greater than two times the mean coverage of corresponding individuals; and (iii) SNPs containing >30% missing data over all samples were excluded from further analyses.

| Historical population demography and gene flow
To distinguish autosomal markers from Z-linked ones, we mapped babbler scaffolds  against the zebra finch (Taeniopygia guttata) assembly (bTaeGut1_v1.p from NCBI) using minimap2 (Li, 2018). We conducted pairwise sequentially Markovian coalescent (psmc, Li & Durbin, 2011) analyses on autosomal markers to infer temporal demographic changes of each of the Taiwanese and mainland Chinese populations. From each population, three individuals with the highest sequencing coverage were included in the demographic reconstructions. After exploratory runs under different parameter settings, we implemented formal runs with commands: -N30 -t5 -r5 -p "4+30*2+4+6+10" which steadily inferred Schiffels & Wang, 2020) analyses. Applied to the same individual birds of each population as in psmc, we ran msmc2 with the same patterning parameter -p "4+30*2+4+6+10" alongside the same nucleotide mutation rate and generation time settings. We used the flag -I to account for the unphased nature of the babbler genomes (Schiffels & Wang, 2020).
The above demographic analyses revealed a shared demographic pattern within Taiwan but not between Taiwan and mainland China (see Results), suggesting a common ancestry of the four Taiwanese populations. We further adopted a coalescent-based approach implemented in fastsimcoal version 2.6.0.3 (Excoffier et al., 2013) to estimate population divergence times and gene flow. We randomly selected 100,000 autosomal SNPs that segregated among the 43 individuals from the two geographical regions and prepared the fastsimcoal entry site frequency spectrum (SFS) using easysfs (https://github. com/isaac overc ast/easySFS). We assumed an isolation with migration (IM) divergence process (Hey, 2010), with gene flow between Taiwan and mainland China restricted to the period spanning from the Taiwan-China split to the most ancient population split within Taiwan. We specified the population tree topology as (China, ((L1, H1), (L2, H2))) based on our recent finding (Kuo et al., n.d.). To determine which of the two IM models (either L1-H1 or L2-H2 presents a more ancient split) was more compatible with the specified population tree, we performed fastsimcoal analyses on these two population pairs separately. According to the resultant divergence time estimates, we specified the full five-population IM model with a more ancient L2-H2 divergence than L1-H1. Based on the full model, we performed 100,000 coalescent simulations, followed by 40 cycles of conditional expectation-maximization computation (ECM cycles). After empirical tests, a -C 3 command was used to collapse low-count entries (≤3) of the SFS into a single category to improve composite likelihood estimation. The above optimization procedure was repeated 50 times, each starting with random initial parameter values. We took 10 runs that showed top composite likelihood values, with which we acquired the mean estimates of the focal parameters.
We then focused on the population structure within Taiwan.
This was assessed with a principal component analysis (PCA) based on 100,000 randomly selected autosomal SNPs from the four Taiwanese populations, carried out with the dudi.pca function in the adegenet package (Jombart, 2008).

| Genomic regions under altitudinally divergent selection
We calculated Nei and Li (1979) nucleotide diversity (π) and Hudson et al. (1992) F ST as measures of within-and between-population genetic variability, respectively. To evaluate variation across rufouscapped babbler genomes in terms of the two diversity measures, we used the R package popgenome version 2.7.1 (Pfeifer et al., 2014) to obtain measures along 10-kb nonoverlapping genomic windows.
Scaffolds of <10 kb (accounting for 1.26% of the genome) were excluded from diversity calculations.
Genomic regions that were subject to altitudinally divergent selection were expected to show large F ST values between populations from different elevations. However, large F ST values may not always reflect divergent selection. Alternatively, purifying selection, which removes deleterious mutations from populations, can also reduce genetic diversity and thus inflate F ST (Cruickshank & Hahn, 2014). Indeed, we found positive correlations between F ST values from cross-elevation population pairs and from same-elevation population pairs ( Figure  We identified windows that showed the top 5% ∆zF ST in each of the four montane-lowland population pairs (i.e., L1-H1, L1-H2, L2-H1 and L2-H2) and regarded the intersection windows as candidate genomic regions under altitudinally divergent selection (hereafter referred to as the "candidate regions" for simplicity). Note that markers obeying different inheritance modes can show varied F ST values simply due to their different effective population sizes rather than reflecting any selection signatures (Wilson Sayres, 2018). As such, we identified candidate regions for autosomes and the Z chromosome separately.
From the candidate regions, we searched for fixed and nearly fixed SNPs (hereafter referred to as "fixed SNPs" for simplicity), which reflected strong selection signals. Such fixed SNPs were each identified based on the following criteria: (i) major allele frequency ≥0.9 and (ii) allele frequency difference between lowland and montane populations ≥0.4 (top 2.5% of values in the candidate regions).
We pooled populations based on whether they were montane or lowland to calculate the allele frequencies of SNPs.

| Regulatory targets of fixed SNPs and their functions
We found SNPs either fixed in montane or lowland populations to disproportionally occur in noncoding genomic regions (see Results).
This implies that genomic regions under strong divergent selection probably most frequently served as regulators of gene expression such as promoters or enhancers. While a recent estimation revealed that, in most cases, the enhancers are <100 kb away from their target genes (Fulco et al., 2019), some enhancer-target interactions occur over longer distances (Vierstra et al., 2014). We identified genes whose start codons were within either 250 kb upstream or downstream of the fixed SNPs, and regarded these genes as putative regulatory targets of the selection-acting mutations. It is worth noting that given the lack of information regarding the sequences or locations of promoters and enhancers in this bird, we could not ascertain how the fixed SNPs influence the putative regulatory targets.
We examined whether the identified genes were enriched in biological function in terms of gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG). This examination was carried out based on orthologous relationships between babbler genes and genes from another five avian species obtained from Ensembl version 95: chicken (GRCg6a), duck (BGI_duck_1.0), flycatcher (FicAlb_1.4), turkey (Turkey_2.01) and zebra finch (tae-Gut3.2.4). Multispecies orthologous relationships among genes were inferred using sonicparanoid version 1.3.0 (Cosentino & Iwasaki, 2018).
Customized Perl scripts (see Data Availability and Benefit Sharing) were used for GO and KEGG enrichment analyses (one-tailed Fisher's exact tests) with false discovery rates (FDRs) controlled at <0.05.
We further conducted phenotype enrichment analysis implemented in modphea (Weng & Liao, 2017) to examine overrepresentations of the identified genes for specific phenotypic traits. Given that modphea only compiled phenotype annotations for six model species-human, mouse, zebrafish, fruit fly, roundworm and budding yeast-we chose mouse for our analysis due to its phylogenetically closer relationship to the babbler and more comprehensive annotations relative to human ones. Phenotype enrichment analyses were carried out based on gene orthology inferred using sonicparanoid that associated babbler genes to Ensembl version 95 mouse (GRCm38.p6) genes.

| Chromosomal inversions in babbler genomes
Inversions have been found to facilitate adaptive evolution and population differentiation of organisms (Kirkpatrick & Barton, 2006). To assess the effect of inversions on divergence between montane and lowland rufous-capped babblers, we first identified chromosomal inversions among the 40 Taiwanese samples using delly version 0.8.1 (Rausch et al., 2012). Among delly-inferred inversions, we focused on those with lengths of >10 kb. In addition, we applied the following filtering criteria to obtain highly confidential inversions: (i) supported with three or more read pairs that each had a mapping quality (MQ) of >20 (i.e., delly's "PASS" criterion); (ii) contained at least one split-read (delly's "PRECISE" criterion) such that the chromosomal breakpoint(s) can be precisely identified; and (iii) had a minor allele frequency of inferred inversion variants that was >0.05.
We examined whether the inversions tended to occur in altitudinally divergent genomic regions (i.e., the 10-kb windows showing high ∆zF ST ) to understand their potential role in rufous-capped babblers' elevational divergence. For this, we first classified windows as inversion-present and inversion-absent ones depending on whether the windows overlapped inversions. We then performed a onetailed Fisher's exact test to evaluate whether the inversion-present windows were more frequently altitudinally divergent ones. In addition, we used one-tailed Wilcoxon rank sum tests to assess differences between inversion-present and inversion-absent divergent windows in levels of genetic differentiation (F ST ) and net sequence divergence between populations (D a , Nei, 1987)-the latter was proposed for detecting selection signatures between populations with low levels of divergence (Noor & Bennett, 2009). Wilcoxon rank sum tests were applied to compare populations from different elevations as well as from the same elevations. To obtain window-wise D a values, we first used popgenome to calculate the total sequence divergence between populations (D xy , Nei, 1987). We then derived D a based on Hudson et al. (1992) F ST , with F ST = D a /D xy. Finally, we performed GO, KEGG and phenotype enrichment analyses as previously described but this time focusing on genes in the inversion-present divergent windows to investigate these genes' biological functions.

| Elevational differences in beak morphology
Two-way ANOVA indicated significant elevational differences in beak surface area, beak length and beak width (p < .05). Tukey post hoc tests revealed significantly longer and wider beaks and larger beak surface areas in lowland birds than in midelevation and montane birds (p < .05; Figure 2). One-way ANOVA indicated significant elevational variation in beak height in female birds but not in males (p < .001 and >.05, respectively). Tukey tests supported higher beaks in lowland females compared to those in midelevation and montane females (p < .01; Figure 2). In all variables, midelevation and montane birds showed insignificant differences (Tukey tests, p > .05). The results of smaller beaks in montane birds than lowland birds were in concordance with the expectation of Allen's rule.

| Characteristics of genomic data sets
We obtained a total of 17,654,632 SNPs from the 40 Taiwanese rufous-capped babblers. Each individual had an average of 5,473,846 ± 147,127 (SD) SNPs with an average sequencing depth of 15.09 ± 2.07× (Table S1). The three rufous-capped babbler individuals from mainland China had an average of 15,018,682 ± 25,259 SNPs and an average sequencing depth of 27.49 ± 0.31× (Table S1).
After including these additional samples, the 43 birds were segregated with a total of 33,714,823 SNPs.

| Historical demography of Taiwanese rufouscapped babblers
psmc analysis revealed contrasting effective population size (N e ) dynamic curves between Taiwanese and mainland Chinese rufouscapped babblers but similar curves within each group (Figure 1b). That is, Taiwanese populations maintained lower N e since 1-2 million years ago (Ma) compared to the mainland Chinese populations. msmc2 analysis revealed a similar pattern with lower N e in Taiwanese birds than in mainland Chinese birds, although less conspicuous mainland Chinese N e changes were inferred by this than the psmc analysis ( Figure S3). The fastsimcoal analysis inferred a divergence time between the Taiwanese and mainland Chinese populations at around 1 Ma, followed by a divergence event that separated populations in Western and Eastern Taiwan at 0.33 Ma (Figure 1c; Table S2). More recently, the divergence between montane and lowland populations occurred at least twice on Taiwan, once on each side of the Central Mountain Range (0.22 and 0.03 Ma for the Western and the Eastern population splits, respectively; Figure 1c; Table S2). fastsimcoal inferred stronger gene flow from mainland China to Taiwan than in the opposite direction (Figure 1c; Table S2). The four Taiwanese populations were connected to each other with more than 0.5 migrants per generation in both directions (Figure 1d; Table S2). Nevertheless, the PCA revealed marked population structure within Taiwan, although to a lesser extent between the lowland and montane populations in Western Taiwan (Figure 1e).

| Genetic variability within and between populations
Averaged over genome-wide 10-kb windows, the nucleotide diversity (π) of the four Taiwanese populations ranged from 0.0035 to 0.0037. Autosomes exhibited higher average π values (ranging from 0.0037 to 0.0039 in the four populations) than the Z chromosome (ranging from 0.0011 to 0.0013; Table S3). The four populations showed pairwise mean F ST ranging from 0.0195 to 0.0487.
Autosomes exhibited lower mean F ST (ranging from 0.0183 to 0.0460) than the Z chromosome (ranging from 0.0346 to 0.0845; Table S4).

| Genomic regions under altitudinally divergent selection
We identified a total of 311 genomic windows as candidate regions under altitudinally divergent selection based on ∆zF ST (cerulean dots in Figure 3). Among these candidate regions, 301 were located in F I G U R E 2 Altitudinal differences among rufous-capped babblers in their beak surface area (upper left), beak length (upper right), beak width (lower left) and beak height (lower right). All measurements are scaled by individual birds' wing lengths to obtain variables relative to body size and are rank-transformed to increase conformity to normality. Comparisons are visualized by violin plots, with dots and bars denoting means and standard deviations, respectively. For beak surface area, beak length and beak width, intergroup differences are assessed with male and female birds pooled because a pilot two-way ANOVA indicates an insignificant interaction between the sexual and altitudinal effects (p > .05). For beak height, individuals are separated by sexes given a significant sex × altitude interaction from the two-way ANOVA (p < .05). The plots show results from Tukey post hoc tests (***p < .001; **p < .01; *p < .05; blank: p > .05).
autosomes, showing mean F ST ranging from 0.1019 to 0.1265 and mean ∆zF ST ranging from 2.3346 to 2.7387 in the montane-lowland population contrasts (Tables S5 and S6). Ten candidate regions were located in the Z chromosome, showing mean F ST ranging from 0.1195 to 0.1664 and mean ∆zF ST ranging from 1.8520 to 2.3095 (Tables S5   and S6). Overall, the candidate regions were more concentrated in autosomes than in the Z chromosome (two-tailed Fisher's exact test, p < .01). We noted that for candidate regions containing both coding and noncoding DNA sequences, the montane-lowland F ST rose further in the noncoding region than in the coding one (two-tailed paired t test, p < .001).
Within the candidate regions, 320 SNPs were identified as fixed in either the montane populations or the lowland ones. Among the fixed SNPs, 174 (54.38%) were in the intergenic regions, and 143 (44.69%) and three (0.94%) were in the introns and the exons of protein-coding genes, respectively. Compared to the constitution of the candidate regions, including 52.9% (1,664,783 bp) of intergenic regions, 44.2% (1,374,851 bp) of introns and 2.9% (90,366 bp) of exons, the fixed SNPs were not more aggregated in intergenic than genic (i.e., intronic plus exonic) regions (two-tailed Fisher's exact test, p > .05). Nevertheless, fixed SNPs were over-represented in introns within genes (two-tailed Fisher's exact test, p < .05). In addition, fixed SNPs were less likely to occur in coding regions (i.e., exons) than in noncoding ones (introns plus intergenic regions; twotailed Fisher's exact test, p < .05).
The above results suggest that the altitudinal adaptation of Taiwanese rufous-capped babblers was predominantly achieved with noncoding genetic variants, probably associated with changes in regulatory regions that modulate gene expression. We then separated SNPs that were fixed in the montane populations (n = 190; blue cross in Figure 3) from those that were fixed in the lowland populations (n = 130; orange cross in Figure 3) and identified genes within either 250 kb upstream or downstream of them as their potential regulatory targets.
We identified 146 genes within 250 kb around the 190 montanefixed SNPs. Two of the 146 regulatory target genes-COL9A1 and SOX11-were probably related to avian beak morphogenesis as they or their paralogues were found to be associated with beak length of songbirds or involved in the facial bone development of chickens (Abzhanov et al., 2007;Bosse et al., 2017;Cheng et al., 2020). In addition, these 146 genes were enriched in several GO terms related to the gamma-aminobutyric acid (GABA) signalling pathway (Figure 4a) as well as in mouse phenotypes associated with survival at high elevations ( Figure S4a; see also Table S7 for related genes). For example, the phenotypes "decrease haemoglobin content," "decrease haematocrit" and "abnormal respiratory mucosa goblet cell morphology" were presumably related to high-elevation hypoxia adaptation, while "thick dermal layer" might be related to thermoregulation.
We found 187 genes within 250 kb around the 130 lowland-fixed SNPs. Similar to the regulatory target genes of the montane-fixed SNPs, these 187 genes were also enriched with GABA-related gene ontologies (Figure 4b). They were further enriched in the "decrease immunoglobulin G2a (IgG2a) level" mouse phenotype ( Figure S4b and Table S8). Neither genes surrounding montane-fixed SNPs nor those surrounding lowland-fixed SNPs were significantly overrepresented in KEGG pathways.

F I G U R E 3
Net genetic differentiation (∆zF ST ) between lowland and montane rufous-capped babbler populations, calculated along 10-kb nonoverlapping genomic windows. Plots are shown for two of the four cross-elevation population comparisons: L1-H1 and L2-H2. Assembled babbler scaffolds are arranged according to their synteny with zebra finch chromosomes. Cerulean dots indicate windows that show top 5% ∆zF ST values in each of the four cross-elevation population comparisons (separately selected for autosomal and Z-linked windows). Orange and blue plus signs above the panel denote windows containing SNPs that are (nearly) fixed in lowland and montane populations, respectively. Genomic windows that contain inversions inferred by delly analysis are shown in yellow.

| Chromosomal inversions
When examining altitudinally divergent genes located in the inversion regions, we focused on 45 high-confidential inversion loci (total length: 105 Mb, 10.2% in the genome; yellow dots in Figure 3) that were segregated over the 40 Taiwanese rufous-capped babblers. We used a two-tailed binomial test to evaluate whether the 45 inversion events tended to occur in autosomes or the Z F I G U R E 4 Enriched gene ontologies (GOs) of the genes that represent the putative regulatory targets of the SNPs that are fixed in either (a) the montane or (b) the lowland rufous-capped babbler populations. The regulatory targets are identified when genes have their start codons located within ±250 kb of the fixed SNPs. Results are shown for GO terms that each acquires a false discovery rate (FDR) of <0.05. chromosome, given the length ratio of these two chromosome types. We found an insignificant tendency regarding their occurrence in the two chromosome types (p > .05), and thus conducted the following analyses without separating autosomal and Z chromosomes. The candidate regions showed a significant tendency to be located in chromosomal inversions (59/311) compared to noncandidate regions (one-tailed Fisher's exact test, p < .001), suggesting that inversions contributed to the elevational divergence between populations. Concordant with this, in all four montanelowland population contrasts, the inversion-present candidate regions showed significantly higher D a than the inversion-absent candidate regions (one-tailed Wilcoxon rank sum tests, p < .05; Figure 5a). Similarly, in two of the four montane-lowland population contrasts, the inversion-present candidate regions showed significantly higher F ST compared with the inversion-absent candidate regions (one-tailed Wilcoxon rank sum tests, p < .05; Figure 5b). Note that the inversion-present and inversion-absent candidate regions showed insignificant differences in either D a or F ST between populations from similar elevations (one-tailed Wilcoxon rank sum tests, p > .05; Figure 5). That is, inversions did not ubiquitously increase interpopulation divergence-they only contributed to cross-elevation divergence.
We identified 24 genes that overlapped both the candidate regions and the inversions. Among these genes, CSDC2 was found to participate in hair follicle cycling in goats (Yang et al., 2017).
Deficiencies in RAMP3 and PSME3 had been shown to cause reduced body weight and retarded postnatal growth in mice, respectively (Dackor et al., 2007;Murata et al., 1999). In addition, RAMP3 plays a role in protecting hearts by reducing cardiac hypertrophy and perivascular fibrosis, and PSME3 affects the immune response by positively regulating NF-κB signalling when facing bacterial infection (Sun et al., 2016). DHX8, RPL27, MLX and POLR3H play roles in the F I G U R E 5 Comparisons of (a) D a and (b) F ST levels between altitudinally divergent windows that overlap chromosomal inversions and those that do not. Comparisons are made for population pairs from different elevations as well as for those from similar elevations and are visualized by violin plots with boxplots inside. Intergroup differences are evaluated by one-tailed Wilcoxon rank sum tests (***p < .001; *p < .05; ns: p > .05).

(a) (b)
regulation of transcription. GO enrichment analyses identified "respiratory basal cell differentiation" and "alveolar secondary septum development" (Figure 6a), both of which influence respiratory ability.
Finally, mouse phenotype enrichment analyses identified "abnormal retinal neuronal layer morphology," "abnormal hippocampus morphology" and "hyperactivity" (Figure 6b). In sum, genes located in the inversion regions could facilitate adaptations of rufous-capped babblers to different elevations by regulating their skin appendages, respiration, retinal morphology, immunity and metabolism.
We noted that montane and lowland birds showed little difference in their inversion allelic compositions: the 45 inversion loci showed differences in allelic frequencies only ranging from 0 to 0.12. In 10 of the 45 loci, montane and lowland birds were nearly fixed (allelic frequency ≥0.9 for each) for the same inversion alleles.
In the other 35 loci, neither of the two alternative alleles was dominant in the local environments. That is, great altitudinal divergence among populations in the inversion regions was not caused by differential usages of alternative inversion alleles between montane and lowland birds.

| DISCUSS ION
Ecological adaptation is one of the main mechanisms shaping biodiversity. High mountains have harsh environments that impose strong selective pressure on the organisms colonizing them, and thus serve as an ideal system for testing the effect of ecological adaptation in the early stage of population divergence. Based F I G U R E 6 (a) GO and (b) mouse phenotype enrichments of altitudinally divergent babbler genes that are located in chromosomal inversions. Results are shown for terms that each acquire a false discovery rate (FDR) of <0.05. on genomic sequences, we reconstructed the colonization history of rufous-capped babblers in Taiwan and suggest that elevationassociated selection has caused a divergence between the montane and lowland populations even with certain levels of gene flow between them. Consistent with Allen's rule, we found that the beaks of montane rufous-capped babblers were smaller than those of lowland ones. This implies that the bird has evolved smaller beaks to reduce heat loss in mountains, although alternative explanations cannot be completely ruled out. We further identified altitudinally divergent genes that may explain beak size variation (COL9A1 and SOX11) and other physiological changes in haematopoietic, metabolic, immune, auditory and vision systems, cerebrum morphology, and plumage development. These morphological and physiological changes may contribute to local adaptations to hypoxia, low temperature and the strong UV light at high elevations. We also found that altitudinally divergent genomic sites tended to occur in noncoding regulatory regions and chromosomal inversions in autosomes, rather than sex chromosomes, shedding light on the genomic architecture of ecological adaptation.

| Multiple high-elevation colonization
Our psmc and msmc2 analyses revealed divergent demographic trends between Taiwanese and mainland Chinese rufous-capped babblers starting around 1-2 Ma and a gentle N e decrease in Taiwan since then (Figure 1b; Figure S3). It is possible that the N e decrease reflects the period of accelerated mountain building on the island Our demographic analyses showed periods of sharp population growth in the recent past, which might have resulted from structured populations with certain levels of gene flow (Mazet et al., 2016).
Indeed, we found genetic differences between the four populations (H1, L1, H2 and L2) along with moderate or high levels of gene flow among one another. Given that both psmc and msmc2 assume panmixia, population structure and changes in gene flow may cause the spurious population increase (Mazet et al., 2016). Other factors, such as recent inbreeding, may also bias the inferences, leading to the recent artificial increase in N e (Lu et al., 2022). However, the mainland Chinese population was inferred with a relatively large N e . In addition, the habitats of rufous-capped babblers in Taiwan were likely to increase due to the birds' occupancy of new, montane niches. Thus, inbreeding might not have frequently occurred in the analysed populations of the recent past. Although it is difficult to ascertain the cause of the spurious N e change in psmc and msmc2, this does not obstruct our attempt to identify the genetic mechanism underlying the altitudinal adaptation of this bird.

| Regulatory regions contribute to elevational changes in the beak
In concordance with Allen's rule, we demonstrated that montane populations of Taiwanese rufous-capped babblers had smaller beaks relative to their body sizes than lowland populations (Figure 2). Note that factors other than thermoregulation may also cause the altitudinal change in this bird's beak size, such as the size of local prey.
For example, studies suggest that insects may have larger body size as environmental temperature increases (Maher & Shelomi, 2022;Ramírez-Delgado et al., 2016; but see Shelomi & Zeuss, 2017;Tseng et al., 2018), and rufous-capped babblers as insectivorous birds may need to have larger beaks to feed on larger insects in lower elevational areas with higher temperature. However, previous studies have shown similarly positive relationships between beak size and environmental temperature in other avian species with diverse diets (Fan et al., 2019;Laiolo & Rolando, 2001;McQueen et al., 2022;Symonds & Tattersall, 2010), suggesting Allen's rule is still a possible explanation for bird beak size. Nevertheless, the genes that are responsible for such temperature-driven beak size changes remain unclear.
Among our regulatory target genes, we propose that COL9A1 and SOX11 played roles in shaping rufous-capped babbler beak morphogenesis related to altitudinal adaptation because their paralogues were found to regulate avian beak size (Bosse et al., 2017;Cheng et al., 2020). Investigating genomic divergence between the longbeaked ground tit (Pseudopodoces humilis) and its short-beaked relatives, Cheng et al. (2020) showed that variation in beak length could be attributed to variation in COL27A1, which also showed signatures of positive selection in the ground tit. Bosse et al. (2017) found that individual SNPs in COL4A5 and SOX6 were significantly associated with the beak length of the great tit (Parus major). Examining cell differentiation during cranial bone formation in chicken embryos, Abzhanov et al. (2007) suggested the final fate of the skeletogenic progenitor cells was either chondrocytes or osteoblasts, partially determined by COL9A1, COL2A1 and SOX9. As a quantitative character, beak morphology is probably determined by multiple genes. Thus, COL9A1 and SOX11 may not be the genes solely responsible for the beak size change in rufous-capped babblers, but are promising candidate genes for further research, including functional validation. In addition, given the importance of noncoding genetic variation to elevational divergence in rufous-capped babblers, gene expression data from montane and lowland babblers could verify our conjecture.

| The function of regulatory target genes
Our results showed that the fixed SNPs were under-represented in coding regions. Therefore, the altitudinal adaptation of rufouscapped babblers may be dominantly fuelled by mutations in regulatory regions (or with regulatory functions), consistent with the findings for another songbird in Taiwan (Lai et al., 2019) and a raptor on the Qinghai-Tibet Plateau (Hu et al., 2022). Noncoding regions have been found to contribute more to ecological divergence at the initial stage of population divergence or speciation than coding regions (Jones et al., 2012). This is because noncoding regions have weaker pleiotropic effects and thus may contribute more to fast or early adaptive evolution than coding genes (Storz & Cheviron, 2016).
In addition, the only three fixed SNPs in the coding sequences found in this study were also synonymous mutations.
Although the potential regulatory target genes of the fixed SNPs did not include those previously found to be related to highelevation adaptation in other animals (e.g., Qu et al., 2013: HIF1AN in the great tit; Ma et al., 2019: EPAS1 in Tibetan pigs), our phenotype enrichment analyses showed enriched functions pertinent to hypoxia adaptation. Specifically, we found regulatory target genes associated with haemoglobin content/haematocrit regulation. This is consistent with studies that suggest when humans and other vertebrates are exposed to high elevations, their arterial haemoglobin concentration increases to counteract the reduction of arterial O 2 content (Siebenmann et al., 2017;Storz, 2021). We also found regulatory target genes functioning to modulate the level of IgG2a and pre-B cell number. In fact, the amounts of immunoglobulins (e.g., IgA, IgM and IgG) and B cells are altered in humans when they experience elevational change (Mishra & Ganju, 2010). Our findings echo previous studies which found selection signals in or varying copy numbers of immune genes across elevations in songbirds and chickens (Qu et al., 2013;Zhang et al., 2016). Xin et al. (2019) recently identified genes that were differentially expressed between yak (a high-elevation dweller) and cattle (from lower elevations) lungs and play a role in blood cell development as well as genes in their muscles affecting the animals' immune activities. Together with our study, these findings suggest the importance of haematopoietic and immunological functions in altitudinal adaptation.
Our GO enrichment analyses identified GABA-related terms ( Figure 4). It has been shown that GABA receptors in the brainstem nucleus tractus solitarius (NTS) neurons are critical for these neurons' function in modulating the ventilation of rats during hypoxia (Chung et al., 2006;Tabata et al., 2001;Tolstykh et al., 2004). Thus, the regulatory target genes associated with the GABA signalling pathway may mediate rufous-capped babblers' adaptation to mountains by changing their respiratory response to hypoxia. Notably, we also found several enriched phenotypes related to the auditory system, including "cochlear outer/inner hair cell innervation pattern," "spiral ligament fibrocyte," "cochlear labyrinth" and "cochlear duct morphology." The cochlea is a delicate organ, and insufficient oxygen supply could cause hearing damage. Indeed, it has been shown in mice that hypoxia treatment rendered reduced oxyhaemoglobin concentration in cochlear blood vessels (Reif et al., 2012). Another study demonstrated decreased spiral ganglion neurons and outer/ inner hair cells in the cochlea of rats after translocating the animals to high elevation for 3-6 months (Fan et al., 2016). It is possible that montane rufous-capped babblers evolved an adaptive mechanism to protect the cochlea system from hypoxia.
Together, these results suggest that altitudinal adaptation can be accomplished by changing nucleotides in the regulatory regions of target genes related to the haematopoietic system, immune response, GABA signalling pathway, auditory system (especially in the cochlearsystem) and beak morphology (see previous section).

| The function of genes in chromosomal inversions
Inversions that capture one or more locally favoured alleles in one genomic block may facilitate local adaptation by reducing gene flow between populations (Kirkpatrick & Barton, 2006). Given suppressed recombination in the genomic regions that contain inversions, divergent selection acting on individual genes within these regions would reduce gene flow over not just individual genic loci but the whole region, rendering large divergence blocks (Rieseberg, 2001).
Consequently, we investigated divergent genes between rufouscapped babblers' montane and lowland populations that were located in the inversions.
Notably, we identified CSDC2-a gene regulating the hair growth cycle in mammals, by means of downstream regulation of two key hair-follicle factors (Yang et al., 2017). Fang et al. (2022) found that the plumage of rufous-capped babblers diverged at different elevations. More specifically, montane populations have higher plumage UV-reflectance and brightness than lowland populations, suggesting plumage adaptation to high elevations that exhibit high environmental UV intensity and background brightness. Although it is unclear whether mammalian hair and avian plumage development share molecular mechanism, it is possible that CSDC2 played a role in the formation of this plumage, but this needs to be verified in future studies. Relatedly, the plumage represents an important agent in intraspecific communication (Maia et al., 2016). Therefore, the difference of UV-reflectance might also fuel adaptive change in the vision system of rufous-capped babblers. Indeed, divergent genes in the chromosomal inversions were enriched in mouse phenotypes related to retinal morphology ( Figure 6b). Alternatively, it has been shown in mammals that retinal morphology and optical neurons are affected when exposed to hypoxic conditions (Fischer et al., 2012;Sutherland et al., 2008).
Rufous-capped babblers may have had their vision systems modified by genes located in the inversions in response to different light environments or hypoxia at high elevations.
Additional enriched mouse phenotypes are related to cerebral morphology. This is probably because montane rufous-capped babblers have adapted to cope with hypoxia that may cause pathological changes in brains such as cerebral oedema (Wei et al., 2017;Zhang et al., 2013). Inversion-associated genes were also enriched for respiration-related functions (GO terms: "respiratory basal cell development" and "alveolar secondary septum development").
We found that neither inversion alleles nor reference alleles dominated the montane or lowland rufous-capped babbler populations. A potential reason for this is that the inversion and reference alleles in inversion regions both capture combinations of some genes favoured and some disfavoured by the local environments, especially when the environmental conditions have changed over time (Roesti et al., 2022). Thus, both the inversion and reference alleles can survive in the populations, stabilizing the inversions at an intermediate frequency and also maintaining the polymorphisms in the inversion regions (Kirkpatrick & Barton, 2006). In line with this explanation, we found significantly higher π values in the inversion-present candidate regions than in the inversion-absent ones (one-tailed Wilcoxon rank sum tests, p < .01; Figure S5). Alternatively, it is possible that when populations locally adapt to different environments in the face of gene flow, divergent selection may not drive fixation of inversions in the populations, but only increase genetic differentiation between inversion and reference alleles (Harringmeyer & Hoekstra, 2022).

| The sex chromosome plays a minor role in altitudinal adaptation
We found that altitudinally divergent genomic regions were significantly deficient in the Z chromosome, even though its overall level of divergence was higher than that of the autosomes. Regarding chromosomal inversions, we found that their occurrence was similar between autosomes and the Z chromosome. Therefore, the differentiation of altitudinally divergent genes between the montane and lowland rufous-capped babblers, which are at the early stage of population divergence, is not enhanced via accumulating mutations in the sex chromosome. Qvarnström and Bailey (2009) also suggest that ecologically divergent genes are mainly located in autosomal chromosomes at the early stage of population divergence, but accumulate in sex chromosomes towards the late stage of speciation. However, we cannot rule out the possibility that altitudinally divergent genes are associated with traits that are subject to sexual selection, such as social signal traits (Edwards et al., 2015;Malinsky et al., 2015), and thus restrict gene flow across elevations. For example, we found altitudinally divergent genomic regions associated with plumage development, vision and auditory systems that may be used for signal conveyance/perception during mate choice, although more experimental evidence is needed to confirm this.

| Implications to adaptation under climate change
Tropical mountains have remarkable biodiversity at a global scale (Orme et al., 2005), and thus altitudinal adaptation has been argued as an important driving force for speciation (Ribas et al., 2007;Smith et al., 2007). However, organisms specific to tropical montane habitats are also particularly sensitive to climate change (Şekercioğlu et al., 2012). Our work provides an insightful view into the role of genomic architecture in adaptation by demonstrating that mutations in regulatory regions and chromosomal inversions contribute to altitudinal adaptation in the initial stage of divergence more than mutations in coding regions. Given that noncoding regions have higher mutation rates, global warming may have less of an impact on traits that are regulated by such mutations in rufouscapped babblers as they may acquire new changes more quickly in response to rising temperatures. However, mutations trapped in chromosomal inversions face lower chances of recombination, and thus their rates of evolution would be lower and their associated traits would be more likely to make the organisms vulnerable to the changing climate.

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare. Acknowledgements. This research also generates genomic and ecological data that are accessible to the whole scientific community.