Genome-culture coevolution promotes rapid divergence of killer whale ecotypes

Analysing population genomic data from killer whale ecotypes, which we estimate have globally radiated within less than 250,000 years, we show that genetic structuring including the segregation of potentially functional alleles is associated with socially inherited ecological niche. Reconstruction of ancestral demographic history revealed bottlenecks during founder events, likely promoting ecological divergence and genetic drift resulting in a wide range of genome-wide differentiation between pairs of allopatric and sympatric ecotypes. Functional enrichment analyses provided evidence for regional genomic divergence associated with habitat, dietary preferences and post-zygotic reproductive isolation. Our findings are consistent with expansion of small founder groups into novel niches by an initial plastic behavioural response, perpetuated by social learning imposing an altered natural selection regime. The study constitutes an important step towards an understanding of the complex interaction between demographic history, culture, ecological adaptation and evolution at the genomic level.

Correlation of per-site F ST between low-coverage whole genome sequencing (WGS) data generated for this study and high coverage (>20×) published RAD-seq data. a, Per-site F ST estimates from a pairwise comparison of WGS data of 10 residents and 10 transients plotted against F ST estimates of the same 547 polymorphic sites from a pairwise comparison from RAD data of 52 residents and 37 transients 19 . b, Distribution of the correlation coefficients (r) of the per-site F ST estimates from a pairwise comparison of WGS data of 10 residents and 10 transients with 1,000 random re-samplings with replacement of 10 residents and 10 transients from the RAD-seq dataset 19 . were further mapped to a reference mitochondrial genome (GU187176.1) and compared with previously published mitogenome sequences from these individuals 11 . The assembled mitogenome sequences were a 100% match with those previously generated for these individuals using targeted sequencing approaches 11 . As previously reported, the mitogenomes of each ecotype clustered in strongly supported mitochondrial DNA clades, with the exception that one type B1 individual sampled at a different geographic location to the other type B1 individuals, had a highly divergent mitogenome haplotype 11 . transient resident type B1 type B2 type C b Supplementary Figure 3b, Distance-based tree generated from the 48 low coverage and one high (~20×) coverage nuclear genome sequences. We generated 100 matrices of pairwise genetic distances, in which pairwise genetics distances were calculated using ngsDist 72 , which takes genotype uncertainty into account by avoiding genotype calling and instead using genotype posterior probabilities estimated by ANGSD. A block-bootstrapping procedure was used to generate 100 distance matrices, obtained by repetitively sampling blocks of the original data set (Supplementary Data 3). Pairwise genetic distances were visualised as a phylogenetic tree using the distancebased phylogeny inference program FastME 2.0 73 . Individuals largely clustered by ecotype indicating that segregating alleles are shared among individuals within each ecotype. The short branches that did not cluster as closely to ecotype were the individuals with the least sites covered at ≥2× (see Supplementary Table 1). We estimated pairwise relatedness due to identity-bydescent (IBD), i.e. genetic identity due to a recent common ancestor, of every possible combination of two individuals using NgsRelate 74 . NgsRelate provides ML estimates of R, where R = (k0,k1,k2) and km is the fraction of genome in which the two individuals share m alleles IBD. NgsRelate provides maximum likelihood estimates of R by finding the value of R that maximizes this likelihood function with an Expectation Maximization algorithm using genotype likelihoods instead of genotypes to account for the inherent uncertainty of the genotypes. NgsRelate has been shown using simulations and real data to provide robust estimates for low-depth NGS data (as low as 1-2×), which are markedly better than genotype-based methods 74 Figure 5. a, Maximum-likelihood graph from TreeMix. The scale bar shows ten times the average standard error of the entries in the sample covariance matrix. b, Residual fit of the observed versus predicted squared allele frequency difference, expressed as the number of standard errors of the deviation. Colours are in the palette on the right. Residuals above zero represent populations that are more closely related to each other in the data than in the best-fit tree, and are candidates for admixture. c, Maximum-likelihood graph allowing a migration event to improve the fit of the tree to the data. residents and the Atlantic sample, with the bottlenose dolphin as the outgroup. This statistic identifies an excess of shared derived alleles between taxa, which could result from introgression or ancestral population structure. The statistic can thus be used to identify departures from 'treeness' of a given topology. For example, if H1, H2 and H3 are taken to denote 3 ecotypes, the test can be used to evaluate if the data are inconsistent with the null hypothesis that the tree (((H1, H2), H3), dolphin) is correct and that there has been no gene flow between H3 and either H1 or H2 or any populations related to them. The definition used here is from reference 75: where nABBA is the number of sites where H1 shares the ancestral allele with the dolphin, and H2 and H3 share a derived allele (ABBA sites); and, nBABA is the number of sites where H2 shares the ancestral allele with the dolphin, and H1 and H3 share a derived allele (BABA sites). Under the null hypothesis that the given topology is the true topology, we expect an equal proportion of ABBA and BABA sites and thus D = 0. Hence a test statistic that differs significantly from 0 provides evidence either of gene flow or the tree being incorrect due to ancestral population structuring. The significance of the deviation from 0 was assessed using a Z-score based on jackknife estimates of the standard deviation of the D-statistics. This Z-score is based on the assumption that the D-statistic (under the null hypothesis) is normally distributed with mean 0 and a standard deviation equal to a standard deviation estimate achieved using the "delete-m Jackknife for unequal m" procedure. The tests were implemented in ANGSD and performed by sampling a single base at each position of the genome to remove bias caused by differences in sequencing depth.
The positive values over the critical value of 3 indicate an excess of ABBA patterns over BABA patterns in terms of the number of shared derived alleles. This indicates that the relationship among these taxa is not fully described by a bifurcating tree model, but that instead ancient admixture occurred between the ancestral populations of the North Pacific transient samples and the North Atlantic samples included here.  Overlaying the two plots, (c) we observe that the two genomes show no convergence in effective population size (N e ) prior to the date that they are estimated by the same authors to have shared a common ancestor 22 or even back as far as 1 MYA. This highlights the bias introduced into the analysis when comparing sequences of low and differing coverage, which results in different rates of false negative detection of heterozygote sites, producing the same effect as using a lower mutation rate for the sequence with lower coverage.
A recent study used PSMC to reconstruct ancestral changes in Ne through time using ~13× and ~20× coverage diploid autosomal genome sequences of a North Atlantic a b c killer whale (accessed ahead of publication by the consortium that generated the data) and a North Pacific killer whale respectively 17 . The authors interpreted the plots as being indicative of a global decline driven primarily by climate change during the last glacial period of the Pleistocene 17 . This study dismissed changes in connectivity having a role in the observed demographic changes as being 'unlikely to generate the specific pattern observed (strong population decline) or the very similar profiles for each ocean 17 . However, PSMC heavily relies on the distribution of polymorphic sites across the genome, and in particular, the length of shared runs of homozygosity, and can be biased when heterozygous sites are wrongly called as being homozygous. PSMC plots of genomes with <20× coverage have been shown not to be directly comparable to higher coverage genomes without first applying a correction for an appropriate false negative rate of detecting heterozygotes 76 . Additionally, the effect of mapping short read data to a reference genome comprised of short contigs may also be problematic, especially when many contigs fall below ~50-kb (the typical size of shared fragment size from 1,000 generations ago in humans), although to the best of our knowledge, the effect of mapping to different quality reference assemblies on PSMC analysis has not been tested to date. Moura et al. 17 mapped short read data from two individual killer whales to a draft assembly of the bottlenose dolphin (assembly turTru1, Ensembl database release 69.1) made up of 0.24 million scaffolds, with a scaffold N50 of 116,287, in which 94% of the scaffolds, comprising approximately 25% of the genome, are less than 50-kb long. Moura et al. 17 thereby generated a 20× average coverage sequence and a 13× average coverage sequence, and did not apply a correction for differences in false negative rate of detection of heterozygotes due to the difference in coverage. The PSMC plots from the two genomes presented in Moura et al. 17 do not converge in effective population size even though the two individuals shared a relatively recent common ancestor (TMRCA estimated at ~150 KYA by the same authors 22 ). Figure 10. PSMC plots of inferred demographic history of the same individual using 13×, 20× and 50× coverage of sequencing data. In order to better understand if these methodological issues (<20× coverage sequence data mapped to a highly fragmented reference) led to erroneous inference of the demographic histories and the underlying processes in this previous study 17 , PSMC was used to analyse down-sampled versions of the high coverage North Atlantic killer whale genome. A 50× bam file was produced by mapping the short read data generated from the North Atlantic killer whale to scaffolds of the autosomal regions of the high quality killer whale reference genome that were greater than 10-Mb in length, totaling 1.5 Gb and which had all repeat regions masked as noted above. The 50× coverage bam file was then down-sampled to produce 13× and 20× coverage bam files. A consensus sequence of each of the three bam files was then generated in fastq format sequentially using: firstly, SAMtools mpileup command with the -C50 option to reduce the effect of reads with excessive mismatches; secondly, bcftools view -c to call variants; lastly, vcfutils.pl vcf2fq to convert the vcf file of called variants to fastq format with further filtering to remove sites with less than a third or more than double the average depth of coverage and Phred quality scores less than 30. The PSMC inference was then carried out using the recommended input parameters for human autosomal data 21 , i.e. 25 iterations, with maximum TMRCA (Tmax) = 15, number of atomic time intervals (n) = 64 (following the pattern (1*4 + 25*2 + 1*4 + 1*6), and

