Reproductive isolation is a defining feature of speciation and yet hybridisation between species is an 45 important general phenomenon in trait evolution [1, 2] . Among birds, the Anseriformes (ducks, geese and 46 swans) show particularly pervasive hybridisation, 41.6% to > 60% of species hybridising with each other 47 [3, 4] . Domestication generates differentiated gene pools and reproductive isolation between domestics 48 and their wild progenitor, but hybridisation between domestic and wild forms has been well demonstrated 49 in both plants [1, 5] and animals [6, 7] . The impacts include genetic and trait enrichment of domestics, for 50 instance, in chicken the acquisition of a yellow skin phenotype is a result of past mating between red 51 junglefowl and grey junglefowl [8] . In geese, the high tendency for hybridisation among wild species is also 52 shown among wild and domestic forms [9,10], creating an exciting opportunity to study the complex 53 dynamics of hybridisation and domestication. 54
The domestic geese of the world (European and Chinese forms) are derived from two different wild 55 amplification of the library used primers complementary to barcoded and common adapters with products 108 purified as above, and the samples were 100 bp SE-sequenced with Illumina HiSeq 2000/2500 at the BRC. 109 c) GBS pipeline and SNP calling 110 Different models of demographic history were tested with fastsimcoal2 ver 2.6 [41] . We excluded all 141
SNPs that had missing data within the whole data set and executed the analyses with a site frequency 142 spectrum (SFS) based on 9,212 SNPs. As there are no estimates of the genetic diversity per base pair for 143 greylags, we estimated the proportions of variable and monomorphic sites in the data. From the BAM file 144 with -depth option in SAMtools 1.7 [42], we estimated 9,801,382 bp covered with GBS tags. We then 145 mimicked the filtering steps done for the biallelic SNPs to reduce the total number of sites in equivalent 146
proportions. We removed the same number of sites that corresponded to the number of SNPs that were 147 removed because they were indels, had more than 2 alleles or had heterozygosity over 0.75. Since some of 148 the SNPs were removed from this analysis due to missing data in some individuals, we removed an equal 149 proportion of sites from the total number of sites as well. The final SFS had 1,681,316 sites. 150
To infer the demographic history, we chose a subset of individuals from both wild-collected greylags 151 and domestic geese to represent the genetic variation in both groups. Therefore, 11 greylags with > 90.8% 152 of greylag ancestry and 15 domestic geese with > 91.4% of European domestic goose ancestry were 153 selected for the analysis. The mutation rate for the simulations was 1.3810 -7 [43]. The parameter 154 estimation for each model tested involved 100,000 simulations and 40 conditional maximization (ECM) 155 cycles. The parameters for each model were estimated with 100 independent runs to obtain the global 156 maximum. The models tested were i) simple divergence of two populations with no migration, ii) 157 divergence of two populations with migration and iii) divergence of two populations with two different 158 migration matrices (Figure 2 ). The best model was selected based on Akaike's weight of evidence as in [41] . 159
For parametric bootstrapping 100 SFS were simulated with the parameter estimates obtained from the real 160 SFS, followed by maximum likelihood estimation with 50 independent runs for each bootstrap SFS. The 95% 161 confidence intervals were obtained from the bootstrap data for each estimated parameter. 162
f) The estimation of genetic diversity 163 (Figure 3a-b) . At K = 2, populations/breeds are clustered based on their status (wild or 177 domestic) and, at K = 3, domestic geese are further separated into European and Chinese. At K = 4, the 178 fourth cluster is admixed within greylag populations but none of the individuals are unanimously assigned 179 to that cluster. The likelihood was highest for K = 3. These results were supported by PCA as the first two 180
PCs out of 14 significant PCs (p < 0.05) were enough to separate the three groups (wild, European domestic, 181
Chinese domestic) from each other (Figure 3a ). Overall the greylag populations showed 3.2% -23.5% 182 admixture proportions with European domestic geese when K = 3. In contrast, not all European domestic 183 geese showed admixture with greylags and the admixture percentages ranged from 0.0 to 8.4%. At K = 3 184 many European domestic goose breeds showed mixed ancestry with Chinese domestic geese (0.0-27.1%). 185
The neighbor-joining tree repeated the major patterns observed with STRUCTURE and PCA revealing 186 a star shaped phylogeny and confirming that the domestic and greylag geese largely form different clades 187 (ESM, Figure S1) . Surprisingly, the Chinese domestic geese were closer to European domestic geese and 188 greylags, than to the swan goose reference genome. In addition, one greylag from Turkey was more closely 189 Although, STRUCTURE implied considerable mixed ancestry from multiple genetic clusters for Dutch and 199
Turkish greylags, the ƒ3 analysis did not confirm admixture for these populations. The two Turkish greylag 200 samples came from the same area as our NW-Turkish domestic population, which among Turkish domestic 201 geese showed highest admixture with greylags (2.2%), but admixture was not confirmed with the ƒ3 test. 202
The ƒ3 analysis confirmed admixture of domestic geese in line with STRUCTURE results. Most notably, 203
the African breed is a hybrid between European and Chinese domestic geese (Z-score -6.399), unexpected 204 as this breed has been assumed to have originated solely from swan goose. The European-Chinese hybrid 8 respectively). The Kholmogor breed also fell halfway between European and domestic geese both in 207 STRUCTURE and PCA, whereas the Steinbacher was genetically closer to European domestic geese in the 208 PCA. However, the Diepholzer breed, which reportedly is also a hybrid, was not confirmed as such in our 209 analysis. Other domestic breeds/groups with admixture status in STRUCTURE were also confirmed to have 210 a European-Chinese admixture when a Z-score threshold of -3 (roughly corresponding to p < 0.01) was 211 used: Sebastopol, Toulouse cross, Domestic NY, Embden, Tufted Roman (Table S2 , Figure 3c ). The Crested 212
Faroese breed gave indication of admixture based on STRUCTURE analysis and the ƒ3 test supported this (Z-213 score of -2.228, p < 0.05). Surprisingly, the Northern Turkish domestic population was not admixed with 214 Chinese domestic geese in STRUCTURE but ƒ3 analysis gave a contrasting signal (Z-score -2.459, p < 0.05). 215
The demographic model that best fit our data suggested divergence of greylag and domestic geese 216 with a recent migration rate change ( 
We studied the dynamics of domestication and hybridisation in grey (Anser) geese using genome-wide SNP 246 data. The results demonstrated genetic divergence between Eurasian wild greylag and European domestic 247 geese with long-term genetic exchange between them. We also inferred temporal changes in the direction 248 of gene flow. The degree of hybridisation between greylag and domestic geese also varied geographically. 249
Surprisingly, several domestic goose breeds also showed a substantial genetic contribution of Chinese 250 domestic geese. We also provide insights about the origin and the timing of goose domestication. domestic geese was not sampled in this study and must, therefore, be sought by ancient DNA approaches. 283
It is, however, tempting to make our own determination of the location of goose domestication from 284 the modern data presented here. Given that genetic diversity would be expected to be highest in the 285 'domestication centre' and reduce with increasing distance from there, the high mtDNA diversity of Turkish 286 domestic geese means the eastern Mediterranean remains a strong candidate for the origin of goose 287 domestication, and our nuclear data accord with that. However, as we have shown, hybridisation between 288 wild and domestic geese can also generate high genetic diversity both within and outside the original 289 domestication location. Table S1 . The colours in a) and b) are associated to different groups as 570 follows: greylags (blue), European domestics (green) and Chinese domestics (red). a) The first three 571 principal components summarising the genetic variation in geese (percentage explained by each PC is 572 shown). b) STRUCTURE assignment plots for K=2, K=3 and K=4. Each vertical line represents one individual 573 with K number of colours indicating proportion of ancestry from the inferred clusters. c) Plot relating to the 574 ƒ3 (Table S2) values obtained for each population. The more negative the ƒ3, the more significant is Z-score 575 in favour of admixture. The ƒ3 values were not calculated for Landes and the Chinese geese, as they were 576 used as source populations, thus they were given an ƒ3 value of 0.