An island- hopping bird reveals how founder events shape genome- wide divergence

When populations colonize new areas, both strong selection and strong drift can be experienced due to novel environments and small founding populations, respectively. Empirical studies have predominantly focused on the phenotype when assessing the role of selection, and limited neutral- loci when assessing founder- induced loss of diversity. Consequently, the extent to which processes interact to influence evolutionary trajectories is difficult to assess. Genomic- level approaches provide the opportunity to simultaneously consider these processes. Here, we examine the roles of selection and drift in shaping genomic diversity and divergence in historically documented sequential island colonizations by the silvereye ( Zosterops lateralis ). We provide the first empirical demonstration of the rapid appearance of highly diverged genomic regions following population founding, the position of which are highly idiosyncratic. As these regions rarely contained loci putatively under selection, it is most likely that these differences arise via the stochastic nature of the founding process. However, selection is required to explain rapid evolution of larger body size in insular silvereyes. Reconciling our genomic data with these phenotypic patterns suggests there may be many genomic routes to the island phenotype, which vary across populations. Finally, we show that accelerated divergence associated with multiple founding steps is the product of genome- wide rather than localized differences, and that diversity erodes due to loss of rare alleles. However, even multiple founder events do not result in divergence and diversity levels seen in evolutionary older subspecies, and therefore do not provide a shortcut to speciation as proposed by founder- effect speciation models.

