Niche variation owing to individual differences in ecology has been hypothesized to be an early stage of sympatric speciation. Yet to date, no study has tracked niche width over more than a few generations. In this study, we show the presence of isotopic niche variation over millennial timescales and investigate the evolutionary outcomes. Isotopic ratios were measured from tissue samples of sympatric killer whale Orcinus orca lineages from the North Sea, spanning over 10 000 years. Isotopic ratios spanned a range similar to the difference in isotopic values of two known prey items, herring Clupea harengus and harbour seal Phoca vitulina. Two proxies of the stage of speciation, lineage sorting of mitogenomes and genotypic clustering, were both weak to intermediate indicating that speciation has made little progress. Thus, our study confirms that even with the necessary ecological conditions, i.e. among-individual variation in ecology, it is difficult for sympatric speciation to progress in the face of gene flow. In contrast to some theoretical models, our empirical results suggest that sympatric speciation driven by among-individual differences in ecological niche is a slow process and may not reach completion. We argue that sympatric speciation is constrained in this system owing to the plastic nature of the behavioural traits under selection when hunting either mammals or fish.
Ecological variation is the raw material by which natural selection can drive evolutionary divergence [1–4]. This variation can take the form of differences between specialist populations, but also individual differences in niche within a generalist population [5–8]. Recent studies have suggested that strong natural selection owing to niche variation within a generalist population (e.g. among-individual differences) can promote adaptive divergence of phenotypic traits and assortative mating [9–12], and niche variation is therefore a potential early stage of sympatric speciation [13–15]. Niche variation would need to be present and stable in a population over evolutionary timescales for disruptive selection to promote directional progress towards speciation [8,15]. Most evolutionary studies typically use either comparisons at a single point in time or over timescales representing one to a few generations [15–17] or infer ancestral states using phylogenetic-based methods . This snapshot of evolution may only be partially informative, as throughout the speciation process, diverging populations may have experienced changes in ecological conditions, geographical distribution and population size that could have influenced the strength of selection and the rate of progress towards speciation [19–22]. However, real-time tracking of niche width in long-lived species over evolutionary and ecological timescales can be achieved by the incorporation of ancient DNA (aDNA) and stable isotope data from subfossil specimens . Here, we use stable isotope and aDNA analyses to track the niche width and evolutionary history of apparent generalist lineages of killer whales Orcinus orca.
2. Study system
Killer whales are marine top predators and as a species, are known to consume a broad range of prey including species of mammal, fish, bird and reptile . However, several studies have identified highly specialized populations that consume a narrow range of prey [25,26]. For example, in the near shore waters of the Northeast Pacific, there are two sympatric ecotypes of killer whales, which are thought to have non-overlapping prey preferences: one ecotype feeds on fish, whereas the other feeds upon mammals [25,26]. These North Pacific ecotypes are reproductively isolated, as indicated by strong lineage sorting of mitochondrial DNA (mtDNA) and strongly bimodal genotypic clustering. They therefore appear to be at a late stage, or have even reached completion, of speciation . Phylogeographic and coalescent analyses using complete mitogenomes suggest that sympatry between these two North Pacific ecotypes arose from secondary contact following an allopatric phase . Given this, one can speculate whether foraging specialization and evolutionary divergence were initiated during allopatry or upon secondary contact.
There is evidence that different killer whale populations from elsewhere in their natural range, e.g. the Northeast Atlantic, provide a suitable model for testing the potential that evolutionary divergence can be initiated and progress in sympatry, and in doing so, provide broader insights into the evolutionary outcomes of niche variation at a more general level. In the Northeast Atlantic waters, some degree of ecological diversification between groups, and relative seasonal specialization within groups, have been reported [29–31]. However, in contrast to the findings from studies in the Pacific, individuals sharing the same mtDNA haplotype were found to have very different δ15N stable isotope values indicating differences in relative trophic position . This niche variation within a lineage could be the raw material needed for adaptive radiation, and therefore the early stages of ecotype formation and ultimately speciation . To fully investigate whether such niche variation has been present over evolutionary timescales, we compared isotopic values of killer whale specimens, covering a temporal span from the early Holocene to the present. The samples were collected from a relatively small geographical area (≈1000 km2) of the Northeast Atlantic, predominantly the North Sea, which is a small proportion of the natural range of killer whales (North Sea killer whales have been photographically recaptured off Iceland , approximately 1000 km away). In addition to using isotopic values as an indicator of the degree of ecological divergence, we further investigated the evolutionary outcomes of any differences in ecology using two proxies of the stage of speciation: the degree of lineage sorting of mitochondrial genomes and genotypic clustering of microsatellite loci.
3. Material and methods
We sampled 23 subfossil killer whale bones and teeth recovered by dredging or trawling from the Southern Bight of the North Sea or from archaeological sites in Southern Scandinavia (figure 1). Four of the subfossils were radiocarbon dated to 2810 ± 75; 3250 ± 95; 3900 ± 40 and 6800 ± 115 14C years BP. Others were dated based on archaeological context, for example presence within Iron Age middens. All were estimated to originate from the Holocene and to be greater than 1000 years old. Species identification was confirmed by PCR-based Sanger sequencing of the hypervariable 5′ region of the mtDNA control region of all samples as per reference . These were compared with 20 bones or teeth from more recent museum specimens (1865–1995), which originated from the North Sea and adjacent waters (figure 1) and were used in a previous study .
(a) Estimating isotopic niche width
Stable isotope ratios from tooth or bone can provide a relative measure of the long-term average consumption of different prey resources and variance in these isotopic ratios among individuals can be a useful proxy of the population's niche width [33,34]. We only included bone collagen and tooth dentin from animals estimated, based on size to be older than the post-weaning age of 4 years  and avoided sampling the enamel crown and outer dentin growth layers deposited during nursing. Nitrogen and carbon isotope analyses were performed simultaneously using continuous-flow isotope ratio mass spectrometry under the same conditions and degree of accuracy as reference . To compare isotope values with more recent samples, we applied a correction to δ13C estimates to account for the ‘Suess effect’ caused by the burning of fossil fuels over the past 150 years. Time-dependent corrections of −0.005‰ per year and −0.022‰ per year were applied to periods 1860–1960 and 1960–present, respectively. Isotopic niche width was estimated using a Bayesian approach based on multivariate, ellipse-based metrics , as this method is robust for comparisons between small and different sample sizes and identifies differences in the niche width of ‘typical’ members of the population and may not incorporate outlier individuals. The analysis was implemented in the R package SIAR  to generate standard ellipse areas (SEAB): a bivariate equivalent to standard deviation and a corrected measure for small sample sizes (SEAC). The area within an ellipse is defined by a subsample (40%) of bivariate data, in this case, the ratios of nitrogen 15N/14N (δ15N) and carbon 13C/12C (δ13C), that best explain the covariance, and by resampling multiple times estimates an error term associated with this value . Statistical significance of differences in SEAC between sample sets was based on the proportional outcome of 106 repeat measures.
The differences in amino acid composition among different tissues can lead to large differences in trophic discrimination . For example, δ13C in bone collagen is typically enriched by 4–5‰ relative to the diet, in comparison to a 0.9–1.9‰ enrichment typical for epidermal tissue . We therefore applied a correction factor of +4 to blood, epidermal and muscle δ13C values in figures comparing with tooth dentin and collagen δ13C values. Repeat sampling of tooth and bone material from the same individual indicated a difference of 0.2‰ in δ13C and 0.3‰ in δ15N and therefore that these two tissues were directly comparable.
(b) Mitogenome sequencing and phylogenetic analyses
Mitochondrial DNA sequencing was performed to assess the degree of lineage sorting based on isotopic niche. A previous study had found a broad range of isotopic values within some North Atlantic lineages . However, these lineages were defined using mtDNA control region sequence data. Subsequent studies using complete mitogenome sequences from a global dataset of killer whales have shown the control region to be a poor indicator of phylogenetic relationships, while the coding genes of the mitogenome evolve in a more ‘clock-like’ manner and provide high phylogenetic resolution [27,39]. We therefore sequenced complete mitogenomes to improve phylogenetic resolution.
DNA was extracted and purified as per reference . Blank extractions were included to every five samples to monitor for contamination. To generate mitogenome sequences from recent (less than 200 years old) and ancient (greater than 1000 years old) samples, we used multiplexed DNA target enrichment hybridization capture coupled to high-throughput sequencing, using a published protocol  with a minor change, using heat instead of NaOH to release the captured DNA from the beads. Genomic DNA extract was built into blunt-ended libraries following Meyer & Kircher . Illumina sequencing libraries were built on the DNA extracts using NEBNext (Ipswich, MA, USA) DNA sample prep master mix set 1. Libraries were subsequently index amplified for 15–20 cycles using Phusion high-fidelity master mix (Finnzymes, Thermo Scientific, Finland) in 50 µl reactions following the manufacturer's guidelines. The libraries were then purified using MinElute PCR purification kit (Qiagen, Hilden, Germany) and pooled equimolarly in to a total of 2 µg, at which point they were ready for subsequent target enrichment.
To generate the bait for target enrichment, high-quality killer whale DNA extract was amplified into four overlapping long-range PCR products using primers LR2.1, LR2.2, LR3 and LR4 encompassing the whole mitogenome, following Morin et al. . Subsequently, these amplicons were converted into bait following Maricic et al. , after which target enrichment proceeded on the pooled libraries following Maricic et al. . The DNA concentration of the library eluted postcapture was measured using a 2100 Bioanalyzer (Agilent Technologies, CA, USA), then sequenced in subpartitions of single channels on an Illumina HiSeq 2000 platform using SR 100 bp chemistry.
Illumina HiSeq 2000 reads were filtered with AdapterRemoval , to remove adapter dimers as well as low-quality stretches at the 3′ ends. Filtered reads were then mapped to a reference killer whale mitogenome (GU187176.1) using BWA v. 0.5.9 , requiring a mapping quality of greater than 25. Clonal reads were collapsed using the rmdup program of the Samtools (v. 0.1.18) suite . Ambiguously mapped reads were also filtered out using Samtools and controlling for XT, XA and X1 tags. Consensus mitogenome sequences were then reconstructed using bam files, which were aligned in GENEIOUS (Biomatters Ltd) .
Phylogenetic relationships based on the sequence data were estimated using maximum-likelihood (ML) methods performed using webserver-based PHYML v. 3.0 , using the HKY + Inv + gamma model selected using jModelTest v. 1.1 . The transition/transversion ratio, the proportion of invariable sites, the gamma distribution and the starting tree, estimated using a BIONJ algorithm were also estimated by PHYML v. 3.0. The reliability of the optimized tree was estimated using 100 bootstrap replicates. To allow a heuristic visualization of lineage sorting based on trophic position, δ15N isotopic values were treated as continuous characters and their ancestral state for these characters was inferred using Mesquite v. 2.75 , with the ‘Trace Character Over Trees’ module applying the parsimony reconstruction method.
(c) Microsatellite genotyping and population structure analysis
To assess whether an evolutionary outcome of niche variation in this study system was assortative mating based on dietary preferences, we conducted an individual-based analysis of population structure. As DNA from degraded DNA sources can be prone to genotyping errors, for example allelic dropout, we used only high-quality DNA from modern skin biopsies for this analysis. DNA was extracted from biopsies of killer whales sampled either while feeding on fish or had stranded and found to have fish remains in their stomach, and from a group of six individuals taken by aboriginal hunters off Ammasalik, East Greenland, whose combined stomach contents included four harp seals and a hooded seal, but which also had worn teeth as observed in fish-eating individuals (see the electronic supplementary material, figure S1 for a map of sample locations). Individuals were screened with 12 polymorphic microsatellite loci (D08, D22, EV1, EV37, FCB4, FCB5, FCB11, FCB12, FCB17, KWM2a, Ttru04 and Ttru11) following the methods, loci and QC approaches, outlined in Foote et al.  and compared with previously published genotypes . A Bayesian model-based clustering algorithm performed by STRUCTURE v. 3.2.4  was used to infer population structure and probabilistically assign individuals to K clusters minimizing Hardy–Weinberg disequilibrium between loci within groups, without a priori knowledge of population units and limits. A series of five replicate independent runs were conducted for each value of K, set between one and five, using the correlated allele frequencies and admixture models. Each run used 106 iterations after a burn-in length of 105 iterations. To check for convergence of the Markov chain Monte Carlo, we compared the consistency of the results of the five replicates at each value of K.
4. Results and discussion
(a) Niche variation determined from δ15N and δ13C
The δ15N values of both subfossil and recent samples spanned a range of approximately 8‰, matching the difference in mean δ15N values of known prey items for North Sea killer whales (figure 2). While the range of δ15N values appears to be similar between time periods, there is a greater spread of δ13C values in the recent samples (figure 3). SEAB calculated using Bayesian inference and a corrected value for SEAC indicate a significant (p < 0.01) increase in isotopic niche width from SEAB = 5.4‰2 (SEAC = 5.7‰2) for the subfossil samples to SEAB = 11.0‰2 (SEAC = 11.6‰2) for the recent samples (figure 3). For comparison, the SEAC based on published isotopic values measured from skin of killer whales sampled in the Norwegian fjords when they are seasonally specializing on herring [29,31] was 0.62‰2 (figure 3); a significantly (p < 0.01) narrower niche width than both the ancient and recent North Sea samples. The 50, 75 and 95% Bayesian credible intervals based on 106 re-samplings are shown in the electronic supplementary material, figure S2. There was overlap of isotopic niche between sampling periods; 3.1‰2 of the 5.7‰2 SEAC of the ancient subfossil overlapped with the SEAC of the modern samples. Therefore, isotopic niche variation has been present in North Sea killer whales over timescales approximating to between 200 and 450 killer whale generations, as estimated from the age of the subfossil samples and assuming a generation time of 15 years .
Given the large geographical ranges across which killer whales typically travel, for example some North Sea killer whale groups have been photographically recaptured off Iceland , it is likely that some isotopic variation arises from geographical differences in foraging locations and that these have shifted during the Holocene. The mean δ15N of Atlantic herring ranges from 11.8 off of Iceland to 13.0 in the North Sea . Consequently, a previous study measuring isotope values of killer whale skin biopsies found that there are significant differences in δ15N between Iceland and Norway . This difference in mean δ15N between the two locations was approximately 2‰, much less than the variation in δ15N values found in this study. The samples all originated from the Holocene: a period when there was relatively little variation in δ13C of the baseline carbon source in the North Atlantic . There was no significant temporal trend in isotopic values between 1865 and 1995, despite changes in agricultural practices that could have led to a shift in baseline nitrogen values (see the electronic supplementary material, figure S3). Therefore, temporal and geographical variation in baseline isotopic values do not appear to explain much of the observed variation in isotope values among specimens in this study.
Among-individual differences in the diet likely explain a large part of the observed isotopic variation . This is consistent with the stomach contents of individuals, wherever available. For example, the two recent specimens from Denmark in figure 1 had stomach contents of either marine mammal remains or fish remains, consistent with their disparate isotopic values indicating differences in mean trophic position. Isotopic niche values based on measurements from tooth and bone tissue represent the cumulative intake of different components of the diet over the years that the sampled bone material and teeth growth layers were deposited [32–35]. Therefore, differences among individuals in isotopic niche indicate there is interindividual variation in the mean proportions of different prey items within the diet, but these individual preferences may not be discrete, and we do not have a measure of how specialized each individual is. Measurement of isotopic values from different growth layers within a tooth may provide an indicator of within-individual variation, e.g. . In the absence of this within-individual data, we refer to these among-individual differences in the long-term mean uptake of prey items, inferred from differences in isotopic niche, as relative specialization. The tooth wear, which is prevalent in the specimens used in this study (figure 1), suggests some overlap in either the diet and/or foraging method .
(b) Lineage sorting of mitogenome sequences based on isotopic values
We generated 29 new mitogenome sequences (GenBank accession numbers KF418372–418393) including nine from recently collected epidermal samples, 16 from tooth and bone specimens in museum collections (less than 200 years old) and four from subfossils (greater than 1000 years old). Following QC and filtering, mean sequence coverage ranged from 9 to 184× coverage and only nucleotide positions with a read depth of greater than 5× coverage were included. Consensus sequences from Illumina reads matched 100% with Sanger sequenced fragments from the same sample. The generated sequences resulted in the discovery of 22 new mitogenome haplotypes, all nesting within clades identified by previous analyses of 144 globally distributed modern samples [27,49] (see electronic supplementary material, figure S4).
The degree of lineage sorting of mitogenomes provides a proxy of the stage of speciation and insights into the transmission over time of foraging preferences. Tracing δ15N and δ13C isotopic values as continuous character traits over the ML phylogeny suggested that there has been multiple diversifications in isotopic niche (figure 4). There was also an indication of relatively stable transmission of isotopic niche along matrilineal lines within some clades, in particular those that were dominated by samples from Norway. However, it should be noted that these Norwegian samples were collected over a small geographical area, over a short-time span and were skin samples, therefore giving a short-term dietary signature. There were clear outlier individuals within some clades, typically those containing tooth and bone samples from the North Sea. For example, the second clade in the phylogeny contains two samples with very similar and relatively high N15 values, one from Shetland and one from the Kattegat that had marine mammals in its stomach. In the same clade are an individual from the Humber Estuary in England that had an isotopic signature consistent with a mainly pisciverous diet and an individual biopsy sampled from a pelagic trawler in the North Sea while it was feeding on mackerel escaping from the nets. Given the short branch lengths, this incomplete lineage sorting is consistent with relatively recent divergences in niche by some of the sampled lineages.
(c) Evidence for assortative mating
A predicted evolutionary outcome of long-term niche variation is assortative mating and potentially, speciation [8,14]. We tested whether there was evidence for assortative mating between a group that had fed on mammals prior to sampling and individuals that had been feeding on fish when sampled. The Bayesian model-based clustering algorithm performed by STRUCTURE v. 3.2.4  assigned the group from Greenland with seals in their stomach unambiguously to a population consisting mainly of herring-eating killer whales distributed from Norway to Iceland, across a range of population estimates from K = 2–5 (figure 5). Incorporating multiple individuals from the same group can bias the analysis and misinterpret social structure as additional population structuring owing to the close relatedness and sharing of alleles within the group . However, even at K = 5 (not shown), the seal-eating group was not distinguished from the other individuals within this predominantly fish-eating population. Therefore, although based on samples from only one group, our results do indicate panmixia (of neutral nuclear DNA markers) between at least some groups that feed on fish and some groups whose diet includes seals.
(d) Generalism, trade-offs and niche variation
Generalist individuals are expected to have higher fitness than specialists as they have access to a wider range of resources, and individual niche width is expected to expand to match the population niche width, but only if there are no biomechanical, physiological, cognitive or other constraints that restrict the variety of different prey resources an individual can efficiently consume [6,7,54]. In nature, these constraints appear to be common and consequently, niche width of generalist populations has been found to be owing to ecological differences among relatively specialized individuals . These constraints may also lead to reduced hybrid fitness and promote assortative mating .
One potential constraint of foraging for all available prey resources in killer whales is that different hunting strategies are required to efficiently hunt either marine mammals or fish [24,26,56–60]. The foraging strategies of mammal-eating killer whales are adapted to hunting large, single prey items with acute hearing, and they therefore typically travel in small groups, with a mean group size of three to five animals and forage in silence [24,26,56–60]. By contrast, the foraging strategies of fish-eating killer whales are adapted to hunting smaller, clustered prey with poor hearing, and therefore form a wide range of group sizes [24,26,60] and typically vocalize at a higher rate than mammal-eating killer whales [56–59]. However, the plastic nature of these phenotypic traits associated with each foraging strategy allows switching between strategies . As noted above, the presence of tooth wear suggests that there is some overlap in the diet of North Sea killer whales. Additionally, if phenotypic traits under ecological selection are transmitted through matrilineal cultural inheritance alone and are not genetically heritable, then there would be no cost to gene flow between relative specialists within these sympatric North Atlantic lineages.
Our results have implications for understanding speciation. Our study system provides a useful example of long-term niche variation and its evolutionary outcome. The shallow ecological cline indicated by isotopic values and apparent overlap in prey items consumed, in addition to the incomplete lineage sorting of mitogenomes (based on isotopic niche) and lack of evidence for assortative mating suggests that any progress towards speciation is still at an early stage in this system . This is in contrast to killer whale ecotypes in the North Pacific which fall either side of an ecological step, and for which complete lineage sorting of mitogenomes and strongly bimodal genotypic clustering indicate that speciation is at a late stage or has reached completion . There is evidence that sympatry in North Pacific ecotypes arose from secondary contact following a long allopatric phase . This proposed period of separation would have allowed the divergence of both genetically and culturally heritable traits, and reproductive isolation may have been consolidated upon secondary contact through reinforcement . The pattern of diversification in killer whales therefore reflects similar findings in other taxa, where progress along the speciation continuum varies in rate and directionality resulting in variation among study sites in the development of discrete or continuous phenotypes [8,21,62].
While our study is the first to investigate niche variation from such a lengthy temporal perspective, our results add to a growing body of literature that suggest sympatric speciation owing to niche variation is difficult to achieve. For example, niche variation has led to assortative mating but no associated phenotypic divergence in threespine stickleback Gasterosteus aculeatus , and niche variation has led to weak divergence of adaptive phenotypic traits (beak morphology) but no association with neutral genetic markers in Darwin's medium ground finch Geospiza fortis .
Perhaps the ecological differences among individuals in these empirical studies are not large enough to drive disruptive selection and/or assortative mating in the rapid manner predicted by some theoretical models . Additionally, in the absence of a strong link between ecology and reproductive isolation, sympatric speciation is predicted to be extremely slow .
This work was supported by a Marie Curie Actions ITN Fellowship ‘KWAF10’ grant and an exchange grant from the European Science Foundation research network ‘FroSpects’ awarded to A.D.F., Lundbeck Foundation grant no. ‘R52-A5062’, the Danish Basic Research Foundation ‘Geogenetics’ grant and the Natural Environment Research Council.
We thank all those who provided samples to this study and previous studies, in particular Henry van der Es of the Natuurhistorisch Museum Rotterdam, Kristian Gregersen and Kim Aaris-Sørensen of the Natural History Museum of Denmark and the owners of HG 333 Isafold, Captain Karsten Mølgård and crew for allowing us to collect biopsies during fisheries. We thank Kristian Gregersen and Kim Aaris-Sørensen of the Natural History Museum of Denmark, Copenhagen, Henry van der Es of the Natuurhistorisch Museum Rotterdam, Lillian Petersen, Cecilie Demring Mortensen and Kim Magnussen at the National High-throughput DNA Sequencing Centre, University of Copenhagen for their help with this study. We also thank two anonymous reviewers for their comments, which helped improve this manuscript.
- Received June 8, 2013.
- Accepted July 22, 2013.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.