Genome resequencing and variant calling

The population structure and altitudinal colonization of the Rufous-capped Babblers in Taiwan were studied with genome-wide DNA polymorphisms generated by resequencing. We sampled individual birds from two low-altitude and two high-altitude populations (Fig. 1A), with blood collected from 10 males each population. We extracted genomic DNA using the Puregene Core Kit A (QIAGEN) following the manufacturer’s instructions. Paired-end sequencing libraries were prepared with the TruSeq DNA Sample Prep Kit v2 (Illumina) and then sequenced on an Illumina HiSeq 2500 platform with 2 × 150 base pair (bp) reads. Blood samples were collected from three additional male Rufous-capped Babblers from Hunan Province, mainland China (29.225911° N, 109.335892° E, 900 m asl) to better characterize the Taiwan-mainland relationship. We extracted genomic DNA from these additional samples using the SQ Tissue DNA Kit (OMEGA). Paired-end sequencing libraries were prepared with the MGIEasy Library Prep Kit V1.1 (MGI) and then sequenced on a BGISEQ-2000 platform with 2 × 151 bp reads. The above procedures resulted in a mean sequence coverage of 15.1 × (range = 12.2–19.6 ×) and 27.5 × (range = 27.2–27.8 ×) for the Taiwanese and mainland Chinese samples, respectively.

We removed adapter sequences in raw reads using the Trimmomatic 0.38 [74] commands ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:1:true. Trimmed reads were mapped to the Babbler draft genome [75] using BWA-MEM 0.7.17 [76]. To ensure that each unique DNA fragment was only counted once, we tagged duplicated reads using the MarkDuplicates tool in the Genome Analysis Toolkit (GATK) 4.0.6 [77]. To improve variant discovery, we conducted sample-wise base quality score recalibration (BQSR) prior to the final SNP calling, as described below. First, variants were identified using GATK HaplotypeCaller [78] and BCFtools mpileup [79]. Second, concordant variants were identified by GATK, with criteria recommended by the GATK team [80] adopted to filter lower-confidence variants (QD < 2, MQ < 40, FS > 60, SOR > 3, MQRankSum <  − 12.5, or ReadPosRankSum <  − 8). Third, BQSR was carried out with GATK BaseRecalibrator based on the retained variants. The above BQSR steps were performed twice before final variant calling by BCFtools mpileup, which resulted in a total of 33,714,823 biallelic SNPs that segregated over the 43 studied Babblers.

Relationships among populations and genetic structure estimation

We performed population genetic analyses based on putatively autosomal SNPs. For this, we mapped the Babbler’s assembled scaffolds [75] onto the (pseudo-)chromosomes of the Zebra Finch assembly (GenBank assembly accession: GCA_003957565.2) using Minimap2 [81]. Exclusion of the putative Z or W scaffolds yielded 31,964,249 putatively autosomal SNPs. We then randomly selected 10,000 SNPs to infer the inter-population phylogeny. We noted substantial divergence between Taiwan and mainland China (FST > 0.45), which rendered allelic fixation an indispensable factor in the drift processes. Accordingly, we adopted a method that takes allelic fixation into account for population tree reconstruction [55]. This method incorporates allelic fixation using a nonreversible model of genetic drift, which makes inferred evolutionary parameters “directional” and thus enables tree root identification without need to specify the outgroup.

To further assess the genetic distinctiveness among the four Taiwanese populations, we generated a second set of 10,000 randomly selected autosomal SNPs, each of which was polymorphic in the 40 Taiwanese individuals. With this data, we performed a principal component analysis (PCA) alongside an admixture model-based clustering analysis to investigate clustering of the 40 Taiwanese individuals. For the admixture analysis, we employed an ALStructure algorithm, which implements likelihood-free estimations to improve computational efficiency without compromising accuracy [82]. To confirm the robustness of the results, we repeated each of the phylogenetic analysis, PCA, and admixture analysis three times, each time with a newly generated random SNP set.

Reciprocal transplant experiment