50X 20X 13X
initial theta ratio (r) =5. For the initial comparison between the 13×, 20× and 50× coverage North Atlantic genomes, a generation time of 25.7 years and a mutation rate of 1.53×10 -8 substitutions/nucleotide/generation were applied, as per reference 17. Comparison of the PSMC inference plots based on the 13×, 20× and 50× coverage files, generated from the same individual, highlighted the impact of coverage on inference of both the magnitude of Ne at any given time and the timing of the changes in Ne, consistent with findings by a previous study 76 . In particular, estimates of Ne in more recent times based on the 13× genome assembly differed markedly to inferred Ne from the 20× and 50× genome assemblies. This is a consequence of a higher false negative detection rate of heterozygote sites in the 13× genome assembly, producing the same effect as a smaller mutation rate would have on the plot. The PSMC plots of the 20× and 50× coverage North Atlantic genome were almost identical both regarding the timing and the magnitude of demographic events. Figure 11. Historical population sizes of a North Atlantic (red) and North Pacific resident killer whale (blue) inferred by pairwise sequential Markovian coalescent analysis (PSMC). All three plots of the North Atlantic killer whale genome (13×, 20× & 50×) in Supplementary Figure 10 are consistent in estimating a marked decline in Ne between 100,000 years and 20,000 years ago. To better infer the process underlying this decline in Ne, PSMC was used to compare equal coverage (20×) assemblies of the North Pacific and North Atlantic genomes. A 20× bam file was generated for the North Pacific killer whale using data from the short read archive (SRP035610) 17 and mapped as above. As with other inference methods based on coalescent theory, PSMC can only infer scaled times and population sizes. To convert these estimates into real time and size, all scaled results need to be divided by the mutation rate. To allow comparison of the relative timing of population splits with a published time-calibrated nuclear phylogeny based on RAD-seq data we scaled the PSMC plot using the same mutation rate as reference 22. Although the two papers by Moura et al. 17,22 were published almost concurrently, they use two different mutation rates for nuclear genomic data for each analysis: 1.53×10 -8 substitutions/nucleotide /generation for PSMC 17 and an estimate almost double this rate for their timecalibrated phylogeny of 2.83×10 -8 substitutions/nucleotide/generation, based on their given rate of 0.0011 substitutions per site per million years 22 and a generation time of 25.7 years as above. We therefore scaled the PSMC plots by a generation time of 25.7 years and a mutation rate of 2.83×10 -8 substitutions/nucleotide/generation. A total number of 100 bootstraps were performed. The combined PSMC plot of both genomes is shown compared with population split times from a previously published time-calibrated phylogeny 22 that supports the inference that changes in inferred Ne by PSMC are at least partially driven by changes in connectivity. The x-axis gives time measured by pairwise sequence divergence and the y-axis gives the effective population size measured by the scaled mutation rate of 2.83×10 -8 substitutions /nucleotide/generation and assuming a generation time of 25.7 years. Thin light lines of the same colour correspond to the 95% confidence intervals of PSMC inferences on 100 rounds of bootstrapped sequences. Insert shows a time-calibrated nuclear marker phylogeny adapted from reference 22, scaled using the same mutation rate and plotted on the same x-axis as the PSMC plots. The branches leading to the populations used in the PSMC plots are coloured accordingly and highlight that population splits are followed by changes in inferred N e .

