Rapid adaptive radiation of Darwin’s finches depends on ancestral genetic modules

Recent adaptive radiations are models for investigating mechanisms contributing to the evolution of biodiversity. An unresolved question is the relative importance of new mutations, ancestral variants, and introgressive hybridization for phenotypic evolution and speciation. Here, we address this issue using Darwin’s finches and investigate the genomic architecture underlying their phenotypic diversity. Admixture mapping for beak and body size in the small, medium, and large ground finches revealed 28 loci showing strong genetic differentiation. These loci represent ancestral haplotype blocks with origins predating speciation events during the Darwin’s finch radiation. Genes expressed in the developing beak are overrepresented in these genomic regions. Ancestral haplotypes constitute genetic modules for selection and act as key determinants of the unusual phenotypic diversity of Darwin’s finches. Such ancestral haplotype blocks can be critical for how species adapt to environmental variability and change.

. Incomplete separation can be seen using genome-wide SNPs, while the three species show separation using the regions of association. Points are colored by species and each point is an individual, which were collected from a variety of islands (Data S1) (c) Average beak length and average beak width for all Geospiza species and island combination (from Grant et al., 1985 (69)) listed in Data S1 and used in comparative genomics analysis. Each point is colored by its species designation, and there are multiple points for species we have sequenced on multiple islands. Selected species discussed in the manuscript are labelled. Note the five species with the largest beak dimensions, also highlighted in Fig. 2, and the extensive phenotypic variation found in G. fortis. For each region, bounded by dashed line in the plots, we calculated Z-transformed FST, between population nucleotide diversity (dxy), nucleotide diversity (π), and Tajima's D. Each possible pairwise population contrast is shown and the genes that fall into the region of association are marked in the upper most plot. On the right, boxplots illustrate the difference between the statistic calculated in the main plot minus the median value across all regions of the genome that do not fall into regions of associations marked in admixture mapping. Values that fall above or below zero are more extreme values for that statistic compared to background values. In the boxplot, the centerline indicates the median, bounded by the 25 th and 75 th percentile, with whiskers extending to 1.5x the interquartile range. Fig. S5. Manhattan plots for admixture mapping of body and beak size variation in tree and ground finches. Top, admixture mapping comparing small, medium, and large tree finches (C. parvulus, C. pauper, and C. psittacula, respectively), classified as 0,1,2, based on beak and body size. Points in red mark those that surpass the significance threshold in the Geospiza admixture mapping results (Fig. 1b). Below, the same admixture mapping result as in Fig. 1b, for all SNPs recovered in tree finches, and the points surpassing the significance threshold set by permutation (see Methods) marked in red. The position of the locus (locus #3) containing HMGA2 is marked in both plots.   Estimated time since divergence, based on between-population nucleotide diversity (dxy), is shown for all three pairwise contrasts at each of the 28 loci involving G. magnirostris, G. fuliginosa, and Certhidea olivacea. The black horizontal line represents the time to divergence between G. magnirostris, and Certhidea olivacea as previously reported (21). Error bars represent 95% confidence intervals.