alleles (Lande & Barrowclough, 1987), or in a severe form could result in inbreeding depression (Charlesworth & Charlesworth, 1987) and contribute to population extinction (Frankham, 2005). The stochastic effects of population founding may also result in changes in the frequency of alleles, including loss of rare variants, increased frequency or fixation of rare variants and decreased frequency of previously common alleles, all of which can result in rapid population genetic divergence (Excoffier et al., 2008). The loss and rearrangement of genetic diversity following population founding also forms the basis of "founder-effect speciation" models (Barton & Charlesworth, 1984;Carson, 1968Carson, , 1975Carson & Templeton, 1984;Matute, 2013;Mayr, 1959;Templeton, 1980Templeton, , 1981Templeton, , 1999. These models, which seek to explain rapid speciation as catalysed primarily by founding stochasticity, remain controversial, partly because it is unclear whether the founder effects upon which these models rely occur frequently in nature (Coyne & Orr, 2004). When a species experiences a new or substantially altered environment, such as following the colonization of an island, the evolutionary trajectory is also influenced by selection in the new adaptive landscape (Price, 2008;Reznick & Ghalambor, 2001), the strength of which is likely to be strongest in the early stages of divergence Ingley & Johnson, 2016;Reznick et al., 1997). However, what is less well understood is how the stochastic effects of population founding interact with selective processes to influence the evolutionary trajectory of a newly founded population.
It is well established that population bottlenecks, especially those that are severe and sustained, can result in loss of genetic diversity (Hunter et al., 2010;Pastor et al., 2004;Roques & Negro, 2005;Weber et al., 2004). However, there are key differences between in situ population bottlenecks versus population founding, including higher potential for population recovery, and continued gene flow from a source population in the latter . As such, loss of diversity is not an inevitable consequence of population founding, and indeed empirical evidence from founded wild populations is mixed. For example, significant loss of diversity has been reported in house finch (Carpodacus mexicanus) (Hawley et al., 2006), Coqui frogs (Peacock et al., 2009); Lupines (Vásquez et al., 2016) and Steller's jay (Cyanocitta stelleri) (Burg et al., 2005), but not in South Island saddleback (Philesturnus carunculatus) (Taylor & Jamieson, 2008), Alpine ibex (Capra ibex) (Biebach & Keller, 2012), white-tailed deer (Odocoileus virginianus) (Fuller et al., 2020), New Zealand sea lion (Phocarctos hookeri) (Foote et al., 2020) or humans (Tabbada et al., 2010) to name a few. Empirical evidence for reduced diversity is also mixed following reintroduction/translocation efforts.
For example, there was no loss of diversity in American martens (Martes americana) (Williams & Scribner, 2010) or Seychelles warbler (Acrocephalus sechellensis) (Wright et al., 2014) but significant reductions in diversity following translocation in Merriam's turkey (Meleagris gallopavo merriami) (Mock et al., 2004) and bridled nailtail wallabies (Onychogalea fraenata) (Sigg, 2006). Several explanations can be proposed for why significant loss of diversity is observed in some studies but not in others. First, the dynamics of a founder event, that is the number of effective founders, the number of generations of small effective population size and the number of founding steps, combine to determine loss of diversity (Nei et al., 1975;Wright, 1931). Where populations are the product of multiple founding steps, for example during a stepwise range expansion, the potential for loss of diversity is expected to be increased and the magnitude of this loss exaggerated . Second, substantial reductions may not be observed in populations that are established from genetically depauperate populations as there is limited diversity to lose Taylor & Jamieson, 2008). Third, conclusions can vary depending on the form of diversity examined and the type of genetic marker used.
For example, allelic richness is more sensitive to diversity loss than heterozygosity (Nei, Maruyama & Chakraborty, 1975) and microsatellites may provide a poorer prediction of genome-wide levels of nucleotide diversity (Roques et al., 2019;Väli et al., 2008) than large numbers of single nucleotide polymorphisms (SNPs) (Kardos et al., 2016).
Empirical systems with sufficient information on both colonization history and the probable selective pressures acting on a population offer a unique opportunity to jointly examine the action of drift and selection at a genome-wide level. One such system, now a classic in the ornithological and evolution literature (Freeman & Herron, 2004;Lack, 1971;Mayr, 1942), is the historically documented, sequential colonization of New Zealand and outlying islands by the Tasmanian subspecies of the silvereye (Zosterops lateralis lateralis) over the last 200 years ( Figure 1a) (Mees, 1969).
These populations have undergone a repeated pattern of phenotypic change in which island populations show a shift towards larger body size and/or more robust bills when compared to their mainland conspecifics ( Figure 1b) (Mees, 1969;. Over this timescale, the degree and rate of these phenotypic shifts can only be explained by invoking directional selection . Recent colonizations by the Tasmanian subspecies are complemented by evolutionarily older silvereye subspecies (Zosterops lateralis chlorocephalus and Zosterops lateralis tephropleurus) that colonized their respective islands (Heron Island and Lord Howe Island) thousands to hundreds of thousands of years ago . These evolutionarily older forms represent some of the largest silvereye subspecies (Figure 1b), and demonstrate that the phenotypic trajectory for an insular silvereye is clearly one of larger body and bill size. Despite the range of island locations, parallel selective pressures could result from exposure to a common suite of biotic conditions often experienced in island populations, namely a combination of reduced predation and a shift in the balance of inter-versus intraspecific competition (Blondel, 2000;Losos & Ricklefs, 2009) that change selective pressures in predictable ways (Grant, 1998). The rare combination of having a suite of recent and evolutionarily older island silvereye populations, with well-documented colonization history, and an understanding of the patterns of selection probably operating within island silvereye populations, provides a powerful system to assess loss of diversity and the accumulation of divergence across the genome following population founding at both neutral and non-neutral loci. Given that genes of large effect have been previously associated with morphological variation in passerines (Bosse et al., 2017;Lamichhaney et al., 2015Lamichhaney et al., , 2016, including within the genus Zosterops (Cornetti et al., 2015;, this system also provides the opportunity to determine whether repeated phenotypic changes following island colonization in the silvereye are produced by a common set of genes or a by suite of alternative genes with similar phenotypic effects ( Barghi et al., 2020;Láruson et al., 2020;Zhu et al., 2018).
Here, using restriction-site associated DNA sequencing (RAD-Seq) (Davey & Blaxter, 2010) of silvereye populations representing F I G U R E 1 (a) The distribution of silvereye forms used in this study. Black circles show sampling locations for the recent colonizer, Zosterops lateralis lateralis, and coloured circles indicate sampling locations for the evolutionarily old forms (Australian mainland; Heron Island; and Lord Howe Island). The timing and direction of the recent sequential colonization by Z. l. lateralis are indicated by dates and arrows. Numbers within black circles indicate the number of founding steps separating the Tasmanian population from derived populations.  2 | MATERIAL S AND ME THOD

| Sample collection
Silvereyes were sampled from nine locations across Australia and islands of the south Pacific ( Figure 1a, Table 1). The populations sam- Unlike other recently founded populations, the Tahitian population is the product of a human-mediated introduction (Thibault & Cibois, 2017). Birds were caught using mist nets or traps and 20-40 µl of blood collected from the brachial wing vein was stored in 1 ml of lysis buffer (Seutin et al., 1991). All sampling was nondestructive, and birds were released at point of capture.

| RAD-sequencing and bioinformatics
DNA extraction and preparation of RAD-Seq libraries was conducted at the University of California, Los Angeles following the protocol outlined by Ali et al., (2016) (Danecek et al., 2011) to remove indels and only include biallelic SNPs where minimum genotype quality = 30, minimum depth = 8, minimum minor allele frequency = 0.01, and where SNPs were called in at least 50% of individuals. Data missingness was then visualized using genoscapertools (https://github.com/eriqa nde/genos capeR tools) ( Figure S1) and the VCF file was further filtered to retain 89,026 SNPs and 145 individuals (those individuals with <30% of data missing). The number of samples retained per population ranged from nine to 21 (Table 1).
As the Z. l. melanops genome is assembled to the scaffold level only (scaffold N50 = 3.5 Mb) (Cornetti et al., 2015), we mapped scaffolds to chromosomes of the Taeniopygia guttata genome assembly version 3.2.4 (NCBI Assembly GCA_000151805.2) (Warren et al., 2010) using satsuma synteny version 3.1.0 (Grabherr et al., 2010). Because chromosomal evolution occurs at an unusually slow rate in birds and synteny between avian genomes is therefore high (Ellegren, 2010), we were able to assign 96.8% of the Zosterops scaffolds to T. guttata chromosomes and determine their order and orientation. Custom

| Linkage disequilibrium
We assessed the extent and patterns of linkage disequilibrium within our RAD-Seq data set by calculating the squared correlation coefficient (r 2 ) between all pairs of SNPs within a physical distance of 1 Mb using poplddecay (Zhang et al. 2019).

| Population structure
We examined patterns of population structure using admixture (Alexander et al., 2009). As the authors recommend avoiding SNPs with high linkage disequilibrium (LD), we used plink (Purcell et al., 2007) to remove one of every pair of SNPs with r 2 > 0 within 100kb sliding windows (LD-pruned data set = 5000 SNPs). We tested K values of 1-9, with 100 independent runs for each value of K, and summarized runs using clumpp (Jakobsson & Rosenberg, 2007).

| Genetic diversity, divergence and the presence of rare alleles
Genetic diversity within populations was measured as average expected heterozygosity (H E ) and average nucleotide diversity (π), both of which were estimated using stacks version 1.4 (Catchen et al., 2013). Significance of pairwise differences in H E and π between populations were assessed using Mann-Whitney tests. Confidence intervals were calculated across loci using bootstrap resampling (1000 iterations). For the recently colonized populations, Spearman's rank correlation was used to test for associations between measures of genetic diversity and the number of colonization steps proceeding population founding.

Genetic divergence between populations was assessed using
Weir and Cockerham's F ST (Weir & Cockerham, 1984). F ST values were calculated, on an SNP-by-SNP basis, between each population using vcftools (Danecek et al., 2011), and 95% confidence intervals were estimated using bootstrap resampling (1000 iterations). As the data are in matrix form, Mantel tests were used to test for associations between the number of colonization steps separating any two populations and mean F ST values, and p-values were estimated using 1000 randomizations.
To assess the effects of population founding on rare alleles, for each recently colonized population, we calculated allele frequencies using vcftools (Danecek et al., 2011). Calculation of allele frequencies was limited to 84,886 SNPs that were present in at least 80% of the individuals within each population. We then identified rare alleles as those SNPs where the minor allele frequency within the original source population of Tasmania was ≤0.1 and counted the number of rare alleles that were lost (allele frequency = 0) and the number of rare alleles that become the major allele (allele frequency ≥ 0.5) in subsequently colonized populations.

| Distribution and accumulation of genomewide divergence
To determine how divergence accumulates across the genome as populations are established via sequential founding steps, we calculated F ST in nonoverlapping 500-kb windows between the Tasmanian population and each subsequently colonized population, using

| Identification of outlier loci and candidate genes
We scanned for outlier loci using pcadapt, a principal componentsbased method of outlier detection with a low rate of false-positive detection (Luu et al., 2017). pcadapt requires the choice of K principal components, based on inspection of a scree plot, where K is the number of PCs with eigenvalues that depart from a straight line. pcadapt then computes a test statistic based on Mahalanobis distance and controls for inflation of test statistics and false discovery rate (FDR). SNPs with high LD were removed using a window size of 200 SNPs and an r 2 threshold of 0.1 and outlier SNPs were then identi- Biomart database was queried using the R package biomart (Durinck et al., 2005(Durinck et al., , 2009. Genes that had previously been associated with body/bill size variation in birds were determined to be candidate genes underlying body size differentiation in silvereyes. Genes known to be involved in craniofacial morphogenesis in nonavian species were also considered candidates. Given the large number of outlier SNPs identified in some comparisons, we performed simulations to assess the probability of candidate genes being identified by chance alone. For each population comparison we selected SNPs at random (without replacement), where the number of randomly drawn SNPs was equal to the number of outlier loci identified using pcadapt, and counted the number of candidate genes that fell within 50 kb of those SNPs. We performed 10,000 simulations per comparison.

| Enrichment analysis
Genes within 50 kb of outlier SNPs were tested for overrepresentation of specific gene ontological terms using the web-based Zebra finch gene ontological (GO) analysis (GOfinch) tool (http://bioin forma tics.iah.ac.uk/tools/ Gofinch). Enrichment was determined using Fisher and hypergeometric tests and a population gene list consisting of 3596 zebra finch genes that were located within 50 kb of all SNPs contained in our RAD-Seq data set. For GO analyses we used human gene ontologies as zebra finch genes are not as well GO-annotated.

| Demographic history
For recently colonized populations, we estimated demographic parameters from the site-frequency spectrum (SFS) using fastsim-coal2 (Excoffier et al., 2013;Excoffier & Foll, 2011). As fastsim-coal2 requires the use of unlinked SNPs and is sensitive to missing data, we filtered the LD-pruned data set to include only SNPs present in all individuals (1551 SNPs). VCFs were converted to SFS format using the R package vcf2sfs (https://github.com/sheng lin-liu/ vcf2sfs) (Liu et al., 2018). As outgroup sequences were unavailable, demographic inference used the distribution of minor allele frequencies. To avoid problems associated with model overparameterization, we tested two simple demographic models. The first (null model: Figure 2a) represents a population of constant size and is defined by a single parameter, the current population size (N cur ).
The second (bottleneck model: Figure

| Linkage disequilibrium
Overall mean LD (measured as r 2 ) among SNP pairs was 0.019 and an average LD of 0.050 was estimated for SNPs less than 1 kb apart.
LD decayed rapidly as the physical distance increased between SNPs ( Figure S2).

| Population structure
Maximum likelihood estimation of individual ancestries using admixture (Alexander et al., 2009) consistently provided the lowest cross-validation error for K = 5 (mean cross-validation error across runs = 0.592). At K = 5, individuals were grouped as follows: (i) Mainland (with some individual exceptions); (ii) Tasmania and sequentially colonized populations excluding Norfolk Island; (iii) Norfolk Island; (iv) Heron Island; and (v) Lord Howe Island (Figure 3).

| Genomic patterns of divergence in recently colonized populations
The distribution of F ST values was most skewed (divergence was localized to few 500-kb windows) for the single-colonization step com-  with divergence between populations separated by three founding steps approaching levels seen between the youngest of the evolutionarily old Zosterops forms (Heron Island) and the mainland population (~1000 generations of separation) (Figure 6a).

| Genetic diversity in recently colonized populations
We found evidence that single founder events often resulted in significant, albeit minor, changes in genome-wide levels of genetic Within the Tasmanian silvereye population 27,950 SNPs were identified as having a rare allele (sites where the minor allele occurred at a frequency ≤ 0.1). Of these, between 8396 and 18,074 were not detected in subsequent colonization steps (Table 2). For the naturally colonized populations (South Island, North Island, Chatham Island and Norfolk Island), the number of undetected rare alleles increased as the number of founding steps from Tasmania increased. The largest loss of rare alleles was detected for Tahiti, a population established via human-mediated introduction. Although the major trend was towards a loss or reduction in the frequency of rare alleles, some rare alleles increased in frequency in subsequent colonization steps (see Figure 6d). In a few cases previously rare alleles became the major allele (allele frequency ≥ 0.5) in the newly founded population (Table 2). Of the 178 rare to major allele transitions, only nine could be attributed to selection, that is transitions occurred at outlier loci identified using pcadapt (Table 2). Values in parentheses indicate the number of rare to common allele transitions that corresponded to the position of outlier SNPs.

| Candidate genes underlying body size differentiation
genome within 50kb of pcadapt outliers identified across comparisons (Table S2). These genes were not significantly enriched for any GO terms (all Bonferroni-corrected p > .05). However, based upon a literature review of gene expression and association studies, 11 genes were previously associated with morphological variation in birds or craniofacial morphogenesis in nonavian species (Figure 5a). These genes were therefore considered candidates underlying the repeated pattern of phenotypic change seen in silvereyes (i.e., larger body/bill size in island populations). The simulations showed that candidate genes were highly unlikely to be identified if SNPs were randomly drawn. Hence the candidate genes were not identified by chance, leading us to conclude that their association with outlier SNPs is meaningful ( Figure 5b). As signatures of selection at candidate genes were inconsistent between population comparisons, this result suggests that a suite of alternative genes may underlie body/bill size variation in the silvereye. However, as we make use of a reducedrepresentation sequencing technique, an alternative explanation is that our SNP data set missed candidate genes that may underlie repeated morphological change in silvereyes. For example, we are unable to speculate on the role of HMGA2 and ALX1 (both thought to F I G U R E 5 (a) Manhattan plot of negative log 10 (p-values) estimated using pcadapt. Points highlighted in pink indicate outlier SNPs identified using FDR = 0.01. The dashed line indicates the threshold above which a SNP is considered an outlier. Genes within 50kb of outlier SNPs and thought to be associated with morphological divergence are labelled. The locations of highly diverged regions (500kb windows with the top 1% of F ST values) are highlighted in blue. Chromosomes are shown with alternating light and dark shading.
Chromosomes are numbered according to the zebra finch nomenclature. Arrows indicate the direction of stepwise colonizations. Numbers within black circles indicate the number of founding steps separating the Tasmanian population from derived populations. (b) Probability that randomly drawn SNPs fall within 50kb of candidate genes identified for a given population comparison based on 10,000 simulations. For each comparison the number of randomly drawn SNPs per simulation was equal to the number of outlier loci identified using pcadapt. Labels indicate the probability of identifying all candidate genes identified for a given comparison by chance alone modulate bill size in Darwin's finches: Lamichhaney et al., 2015Lamichhaney et al., , 2016 as these two genes were not covered by our data set.

| Demographic history
Statistical comparisons of the two demographic models showed strongest support (lowest AIC values) for the bottleneck model

| DISCUSSION
A genome-wide comparison of a historically recent, sequential colonization sequence with more ancient island populations of silvereyes in Australia and the south Pacific has provided both novel insights into the dynamics of population founding at the level of the genome, and valuable empirical validation of some key population genetic principles. We provide empirical evidence that genome-wide divergence accumulates with sequential, natural population founding, and demonstrate that detailed patterns of divergence are highly idiosyncratic in nature, varying across founding steps. Genomic variation is eroded with increasing founder steps, and by tracing the fate of rare alleles across the colonization series, we demonstrate the vagaries of chance, with frequent loss and, much more rarely, increases in frequency of rare alleles. Our results highlight the difference in the potential for founder effects between natural and human-mediated colonizations. In the latter case, a single founder event generated reduced diversity and increased divergence that surpassed the F I G U R E 6 (a) Mean genetic divergence of silvereye populations as measured by pairwise F ST . Error bars indicate 95% confidence intervals calculated across comparisons. The number of founder events is the number of island colonizations separating any two populations. Ages of divergence among subspecies and species were based on molecular and palaeobiological evidence (Degnan & Moritz, 1992;Hopley, 1982). Genetic diversity of silvereye forms as measured by (b) expected heterozygosity and (c) nucleotide diversity. Arrows indicate recent colonization sequence. Error bars indicate 95% bootstrap confidence intervals across SNPs. (d) Changes in the minor allele frequency (MAF) of 100 randomly selected rare alleles for the Tasmania to Norfolk Island colonization sequence. Rare alleles were identified in the Tasmania population (MAF < 0.1) and tracked in subsequent colonization steps. Arrows indicate the direction of stepwise colonizations cumulative effects of the natural colonization series. Eleven genes previously associated with body size variation in vertebrates were also identified as potential candidates underlying larger body size in insular silvereyes, but patterns of selection at these candidates were highly variable, with no single candidate gene under selection across all populations.
Over the past decade a major challenge in speciation genomics has been to understand how different evolutionarily processes shape patterns of divergence at the genomic level (Ravinet et al., 2017;Stankowski et al., 2019;Wolf & Ellegren, 2016). cumulative effects of sequential founder events, supports theoretical (Le Corre & Kremer, 1998) and experimental (Bryant & Meffert, 1993) expectations, along with results from a small number of microsatellite studies of sequential colonization in the wild (Lambert et al., 2005;Michaelides et al., 2018;Pruett & Winker, 2005;Thulin et al., 2011), including a microsatellite study of the silvereye . The greater power and precision provided by the large number of SNPs used here convincingly shows that the erosion of diversity with population founding stems primarily from loss of rare alleles during the colonization process. Following the fate of rare alleles also demonstrated the less common cases of increases in frequencies of rare alleles, as can occur under stochastic or selective processes. Genome-wide F ST increases with successive founder events showed that average F ST was approximately five times higher for populations separated by three founding steps compared to one.
Despite these marked cumulative effects on both diversity and divergence, multiple founding steps were insufficient to produce reductions of diversity and levels of divergence seen in even the youngest of the evolutionary old forms (Heron Island silvereyes).
If a single founder event is more extreme, then stronger effects on genetic diversity and divergence may be observed . This was the case for Norfolk Island  Estoup & Clegg, 2003). In contrast, the founding size of the introduced silvereye population has been described as "a handful of individuals" (Sendell-Price, Ruegg & Clegg, 2020; Thibault & Cibois, 2017), which is supported by our coalescence based demographic inference with only 17 founding individuals estimated, and the isolated nature of the archipelago precludes further gene flow from the source. However, the estimated number of generations to population recovery was low, as is typical for silvereyes. Consistent with previous inferences from microsatellites Estoup & Clegg, 2003), the establishment of the Norfolk Island population appears to have involved fewer founders than other natural colonizations, estimated at 48 individuals using our coalescent-based demographic inference.
While this could contribute to this population's high level of divergence, proposed hybridization with an endemic species, Zosterops tenuirostris (Gill, 1970), may also be a factor. Despite continued and long-running interest in the idea that population founding can be a radical short cut to species generation (founder-effect speciation) (Barton & Charlesworth, 1984;Carson, 1968Carson, , 1975Carson & Templeton, 1984;Matute, 2013;Mayr, 1959;Templeton, 1980Templeton, , 1981Templeton, , 1999, our results from sequential silvereye founding events, in two cases accompanied by relatively narrow demographic bottlenecks, lend empirical weight to the view that founder effect processes are a relatively uncommon mechanism of speciation (Coyne & Orr, 2004), at least within the silvereye system.
We identified 11 candidate genes potentially underlying body/ bill size variation in silvereyes. Candidate genes within 50 kb of outlier SNPs putatively under selection included: BMP7, FOXI1 and TRIM37-all associated with craniofacial defects in vertebrates (Edlund et al., 2014;Graf et al., 2016;Hämäläinen et al., 2004); GSC and NALCN-both of which are associated with bill size variation in Darwin's finches (Lamichhaney et al., 2015;Lawson & Petren, 2017); ID2-which modulates bone morphogenetic protein (BMP) signalling during craniofacial development (Nimmagadda et al., 2015); LIN28Awhich may promote HMGA2 expression (Vignali & Marracci, 2020); NFIA-associated with bill length in the house sparrow (Passer domesticus) (Lundregan et al., 2018) and craniofacial abnormality in humans (Rao et al., 2014); OTX2-which plays a crucial role in craniofacial development across jawed vertebrates (Kimura et al., 1997), and has been associated with bill length in Berthelot's pipit (Anthus berthelotii) (Armstrong et al., 2018); and SOX6-a transcription factor previously associated with bill length in great tits (Parus major) (Bosse et al., 2017). Given that body/bill size follows a highly repeated pattern in the silvereye towards larger size in insular populations Leisler & Winkler, 2015) and that several genes of large effect have been associated with morphological variation in passerines, for example ALX1 and HMGA2 in Darwin's finches (genus: Geospiza) (Lamichhaney et al., 2015(Lamichhaney et al., , 2016, we had expected some continuity in the position of outlier loci across comparisons and for the position of outlier loci to tightly correspond to candidate genes. However, as outlier loci rarely occurred at the same location across comparisons, the stochastic influence of founder events on the genome noted here raises the possibility that if genes of large effect do explain some variation in bill and body size in silvereyes (as opposed to many genes of small effect), then a suite of alternative genes with similar phenotypic effects may be in operation. However, given the marker density of our RAD-Seq data set (which captures ~2.35% of the Z. lateralis genome) we are not yet able to fully describe the genetic basis of morphological change in this species.

| CONCLUSIONS
Our study is the first genome-wide treatment of sequential colonizations in natural populations, and therefore provides an unparalleled empirical demonstration of changes to the landscape of divergence following population founding. Cumulative founding steps generated genome-wide divergence, but this divergence was highly idiosyncratic at the SNP level, emphasizing the role of stochasticity affecting both neutral and putatively selected loci. We did not identify candidate loci of large effect that would explain repeated phenotypic changes observed, though the possibility remains that these will be revealed by whole genome sequencing. Finally, our study also provides the highest precision yet achieved in quantifying not only the pattern but also the magnitude of cumulative decreases in population-level diversity and increases in divergence. The efficacy of sequential founder events in producing these changes was confirmed, though even multiple founder events are not sufficient to produce values seen in evolutionarily older island silvereye subspecies. This furthers the view that founder-event speciation may be rare in nature.

ACKNOWLEDG MENTS
This work was funded by a Natural Environment Research Council and B.C.R.

DATA AVA I L A B I LIT Y S TATE M E NT
Resequencing data produced for this study have been submitted to