Supplementary
Supplementary Figure 12. Effect of population structure and divergence on PSMC. Finally, to further understand the relationship between connectivity and PSMC inference we simulated data consistent with a population split using scrm 77 and plotted the changes in effective population size inferred by PSMC. We considered two models of population divergence where an ancestral population of size N A =20,000 splits 6,000 generations ago into two descending populations of size N 1 =N 2 =10,000. In a) populations become isolated after the split event and diverge without gene flow, representing a scenario of a sudden change in connectivity; whereas in b) there is a period of symmetric gene flow (2Nm=20) until 1600 generations ago, mimicking a gradual change in the connectivity of populations. We considered models with no changes in the size of populations, which in total remains 20,000, and hence these models could represent scenarios of divergence due to vicariance without founder events. As can be seen in the plots on the right, reductions in the connectivity of populations due to population divergence can lead to changes in the effective sizes inferred with PSMC. For both models, even though we did not simulate any demographic changes due to population bottlenecks or expansions, PSMC infers a gradual decline on the effective sizes after populations become more isolated, coinciding with the reduction in connectivity among populations. Note that these simulations are not intended to capture the recent demographic history of the killer whales, but rather illustrate that changes in population structure due to population divergence can lead to changes in the inferred PSMC effective sizes. We simulated and analysed with PSMC eight independent datasets for each model (corresponding to the different lines in the PSMC plots), assuming a mutation rate m of 2.34x10 -8 per site per generation, an arbitrary recombination rate set to 0.80*m , and a generation time of 25.7 years. Each datasets consisted of 10 blocks of 100 Mb (total of 1,000 Mb) sampled from a single diploid individual from population 1, generated using scrm with the following command lines: a) ms 2 10 -t 93600 -r 74880 100000000 -I 2 2 0 0.0 -n 1 1 -n 2 1 -G 0.0 -ej 0.15 2 1 -en 0.15 1 2; b) ms 2 10 -t 93600 -r 74880 100000000 -I 2 2 0 0.0 -n 1 1 -n 2 1 -G 0.0 -m 1 2 0 -m 2 1 0 -em Supplementary Figure 13. Population-specific allele frequency changes estimated as F ST -based branch lengths in 50-kb sliding windows. The mean population branch statistic (PBS) is highest along the two branches that were inferred by TreeMix to have undergone the highest amounts of drift: the branch to the common ancestor of the Antarctic types and the branch to the resident ecotype. The green and red lines indicate the 99.5 and 99.9 percentiles respectively. Several distinct peaks are seen in the Manhattan plots for each branch, these are further explored using PBS estimates at the exon level to identify candidate loci that have potentially evolved under selection.