Fig. S9
. Allele age distribution at top associated SNPs. Allele age estimates for the top 100 associated SNPs at each of the 28 loci. Allele ages were estimated in the software GEVA (70), which uses a Hidden Markov Model approach for unbiased estimates of allele age at each position. The average estimate for time to most recent common ancestor (in generations) for each allele was estimated using both the mutation and recombination rate. Time in generations was converted to allele age using an estimated generation time of 5 years. On the bottom, the "Random" distribution is of 889 variants sampled across the genome every 10,000 variants. Peaks are sorted by the mean of all variants in the locus and beneath each peak is a rug plot of all positions at that locus. A vertical line in the distribution marks the 97.5% percentile cut-off and each distribution is colored by the time estimate. Note the multimodal distribution of several loci, likely representing both old and recent mutations. On the right, the 97.5% percentile cut-off for time to most recent common ancestor is listed for each locus, the average of which is presented in the main text. This value represents the estimated upper age estimate of haplotypes.  S10. Enrichment analysis. The software GREAT (41) was used to screen for enrichment of gene annotation terms associated with the 28 differentiated regions. Annotation term sharing between significantly enriched terms from Mouse Genome Informatics (MGI) gene expression databases are shown as a heat-map indicating the degree of gene sharing between each significant category. Gene sharing was defined for category pairs, as the proportion of genes from the category with the fewest genes which also occur in the category with the highest number of genes. Annotation terms along the x-axis occur in the same order as on the y-axis and the labels on the y-axis indicate term names. Also shown for each term are -log(Bonferroni P-value), degree of enrichment (E) and the number of genes in each category (N).   Table S1: List of all tissues used for RNA-seq for genome annotation. The type of library preparation and kit are listed, as well as the accession numbers for raw data on the European Nucleotide Archive (ENA) database. (Excel file not included in pdf) Table S2: Summary of the number of samples by species and island re-sequenced and used for population genetic analysis. For a complete list by sample, see Data S1.
Watterson's theta estimate of 0.00126 (21). Last, we used the pade command to generate Padé coefficient tables using the same estimate of Watterson's theta and the suggested 11 coefficients. We ran LDhelmet in 50 bp windows for 10,000,000 iterations (discarding the first 200,000 as burn-in), and the following additional flags: --max_lk_end 100 --prior_rate 0.05 -b 5.0 A previous study identified a block penalty of 5, as used here, as the ideal value for demonstrating fine scale recombination patterns (35). We extracted ρ/bp, representing the population scaled recombination rate, using post_to_text within LDhelmet. We summarized ρ by 20 kb sliding, non-overlapping windows (using the R package windowscanr, https://github.com/tavareshugo/WindowScanR), and report ρ/kb in comparison with diversity statistics. We additionally converted ρ to recombination fraction following r = ρ/2cNe (68) and Ne = 33,861 (21). Recombination fractions were converted to cM/Mb using the Morgan function r * 100 (for all plots in Fig. S2).
Genotype phasing. We inferred haplotypes using a combination of WhatsHap (82) (v0.18) and SHAPEIT4 (83) (v4.1.3) to generate a phased variant callset of biallelic SNPs. We used WhatsHap under default settings to identify phase informative reads (i.e., a read spanning at least two heterozygous sites) and the resulting VCF file was used as input for SHAPEIT4. WhatsHap was run separately for each individual, after which all individual VCFs were merged using the bcftools merge command. We ran SHAPEIT4 with a constant recombination rate of 1 cM per Mb and an expected error rate in the phase informative read sets of 0.0001 (the default setting).
Admixture mapping. We used GEMMA (37) (genome-wide efficient mixed-model association 0.98.1, 2018-12-10) to search for loci associated with the scaling up in beak and body size from G. fuliginosa (n = 26), G. fortis (n = 34), and G. magnirostris (n = 15) by recoding each species as 0,1,2, respectively, and estimating the strength of association for each SNP under an additive model. In order to account for relatedness among samples, we first calculated a relatedness matrix. As input for the relatedness matrix, we pruned autosomal SNPs by retaining only 1 randomly selected SNP per 20 kb in bcftools +prune (v1.11-54-gaf54707) to reduce the effect of linkage disequilibrium on the relatedness calculation. We converted the resulting pruned VCF file to bimbam_dosage format using qctool (https://www.well.ox.ac.uk/~gav/qctool_v2/) and calculated the relatedness matrix using the GEMMA flag -gk 1. A separate relatedness matrix was calculated for SNPs located on chromosome Z using identical parameters. We ran linear mixed-effects models (LMM, -lmm 2 option in GEMMA) with a maximum likelihood estimation (MLE) using the relatedness matrix calculated previously as a covariate and separately for all autosomes and chromosome Z. The same analysis was run for small, medium and large tree finches (C. parvulus n = 25, C. pauper n = 10, and C. psittacula n = 11). For Geospiza, we used the code implemented in gemma-wrapper to run permutations on 100 independent runs of shuffled phenotype values to establish a baseline significance threshold for SNPs exceeding the 67 th percentile of the distribution of P-values. This resulted in thresholds of -log10(P-value) = 7.3 (P-value < 5x10 -8 ). Using custom R scripts (v4.0.3) (84), we identified genomic regions that contain SNPs exceeding this threshold by extending peaks from the top SNP of association down to an arbitrary log10(P-value) threshold of 3.6. Within these regions, we merged all SNPs within 75 kb into a single genomic interval using the GenomicRanges package and requiring each interval to have more than one SNP. This resulted in 28 genomic regions of association on autosomes that exceeded the threshold set by permutation. In order to extract SNPs tagging the haplotypes of association, we selected the top SNPs of each region use the `slice` function in the R package dplyr to choose the top 100 SNPs based on their strength of association. In some genomic regions, many SNPs are linked and those with identical P-values were retained for analyses including the top SNPs, meaning that in 15 genomic regions shown in Data S5 include > 100 SNPs. In four regions, <100 SNPs exceeded the thresholds defined here.
Genetic diversity and divergence analyses. We calculated unbiased estimates of FST, dxy and π using the software pixy (v1.0.0.beta1) (85). For this analysis, we used the VCF that included both high confidence invariant sites and biallelic SNPs (see methods under Shortread variant analysis). We calculated FST and dxy for all three possible combinations of G. fuliginosa (n = 26), G. fortis (n = 34) and G. magnirostris (n = 15) and nucleotide diversity (π) for each species in non-overlapping 20 kb windows and by applying a minor-allele frequency filter of 0.05 within the pixy command module. We additionally removed windows with < 100 SNPs. Time to divergence was estimated by the conversion of genetic divergence (dxy) to time using T = dxy/(2μ) and a mutation rate of 2.04 X 10 -9 per site per year (21).
Tajima's D was calculated using the same input for 20 kb genomic windows using vcftools (v0.1.16) (81). In order to overlap divergence statistics, which were calculated in windows, with per-SNP admixture mapping results, we used the "windowscanr" R package (https://github.com/tavareshugo/WindowScanR). windowscanr: Apply functions using sliding windows. R package version 0.1.) to calculate mean -log10(P-value) in 20 kb windows and overlap these values with the divergence statistics using the "fuzzyjoin" R package (https://github.com/dgrtwo/fuzzyjoin) with a max distance of 10 bp. We calculated correlation coefficients between FST and -log10(P-value) using the base R command "cor" and report R 2 as the square of this output. We calculated R 2 between recombination rate (r) and nucleotide diversity using the same method.
Allele age estimation. We estimated the age of the top 100 SNPs in each locus of association using the probabilistic method implemented in the software GEVA (v1beta) (70). GEVA estimates the age of individual variants by generating a probability distribution of the time to most recent common ancestor for all haplotypes in a dataset using Hidden Markov Models.
The strength of this approach is that it incorporates both mutation rate and recombination rate into the estimate of allele age and does not make assumptions about historical demography and selection. We were interested in the age of the haplotypes that segregate amongst G. fuliginosa, G. fortis, and G. magnirostris and first filtered the phased variant files generated by SHAPEIT4 for these 75 individuals. We next converted this variant file and the recombination map (described above) to the GEVA input binary format using the conversion function provided by the software. In doing so, we assumed that the alternate allele called relative to the C. parvulus reference genome represents the derived state. We ran GEVA for the top 100 SNPs in each locus separately with a mutation rate of 1.02 X 10 -8 per site per generation (21), 100 max concordant and 100 discordant pairs, and an estimated effective population size of 40,000. We identified this effective population size as the approximate ancestral effective population size based on previous ancestral effective population size reconstructions (21). To compare allele age estimates to the genomic background, we excluded the genomic regions containing the 28 loci and estimated allele age for every 10,000 SNPs across the genome (n = 889 positions). We estimated the upper limit of allele age at each locus using the quantile function in R.

Analysis of introgression.
We first used the software Dsuite v0.4 (86) to calculate ABBA and BABA sites for all species and islands of G. magnirostris, G. propinqua, and G. scandens using T. bicolor and L. noctis as outgroups to determine the derived state for each allele. We selected these three ingroup taxa, because G. scandens and G. propinqua share a most recent common ancestor more recently than either does to G. magnirostris (21). We first established that G. propinqua, which is only found on Genovesa and sister to the allopatric common cactus finch G. scandens, shares the largest number of derived alleles with the population of G. scandens on Pinta. We subsequently used G. scandens samples from Pinta as P1 in a phylogenetic arrangement as shown below: We next used Dsuite to calculate the fraction of introgression in sliding windows of 50 SNPs (by 25 SNPs) across the genome as measured by the df which is distinct from other measures of introgression in that it accounts for divergence between each comparison in the trio, to better assess the difference between ILS and introgression (39).
We additionally compared dxy values in the 28 regions between G. propinqua and G. magnirostris using the same methodology in Pixy as described above.
Gene enrichment analysis. The software GREAT v3.0 (41) was used for annotation term enrichment analysis. Significant loci from admixture mapping were converted to chicken genome assembly (Galgal6, NCBI accession GCA_000002315.5) coordinates using Progressive Cactus v0.1-c4bed56(87) to first perform whole genome alignment between the C. parvulus genome assembly and Galgal6. Next, halLiftover from the halTools package (88) was used to obtain the corresponding chicken coordinates. We then used the UCSC chicken to human chain file to extract syntenic human (hg19) coordinates. Extracted human coordinates were merged and then submitted to the online GREAT tool v3.0 (http://great.stanford.edu/public/html/) using standard settings. Enriched terms (Bonferroni Pvalues < 0.05) from hypergeometric tests were retained and, as a conservative measure, only terms showing enrichment values ≥2 were retained for further analyses and plots.
Phylogenetic reconstruction. Neighbor-joining (NJ) trees were generated by first calculating per-species average allele frequencies from genotypes recorded in VCF files, then reformatting allele frequencies to serve as input for PopTree (89), where NJ trees were generated using the "poptree" command with FST as the distance option and 100 bootstrap iterations. For the genome-wide species tree the input VCF was filtered to include all SNPs with a minor allele frequency ≥0. 15. The topology of this tree is similar to previous studies (21). In order to create a cladogram for display in Fig. 3A, we built a chronogram for the species tree using the phytools commands makeChronosCalib and chronos, setting a maximum age at 1 million years following Lamichhaney et al. (21). The species ages from this tree were not used for analysis in this manuscript and we refer the reader to Lamichhaney et al. (21) for a more comprehensive species tree analysis, although we note the estimated node ages are comparable to our estimates from the NJ species tree.
Separately, we created a maximum likelihood phylogenetic tree for all G. fuliginosa, G. fortis, and G. magnirostris samples included in the admixture mapping analysis (Fig.  S3a). For this analysis, we used the pruned genotype data that was used for generating the relatedness matrix in GEMMA (see Methods above). We used the BEDtools (v2.29.2) (90) intersect command with the `-v` flag (while supplying a bed file of genomic regions) to retain only genomic regions laying outside the 28 loci identified in admixture mapping and converted the VCF file to fasta format using vcf2phylip v2.4 (91). We included a randomly selected allele at each heterozygous position by supplying the -resolve-IUPAC flag in vcf2phylip. We ran FastTree (92) to generate a maximum likelihood phylogenetic reconstruction with the following settings: FastTree -nt -log finch_fasttree_log_file -gtr -fastest finches_geospiza.fasta We used plink2 (93) (v2.00a2.3LM) to generate a covariance matrix of all genotypes either within or outside of the 28 genomic regions identified in the admixture mapping analysis. As input we used the pruned dataset described above that was used to generate the maximum likelihood phylogeny of the three Geospiza species. We also created a genotype matrix (VCF file) of only the 28 regions using a bed file of these regions supplied to BEDTools intersect. We ran the plink -pca command under default settings for both the genotype data within and outside the 28 regions of interest, and plotted the resulting covariance matrix in R (Fig. S3b).

Analyses of allele frequencies.
Per-species average allele frequencies for 18 species of Darwin's finches, the Big Bird lineage and two outgroup species from genotypes recorded in VCF files were calculated using custom python scripts.
Topology weighting. We used topology weighting to compare local genealogies across the genome to those in regions we identified using admixture mapping. For simplicity, we limited our analysis to 4 taxa (P. inornata, C. parvulus, G. magnirostris, and G. fuliginosa). We created neighbor-joining trees in non-overlapping 100 SNP windows using phyml v.3.3.2 (94) and calculated topology weightings using twisst (95).          Fig. 11 a b T. guttata RUNX2 T. guttata ALX1 Supplemental Fig. 12