Another 40 male adult birds were caught in 2016 from three Taiwanese populations—10 from each of the two lowland populations, L1 and L2 (< 300 m asl), and 20 from the montane population, H (2200–2600 m asl; Fig. 1A)—and used for a reciprocal transplant experiment. This experiment did not include the H′ population due to logistical difficulties for shipping birds quickly enough to our high-altitude common garden site (h) in another mountain (Fig. 1A). Half of the captured birds (n = 20) from each population were randomly assigned to a common garden (l) at a low-altitude research station (Jiji, 250 m asl) and the other half (n = 20) to a high-altitude one (h, Hehuanshan, 3,000 m asl) belonging to the Endemic Species Research Institute. The common gardens controlled for altitude-related environmental differences such as those in temperature and oxygen pressure. Individuals were acclimated to the low- or high-altitude condition for 35–94 days (between October, 2016, and January, 2017)—except for one individual, which was acclimated for 260 days (median = 64 days; Additional file 1: Table S1 and Fig. S2)—before being sacrificed. In the common gardens, the birds were kept in 1 × 1 × 1 m cages located in rooms with many windows kept open. This setup maintains the temperature close to the outdoor temperature as well as reduces exposure of the birds to bad weather and predators. To control for potential circadian and seasonal differences in expression plasticity, all birds were sacrificed between 10 AM and 2 PM on each sampling day from January 7–20, 2017 (the coolest month of a year; Additional file 1: Fig. S2), and the entire brain, liver, and other organs were removed from each individual within 7 min. Collected tissues were immediately placed in RNAlater (Invitrogen) and incubated at 4 °C overnight, followed by storage at − 20 °C until RNA extraction.

RNA sequencing and gene expression quantification

RNA extraction was carried out with the NucleoZOL Kit (MACHEREY–NAGEL). We identified four brain and one liver samples that showed RNA integrity number (RIN) values of < 8 in Bioanalyzer 2100 (Agilent) and excluded them from subsequent analyses. Libraries for polyA-enriched transcriptomes were constructed using TruSeq RNA Library Prep Kit v2 (Illumina); they were then sequenced on Illumina HiSeq 2500 for 2 × 125 bp paired-end reads. We sequenced 36 brain and 37 liver samples to yield 19,849,176 to 25,221,037 (mean = 22,554,647) read pairs per sample. We used the same Trimmomatic commands as for DNA resequencing to remove adapter sequences, and the command MINLEN:40 to eliminate short reads.

Two additional liver samples were sequenced in a second batch, with 2 × 201 bp reads generated on the same sequencing platform described above. From the resulting 40,139,534 and 45,645,423 read pairs, we randomly selected 22,000,000 per sample with seqtk [83] for downstream analyses to maintain similar sequencing depths across sequencing batches. When trimming these additional samples, we added the command CROP:125 in Trimmomatic to make each read ≤ 125 bp long. We showed in a sample clustering dendrogram (see the next section for details) that the two additional samples were subject to little batch effect.

For all samples, we aligned trimmed reads to the Babbler genome [75] using the splice-aware aligner HISAT2 2.1.0 [84]. When building the genome index for alignment, we incorporated splice site and exon annotations using python scripts hisat2_extract_splice_sites.py and hisat2_extract_exons.py from the HISAT2 package. We then quantified gene expression levels by counting the numbers of fragments (i.e., sequences each bookended by a pair of reads) mapped to exons of genes using featureCounts [85] with the default settings.

Prefiltering of genes subject to confounding effects

We examined clustering of samples based on their transcriptomic profiles to identify main factors affecting gene expressions. To this end, we used the build-in method of the R package DESeq2 [86] for sample normalization. We then summarized inter-sample expression differences using Euclidean distances and performed an average-linkage analysis for sample clustering. Prior to Euclidean distance computations, we applied a variance-stabilizing data transformation [87] so that genes with contrasting expression levels would contribute approximately equally to sample clustering. The result revealed samples clustered by tissue types (Additional file 1: Fig. S3). In addition, the two liver samples sequenced in the second batch did not form their own cluster, but were well grouped with the other liver samples, indicating little batch effect.