Supplementary Figure 14.
Correlations between (50-kb) window-based estimates of genome-wide divergence (Dxy), differentiation (F ST ) and nucleotide diversity (π) for the pairwise comparison between ecotypes shown in Fig. 4A. Red circles indicate a positive relationship, blue a negative one; colour intensity and circle size are proportional to Spearman's correlation coefficient. Shared regions of high differentiation, but low diversity and divergence are likely to have been regions on linked selection on the ancestral form, which would remove diversity and result in rapid lineage sorting of allele frequencies due to differential drift and selection in the derived forms. The individual was travelling alone and moving slowly. The skin was peeling off as evidenced in the photograph. Those of us that have conducted field studies over many years have not previously encountered any killer whales with a skin condition such as this. A previous study 42 has detailed that Antarctic killer whales make rapid round-trips of up to 9,400 km from Antarctica (less than 60° South) to subtropical waters (30-37° South) where sea surface temperature was approximately 20-25°C warmer using satellite tag data. Durban & Pitman 42 additionally documented that the same individuals were encountered with the different amounts of accumulation of diatoms on their skin at different times. They hypothesised that these rapid migrations to warm water may allow for epidermal tissue regeneration without the thermal cost that would be incurred if skin regeneration took place in Antarctic waters 42 . The type B1 individual photographed by Conor Ryan in Antarctic waters highlights that skin regeneration may be a strong selective force. The results of this study highlight genomic adaptation of a gene (FAM83H) associated with epidermal regeneration along the branch to the ancestral Antarctic lineage that include four non-synonymous substitutions. To make the autosomal plot and X-chromosome synchronise in the inferred timing of demographic change, requires scaling the X-chromosome plot by a mutation rate of 1.0×10 -8 substitutions/nucleotide/generation if the autosomal mutation rate is assumed to be 2.34×10 -8 substitutions/nucleotide/generation 20 . This requires a male-to-female mutation rate ratio to be >100, making it seemingly biologically unrealistic.