The unequal lengths in individual birds’ acclimation time could tilt gene expressions and thus confound changes caused by factors of interest (described in the next section). We used likelihood-ratio tests implemented in DESeq2 to identify genes affected by such non-identical acclimation durations. Provided that the effect could differ when acclimating an individual to a native or a non-native environment, we identified genes under four conditions separately: lowland birds (pooled samples of L1 and L2) acclimated to low altitude, lowland birds acclimated to high altitude, montane birds (samples of H) acclimated to low altitude, and montane birds acclimated to high altitude. We binned acclimation durations into intervals of > 30, > 60, and > 90 days, and used the likelihood-ratio tests to compare models with and without including the acclimation duration as a covariate. We used the independent filtering step of DESeq2 [88] to enhance detection power and applied multiple testing corrections to control the false discovery rate at < 0.05. We then excluded genes identified in any of the four test conditions from downstream analyses for the corresponding tissues (the liver or the brain).

Identification of adaptation-associated genes

We identified genes associated with the Rufous-capped Babbler’s high-altitude adaptation as those fulfilling two conditions: (1) genes that exhibited large and directionally concordant expression differences in the H vs. L1 and H vs. L2 population contrasts when samples were acclimated to the high-altitude garden with the montane environment, where high-altitude adaptation occurred, and (2) genes that exhibited small expression differences in the L1 vs. L2 population contrast when samples were acclimated to the low-altitude garden with the lowland environment, to which both low-altitude populations have already adapted. Given that remarkable gene expression differences between populations might represent divergence irrelevant to altitudinal adaptation or resulting from genetic drift, we incorporated the second condition to penalize genes that exhibited non-altitudinal divergence. Genes that fulfilled both conditions, referred to as the “altitudinally concordant differential expression (ACDE)” genes, were identified based on the π-values [57]. Given π =|log-fold change|× –log10(P-value), the π-values were always non-negative and increased as increasing magnitude and statistical significance of between-group expression differences. We derived π-values based on log-fold change and the P-value estimates from DESeq2. For each gene, we then calculated ∆π1 = π(H vs. L1) – π(L1 vs. L2) and ∆π2 = π(H vs. L2) – π(L1 vs. L2). We identified the ACDE genes as those at the top 5% for both ∆π1 and ∆π2, and meanwhile showing the same regulation directions (i.e., + /– sign) in the H vs. L1 and H vs. L2 contrasts.

Gene function annotations and enrichment analyses

Babbler genes were functionally annotated based on their orthologous relationships with genes of Chicken, Duck, Collared Flycatcher, Turkey, and Zebra Finch. To this end, we acquired the babbler proteome predicted from the genome assembly of this species [75] alongside proteomes of the other five avian species from Ensembl v95 [89]. A longest peptide isoform per gene was kept for each species, with which multi-taxon gene orthology was inferred using the reciprocally blasting and gene clustering algorithms in SonicParanoid 1.3.0 [90]. Based on the orthology, we associated Babbler genes with gene ontology (GO) annotations of the other five species acquired from Ensembl v95 and merged as many non-replicated annotations from the latter five species as possible onto Babbler genes. Likewise, we associated Babbler genes with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations of the other five birds released in June 2015.

For each tissue type studied, we examined enrichment of ACDE genes in specific GO or KEGG terms compared to the other genes expressed in the tissue, with statistical significance evaluated by one-sided Fisher exact tests. Multiple testing corrections that control false positive rates were applied under the Biological Process, Cellular Component, and Molecular Function GO main categories, separately. Within each main category, corrections were applied in a level-wise manner along the ancestor-offspring hierarchical relationships of individual terms. One false positive rate control was applied to all KEGG terms, which lacked hierarchical structure. We performed these above enrichment analyses using customized Perl scripts.

Expression plasticity estimation for adaptation-associated genes

To infer reinforcing or reversing plasticity for each ACDE gene, we identified each gene’s directions in two expression change quantities—PC and GC (visualized in Fig. 2A). PC and GC measure the plastic response of ancestral populations and the genetic divergence between ancestral and descendant populations, respectively. We identified directions of these quantities by contrasting samples from the Rufous-capped Babbler’s extant populations (shown in Fig. 2A). Reinforcing and reversing plasticity were then inferred when obtaining consistent and opposite PC and GC directions, respectively.

Note that a gene’s expression level in the low-altitude population was used to determine both of PC and GC. Mallard et al. [58] recently demonstrated that reversing plasticity could artificially appear more prevalent than reinforcing plasticity when random sampling errors associated with the shared element falsely led to both PC ≠ 0 and GC ≠ 0. To correct for such a putative artifact, we adapted Ho and Zhang’s [24] parametric bootstrap method to our cases, which aimed to identify genes with PC ≠ 0 and GC ≠ 0 resulting from genuine differences rather than random sampling errors and was implemented as follows. For each gene, we generated two Gaussian distributions to simulate its PC and GC, respectively. The PC and GC Gaussian distributions each had a mean and a standard deviation equal to the empirical log-fold change and the associated standard error, respectively, of the PC or GC estimated by DESeq2. From the two Gaussian distributions, we drew PC and GC random samples iteratively (i.e., bootstrap replicates), with reinforcing or reversing identified by a GC-PC sample pair each time. We performed 1000 such bootstrap replicates per gene, and we concluded bootstrap-supported reinforcing or reversing once either plasticity type was obtained in ≥ 950 replicates. We conducted Gaussian distribution generations and subsequent random sampling using the rnorm function of the base R 3.6.0. We evaluated unequal frequencies between reinforcing and reversing plasticity using two-sided binomial tests against a null proportion of 0.5.

To investigate whether ACDE genes more frequently had their plasticity persisting or evolved during the bird’s high-altitude adaptation, we obtained the measure of a third expression change quantity—GCb—based on the population contrast shown in Fig. 3A. A binary classification between plasticity persistence and plasticity evolution was then carried out by comparing GCb against GC: plasticity persistence was inferred with a GCb-to-GC ratio between 0.5 and < 1.5 while plasticity evolution was inferred with a ratio between − 0.5 and < 0.5 (visualized in Fig. 3B). We used two-sided binomial tests to evaluate unequal prevalence between plasticity persistence and plasticity evolution in the ACDE genes.

Examining factors associated with evolution of expression plasticity

To examine the factors that were associated with plasticity evolution, we first used the two-sided Fisher exact test to evaluate unequal frequencies of plasticity evolution between ACDE genes exhibiting reinforcing and reversing plasticity. Secondly, to test another hypothesis that the ancestral plasticity magnitude (|PC|) determined whether the plasticity evolved, we used the Kruskal–Wallis test to assess the difference in |PC| between groups with evolved and persisting plasticity.

To dissect the above two patterns in-depth, we firstly adopted two other categorization schemes of the ACDE genes regarding the extent to which their plasticity evolved. Specifically, we used three and four equal-interval bins, respectively, to group genes according to their degree of plasticity evolution, with the gradient extending from a GCb-to-GC ratio ~ 1 (low degree of plasticity evolution) to a ratio ~ 0 (high degree of plasticity evolution). In the three-binned categorization, we grouped genes by GCb-to-GC ratios between 0.75 and < 1.25, between 0.25 and < 0.75, and between − 0.25 and < 0.25. In the four-binned categorization, we grouped genes by GCb-to-GC ratios between 0.83 and < 1.17, between 0.50 and < 0.83, between 0.17 and < 0.50, and between − 0.17 and < 0.17 (both categorizations are visualized in Additional file 1: Fig. S4A). Secondly, we tested the two focal patterns with genes’ plasticity evolution measured on a continuous scale. To this end, we used DESeq2 to quantify the absolute difference between the lowland and montane populations in their reaction norms, estimated by the population × acclimation environment interaction in the linear model (visualized in Additional file 1: Fig. S5A).

We also examined whether the hypothesis that |PC| determined the evolution of plasticity held in the rest of the Rufous-capped Babbler’s genes as in the ACDE genes. We pooled samples from the two low-altitude populations for acquiring |PC| and the above plasticity evolution measures for each gene. The sample pooling was to be consistent with the expression module analyses described in the following section.