A brief natural history of the study species
The killer whale is emerging as a useful organism for studying adaptation and speciation, as the phenotype, biogeography and ecology underlying evolutionary divergence and the genetic outcome in terms of neutral genetic differentiation are well described 13,78 . Dietary differences have been studied through: direct observation of naturally marked, site-faithful individuals over many years; multi-chemical markers such as stable isotope and fatty acids; and molecular and visual identification of prey remains from predation events, faecal samples and stomach contents 7-10,79-89 .
Morphology has been described qualitatively and quantitatively. For example, body length has been measured directly from stranded and captive specimens or those taken by whaling operations, or from free-ranging live specimens using laser-metrics and Long-term studies of naturally marked individuals have detected no dispersal between North Pacific ecotypes and no dispersal from the natal matrilineal social in the resident ecotype 99,104 . There is some dispersal of transients from the natal group 102,105 .
The lack of dispersal between ecotypes appears to have been maintained over longer timescales than the field studies based on lineage sorting of mitochondrial genomes and significant differentiation between ecotypes based on microsatellite allele frequencies 11, [106][107][108][109] . Better resolution on the timing and extent of any gene flow between North Pacific ecotypes is needed. Phylogeographic analyses based on mitogenome sequences indicate that the resident and offshore ecotypes share a more recent common ancestor with lineages of Northeast Atlantic killer whales than with the transient ecotype and are consistent with sympatry between the North Pacific ecotypes arising from secondary contact following an allopatric phase 110 . However, phylogenies based on nuclear loci are to some extent discordant with the mitochondrial phylogeny 22 , possibly due to low levels of gene flow between ecotypes within the same ocean basin, which may have occurred upon primary or secondary contact, making it difficult to discern between these two scenarios 111 . There are several genetically differentiated populations of the resident and transient ecotype in the North Pacific 106-108 , but to date only one population of the offshore ecotype has been identified 106,107 .
Killer whales in the waters around the Antarctic continent have diversified into several distinct morphotypes partially overlapping in their ranges 9,10,81,82 . Killer whales in Antarctic waters with the pigmentation patterns that most closely resemble the common killer whale colouration are morphologically classified as type A 9 , but this classification does not infer genetic or ecologically cohesiveness. The Antarctic morphotypes included in this study (types B1, B2 & C) differ from type A as they have a discernable dorsal cape 9,10 . Types B1 and B2 have a large eye patch 9,10 , whereas type C has a smaller forward-slanted eye patch 9 . Body size also varies, with photogrammetry measurements indicating that type B2 is smaller than type B1 10 , and that type C grow up to just 5.6 meters in length, making this the smallest form of killer whale measured to date 94 .
Field observations and stable isotope measurements indicate that there are differences in the preferred habitat and prey of each of the Antarctic types 9,10,81,82,85 . Type B1 is commonly observed in the pack-ice hunting Weddell seals (Leptonychotes weddellii) 10,82 , whilst type B2 forages in more open water, observed killing and eating penguins 10,83 ; type C is most commonly observed in the dense pack-ice and its diet is known from observations to include Antarctic toothfish (Dissostichus mawsoni) and based on stomach contents from Soviet whaling data is thought to be primarily pisciverous 9 . Most observations of these Antarctic types have been made during the Austral summer, however, there are some observations of type B and type C in the Antarctic pack-ice during the Austral winter 9 . There have also been occasional sightings of these Antarctic types at higher latitudes 9 and satellite-tagging data indicates that they make rapid round-trip movements from Antarctic waters to subtropical waters and back, hypothesized to be for skin generation in warmer waters 42 .
Genetic differentiation based on microsatellite allele frequencies is relatively low between type B and type C (G´S T = 0.11) compared with differentiation between the North Pacific resident and transient ecotypes (G´S T = 0.28) 112 . Mitogenome phylogenetic analyses indicate that type B (both B1 and B2) and type C are reciprocally monophyletic, suggesting that there is little or no permanent dispersal between types 11,112 , with the exception of a single sampled type B1individual 11 .
Comparison of amino acid substitutions across the mitogenome identified two nonsynonymous changes resulting in localized changes in polarity that putatively occurred under natural selection within the cytochrome b gene 113 . One change had reached fixation in type B killer whales (both B1 and B2) and the other was close to fixation in type C 113 . The changes were at different sites and in the opposite direction in each type, suggesting divergent evolution since types B1, B2 and C diverged from their most recent common ancestor. All other substitutions across the mitogenome in a global killer whale dataset appeared to have evolved under neutrality 113 .
The published data cited here indicate that these populations of killer whales are ecologically, morphologically and genetically divergent, and thus generally meet the criteria for most of the many definitions of the term 'ecotype' that have been proposed over the years 114 . The term ecotype is therefore adopted here. The behavioural adaptations that each ecotype uses to exploit an ecological niche are thought to be passed on from one generation to the next by social learning within matrilineal groups 12,13 . These behavioural adaptations include: coordinated 'wavewashing' behaviour by type B1 killer whales in Antarctica to dislodge seals from ice floes 82 ; 'carousel-feeding', whereby killer whales in some North Atlantic populations are reported to co-ordinately herd herring schools into a tight ball by encircling them, flashing their white undersides, emitting large bubbles and producing a low frequency pulsed call, prior to tail-slapping the herded herring to stun them 115,116 ; and intentional stranding on to the beach to catch seals performed by some killer whale social groups at the Crozet Archipelago in the Southern Ocean 117,118 . Perhaps due to the complexity and cumulative nature of human culture and due to it being the focus of study in a range of fields from anthropology to zoology, there is no clear definitional consensus