We used the 2 × C Fisher exact test to evaluate inequality among plasticity evolution categories in their proportions exhibiting reinforcing and reversing plasticity (C = 3 and 4 in cases of three and four evolution categories, respectively). This was followed with two-sided 2 × 2 Fisher exact tests for post hoc comparisons between pairwise categories. Similarly, we used the Kruskal–Wallis test to evaluate inequality among evolution categories in |PC|, followed with post hoc pairwise comparisons using two-sided Dunn tests. We implemented Fisher exact tests in the R package rstatix [91], Kruskal–Wallis tests in the base R, and Dunn tests in the R package FSA [92]. Multiple testing corrections were applied to post hoc Fisher and Dunn tests. In cases where plasticity evolution was scaled continuously, we used the Kruskal–Wallis test to evaluate the difference between reinforcing and reversing plasticity genes in their plasticity evolution; we used the one-sided Spearman correlation test in the base R to examine the relationship between |PC| and plasticity evolution.

Transcriptome-wide co-expression analyses

We demonstrated the association between |PC| and plasticity evolution to be transcriptome-wide rather than limited to the ACDE genes (see “Results”). We hypothesized both |PC| and plasticity evolution to be dependent on the level of gene–gene interactions in the ancestral population, leading to the observed positive relationship between |PC| and plasticity evolution. To test this hypothesis, we first applied the WGCNA analyses [61, 62] to genes of the whole transcriptomes to delimit groups of expressionally interacting genes (co-expression modules) in the Rufous-capped Babbler’s lowland populations L1 and L2. We pooled samples from the two low-altitude populations to fulfill the minimum sample size for this analysis [93]. Briefly, the WGCNA grouped genes based on their correlations in expressions. By raising the absolute values of the correlations to a power β > 1 (referred to as “gene adjacency”), the analysis gives weight to gene clustering with high correlations. Prior to the analysis, we filtered out genes with many zero RNA fragment counts because such genes would show spuriously high correlations with one another; specifically, we removed genes showing a median absolute deviation value = 0. We quantified correlations between genes using the biweight midcorrelation, which is robust to the presence of outlier samples [94]; we set 0.1 as the maximum percentile of data that can be considered outliers. We performed “signed” WGCNA analyses such that genes were clustered based on positive correlations. We chose power β values under the following considerations: (1) the fit to the scale-free topology model > 0.8; (2) a large β value to avoid clustering based on spurious correlations between genes; (3) clustering supported by a reasonable mean gene connectivity value (~ 100), which necessarily decreases with increasing β. Consequently, we selected β = 12 and 18 for the brain- and liver-based WGCNA analyses (Additional file 1: Fig. S13A). We performed average-linkage clustering of the genes based on their topological overlaps derived from the adjacency measures. We then applied the dynamic tree cut algorithm for module delimitations with the minimum module size and the deep-split parameters set to 30 and 2, respectively.

To test the dependence of plasticity evolution on gene–gene interactions, we conducted three examinations. Firstly, we used the one-sided Spearman correlation test to evaluate correlations between the module size and the proportion of genes with plasticity evolution in the module. However, genes may differ in the levels of gene–gene interactions even when they occur in modules with similar sizes. Considering this possible scenario, we obtained measures of the intra-modular connectivity (kIM) of individual genes, calculated for each gene as the sum of its adjacencies. For the second examination, we subdivided genes’ plasticity evolution into two, three or four levels and then used Kruskal–Wallis tests to evaluate inequalities among evolution categories in the kIM values. In cases with three or four evolution categories, Dunn tests with multiple testing corrections were used for post hoc pairwise comparisons. For the final examination, we measured genes’ plasticity evolution on a continuous scale and used the one-sided Spearman correlation test to evaluate the hypothesized negative correlation between kIM and the plasticity evolution. In the second and third examinations, we adopted both raw and scaled kIM. We scaled kIM by its maximum value per module to render values all between zero and one; we aimed to confirm the resultant patterns insensitive to different kIM scales across modules of various sizes.

ACDE vs. non-ACDE genes

We used Kruskal–Wallis tests to evaluate differences between ACDE and non-ACDE genes in (1) kIM, (2) ancestral plasticity magnitude (= plasticity magnitude in the lowland populations), and (3) the extent to which plasticity evolved. For (2) and (3), we derived focal quantities with the two lowland populations L1 and L2 pooled. We used the continuous measure of genes’ plasticity evolution for (3).

Statistical significance

Results from all statistical tests were regarded significant whenever obtaining nominal P-values (or multiple testing adjusted P-values) < 0.05. We applied Benjamini and Hochberg’s [95] procedure for multiple testing corrections throughout this study.