“Type D” killer whale genomes reveal long-term small population size and low genetic diversity

Abstract

Genome sequences can reveal the extent of inbreeding in small populations. Here, we present the first genomic characterization of type D killer whales, a distinctive eco/morphotype with a circumpolar, subantarctic distribution. Effective population size is the lowest estimated from any killer whale genome and indicates a severe population bottleneck. Consequently, type D genomes show among the highest level of inbreeding reported for any mammalian species (F_ROH ≥ 0.65). Detected recombination cross-over events of different haplotypes are up to an order of magnitude rarer than in other killer whale genomes studied to date. Comparison of genomic data from a museum specimen of a type D killer whale that stranded in New Zealand in 1955, with 3 modern genomes from the Cape Horn area, reveals high covariance and identity-by-state of alleles, suggesting these genomic characteristics and demographic history are shared among geographically dispersed social groups within this morphotype. Limitations to the insights gained in this study stem from the nonindependence of the 3 closely related modern genomes, the recent coalescence time of most variation within the genomes, and the nonequilibrium population history which violates the assumptions of many model-based methods. Long-range linkage disequilibrium and extensive runs of homozygosity found in type D genomes provide the potential basis for both the distinctive morphology, and the coupling of genetic barriers to gene flow with other killer whale populations.

Tuhinga whakarāpopoto (abstract) te reo Māori

E mihi ana mātou ki a Tangaroa rāua ko Hinemoana, ko te kāinga o te kākahi/maki (ko ētahi o ngā kupu Māori mō te "killer whale", Orcinus orca). Waihoki, ka tika te mihi ki a Tūmatuenga, nā te mea, ki ōna wahi e pae ai ngā ika moana katoa. Ka mihi hoki mātou ki ngā whānau, ngā hapū, me ngā iwi e tiaki i tēnei taonga. Hoki atu ki te kākahi, nāna te moana hei rangatira me te aro o te rangahau nei. Mā te whakamātau i ngā raupapa ira o ngā tātai whakapapa o tēnei ika moana, e hura i te maha o ngā nohonga whakamoe kei roto i tēnei rāngai whāiti. Ko ngā whāinga rangahau e hora nei, ko te whakaraupapatanga tuatahi o ngā ira whakapapa o te momo kākahi ‘Type D’. Ko tōna āhuahanga he kiri kōheru, he mangu. Ko ngā moana matāo o Te Tiri o te Moana tōna kāinga. Mā te whakamātau i ngā raupapa ira nei, kāhore e huhua ana te maha o tēnei momo kākahi, hāunga e ērā atu momo. Koia e tohu ana i te iti o ngā uri o te rāngai ‘Type D’. E tohu ana ngā raupapatanga ira, he mea whakamoe tēnei rāngai kākahi ki a rātou anō (FROH ≥ 0.65). Ko ēnei ngā inenga whakamoe nui rawa puta i ngā uri katoa o ngā kararehe whāngote. Ko tēnei te tohu e tūpono ana ko tēnei iwi kākahi (type D) nei te tokoiti rawa o ngā momo kākahi kua rangahaua e te tangata, puta katoa i te ao. Kua tātari whakahāngai ngā raupapa ira o tētahi whakapapa kākahi (type D) o Paraparaumu (i mate whakauta i te tau 1955), ki ngā kāwai e toru o Cape Horne. Ahakoa he motuhake, he noho wehe ngā kāwai e rua nei, he hanga taurite te āhuaranga o ngā tātai heke, me nga tātai hono i waenga i ngā ira o ngā kāhui kākahi nei. He uaua rawa te kite i te motuhaketanga o te kāwai heke o ngā kākahi o Cape Horne nā te nui o ngā uri kua uhia ki te ao. He tohu pea tēnei i heke mai ēnei kāwai rerekē i tātahi puna puta-tahi. Engari rā, he nui noa atu nga āhuatanga hītori o enei kākahi e noho pohēhē ana, i runga i te uaua o ngā mahi hei whakapono i aua tatau. Koia pea tētahi āhuatanga i ngoikore ai ngā kōhanga tatau o te rangahau nei. Ahakoa kei te kākahi ‘Type D’ nei ōna ake āhuatanga, me tōna ake motuhaketanga, e kite hoki ana i ngā hononga-a-ira, ngā hononga kikokiko i waenganui i ngā kāwai whakapapa o ērā atu momo kākahi. He tohu tēnei i heke mai rātou i te tūpuna ōrite. Ahakoa tērā, kei te ngoto tonu ētahi tūahua hei ārai, hei aukati i te hono me te toha i waenga i ngā kāwai whakapapa nei.

Graphical Abstract

effective population size, inbreeding, killer whale, morphotype, runs of homozygosity

Introduction

Killer whales (Orcinus orca) are currently considered a single polytypic species (Committee on Taxonomy, 2021), although 5 distinct forms have been described from the Southern Ocean alone. These different forms, known as ecotypes and/or morphotypes (hereafter eco/morphotypes), have distinct habitat and prey preferences, as well as type-specific coloration and body sizes (Pitman and Ensor 2003; Pitman et al. 2007, 2011; Durban et al. 2017). Analyses of whole-genome sequences have shed new light on the demographic and evolutionary processes underlying the formation of eco/morphotypes (Foote et al. 2016, 2019, 2021). These include genetic adaptation to ecological niches and signatures of founder effects (i.e. much of the variation within each type coalesces in a recent common ancestor; Foote et al. 2016, 2019, 2021). However, there is variation among killer whale eco/morphotypes in their demographic histories. For example, the sympatric mammal-eating “transient” (also known as the Bigg’s killer whale) and fish-eating “resident” ecotypes found in the North Pacific are strongly differentiated across the genome and share ancestry with different outgroup populations (Foote and Morin 2016; Foote et al. 2016, 2019, 2021). These patterns are most parsimonious with the resident ecotype and transient ecotype each descending from distinct ancestral lineages, and the present-day sympatry between the 2 ecotypes resulting from secondary contact (Foote et al. 2019). In contrast, the Antarctic eco/morphotypes, types B1, B2, and C, despite their ecological and morphological differences, and at least partial geographic separation (Pitman and Ensor 2003; Pitman et al. 2007; Durban et al. 2017), share much of their ancestry and are only weakly genetically differentiated from each other (LeDuc et al. 2008; Morin et al. 2015; Foote et al. 2016, 2019, 2021). There is some evidence for erosion of genetic differentiation between ecotypes over time through sporadic geneflow, potentially via intermediary populations (Hoelzel et al. 2007; Foote et al. 2019, 2021).

Here, we report the results of the first nuclear genomic analyses on arguably the most distinctive killer whale morphotype described to date, the sub-Antarctic “type D” killer whale (Pitman et al. 2011). Type D killer whales have a noticeably bulbous head that appears more similar to a pilot whale than to other killer whale types, and a reduced eyepatch, that may even be absent in some individuals (see Fig. 1; Pitman et al. 2011). The preference of the type D killer whale for the offshore waters of the Southern Ocean makes at-sea sightings and on-shore stranding rare events (Pitman et al. 2011, 2020). The turbulent latitudes of the Roaring 40s and Furious 50s also hamper efforts to collect biopsy samples for DNA extraction (Pitman et al. 2011, 2020). Due to these challenges, type D killer whales have (up until this study) eluded biologists collecting tissue samples in the field, and thus, also those conducting genomic analyses in the lab. With these first fresh tissue biopsies we reconstruct the demographic history of this morphotype, and in doing so provide new insights into the processes that potentially underlie the unique morphology of type D killer whales.

Fig. 1.

The known and inferred distribution of type D killer whale in the Southern Ocean based on the original stranding record from New Zealand in 1955 and at-sea sightings. The latter are mostly observations of depredating animals around fishing vessels (e.g. southern Indian Ocean) and reports from tourism vessels (e.g. around the Scotia Arc, off the tip of South America). All sightings (except the NZ stranding) fall between 46°S and 62°S (total 48 listed), shown as a white, circumpolar band. Top photo is of the 1955 New Zealand stranding, (provided by A. van Helden); middle photos compare type A (photo credit: R. Pitman) and type D killer whales (photo credit: P. Tixier); lower photo is of one of the type D killer whales biopsy sampled in 2019 from which the genomes used in this study were sequenced (photo credit: J. Totterdell).

Materials and methods

Sample collection

Following up on information from local fishermen that type D killer whales were regularly depredating their fishing lines in Chilean waters along the continental shelf edge south of Cape Horn, our field team visited the area and encountered a group of approximately 30 individuals in January 2019 (Pitman et al. 2020). From this 1 group, photo-identification data, underwater video, and, key to this present study, biopsy samples were collected from 3 individuals (hereafter referred to by their NOAA Southwest Fisheries Science Center lab ID numbers: 202035, 202036, and 202037). Biopsy samples were taken using a 67-kg draw weight crossbow, firing floating darts with cutting tips 3.5 cm long and 0.6 cm in diameter. Samples were transferred to the SWFSC under CITES import permit (US774223/9) and export permit (19CL000006WS).

Modern DNA extraction, library build, and sequencing

DNA was extracted from skin biopsies preserved in salt saturated 20% dimethyl sulfoxide using the NucleoMag Tissue extraction kit (Macherey-Nagel, Allentown, Pennsylvania, USA). Genomic libraries were prepared using a blunt-end ligation method with double indexing (Meyer and Kircher 2010; Kircher et al. 2012). The following steps were done by Azenta Life Sciences commercial sequencing service. DNA libraries were quality checked using a D1000 ScreenTape on the Agilent TapeStation (Agilent Technologies, Palo Alto, California, USA), and were quantified using Qubit 2.0 Fluorometer. The DNA libraries were also quantified by real time PCR (Applied Biosystems, Carlsbad, California, USA). The samples were pooled (equimolar) and clustered on 2 lanes of a flow cell on the Illumina HiSeq instrument (4,000 or equivalent) according to the manufacturer’s instructions. The samples were sequenced using a 2 × 150 paired-end configuration. Image analysis and base calling were conducted by the HiSeq Control Software (HCS) on the HiSeq instrument. Raw sequence data (.bcl files) generated from Illumina HiSeq were converted into FASTQ files using Illumina’s bcl2fastq 2.17 software. Libraries for samples 202035 and 202036 were repooled and sequenced on an additional HiSeq lane to increase average depth of coverage to over 15× each.

Museum sample DNA extraction, enrichment capture, library build, and sequencing

On the 13 May 1955, a group of 17 killer whales stranded on Paraparaumu Beach, New Zealand—Aotearoa (Baker 1983). The distinctive shape and size of the eye patch of these killer whales, when compared with others in New Zealand (see Fig. 1), was highlighted by Visser and Makelainen (2000). These features were then found in photographs of animals at sea, resulting in the first formal description of the type D morphotype (Pitman et al. 2011). The museum specimen #1077 from the Museum of New Zealand—Te Papa Tongarewa, Wellington originated from the 1955 mass stranding of type D killer whales (Foote et al. 2013). Sequencing of nuclear genomic DNA from the New Zealand stranded specimen offers the opportunity to compare 2 social groups of type D killer whale. This is important as killer whale social groups are typically matrifocally philopatric and therefore comprised of close kin (Ford 2009).

Under Te Tiriti o Waitangi (the Treaty of Waitangi) Māori have sovereignty (rangatiratanga) over taonga (treasured) species and the natural environment. The Tiriti carries the obligation for scientists to resource consultation with relevant iwi (tribes) about the potential utilization of taonga species samples. In this case, the sample was exported and sequenced in 2016, before these important conversations occurred. We therefore sought guidance on how to redress this oversight and maintain the rights of Ngāti Toa Rangatira, the iwi from the area where the sample was taken, with regard to research on this taonga specimen.

DNA was extracted from a piece of dried soft tissue on the skull and from approximately 0.5 g of powdered tooth from the New Zealand stranded specimen. DNA extraction and purification from powdered tooth followed the protocols set out in Foote et al. (2013). Extraction and purification of DNA from the dried tissue was conducted using the Qiagen DNeasy kit following the manufacturer’s guidelines. Four individually indexed DNA libraries (2 for each tissue type) were built on 5 µL of extract each using the blunt-end single-tube (B.E.S.T) library build method (Carøe et al. 2018) and individually amplified for 20 cycles with an annealing temperature of 60 °C using AmpliTaq Gold DNA Polymerase (Applied Biosystems). Libraries were subjected to 2 rounds of whole-genome enrichment capture using genome-wide biotinylated RNA baits (2 reactions) built from modern killer whale DNA by Arbor Biosciences (see Enk et al. 2014), with a hybridization period of 24 h at 55 °C. Libraries were then pooled and 100 bp single-end read sequenced by the Danish National High-Throughput DNA Sequencing Centre in 1 lane of an Illumina HiSeq2500 sequencer. Blank extracts and library builds were included in this pipeline to monitor for contamination.

Mapping and filtering of sequence data

Reads from each individual were processed with AdapterRemoval2 (Schubert et al. 2016) to trim residual adapter sequence contamination and to remove adapter dimer sequences as well as low-quality (Q < 30) stretches at the read ends. Filtered reads >30 bp were then mapped using BWA-MEM (Li 2013), requiring a mapping quality greater than 30, to both the original killer whale reference genome Oorca1.1 (GCA_000331955.2) (Foote et al. 2015), and to a high-quality, chromosomal O. orca reference assembly generated from a North Pacific transient killer whale (Kardos et al. in press; CNGBdb accession: CNP0002439). The latter was used in analyses that required high contiguity, such as the estimation of runs of homozygosity (ROH). To assess the bias generated by the reference genome, analyses were rerun using an updated high-quality, chromosomal O. orca reference assembly mOrcOrc1.1 (GCA_937001465.1) (Foote and Bunskoek 2022) generated from a Norwegian killer whale. Clonal reads were collapsed using the markdup function of SAMtools v1.13 (Li et al. 2009). Repetitive elements were identified using RepeatMasker v4.1 (Smit et al. 1996) and the Cetartiodactyl repeat library from RepBase (RepeatMaskerEdition-20181026; Kohany et al. 2006) and then masked using BEDtools v2.26 (Quinlan and Hall 2010). Coordinates of repeat regions in BED format are provided in the Dryad repository (https://doi.org/10.5061/dryad.wpzgmsbr7). Only autosomal chromosomes were included, selected using SAMtools view. Patterns of postmortem DNA damage in the genome of the type D museum specimen were estimated and visualized using mapDamage2 v2.1 (Jónsson, et al. 2013).

Additionally, reads were mapped using BWA-MEM to the previously generated mitogenome sequence of the New Zealand type D specimen (accession KF164610; Foote et al. 2013). In a previous study of 139 killer whale mitogenome sequences (Morin et al. 2010), the inclusion of intra- and inter-lab PCR, library build, and sequencing replicates identified inconsistencies in the assembly of polynucleotide repeat regions: one of between 9 and 14 Cs in a row (positions 1130 to 1144 in the original alignment), and another region of 7 to 8 As in a row (positions 5210 to 5217). Morin et al. (2010) therefore shortened these to a fixed set of 9 Cs and 7 As, respectively, to avoid introducing potentially erroneous variation into phylogenetic analysis in that and subsequent follow-up studies (e.g. Foote et al. 2013; Morin et al. 2015). We follow this same conservative approach in our mapping of mitogenomes and identification of sequence variation in this study.

Inferring matrilineal relationships

As the mitochondrial genome of the New Zealand stranding had previously been sequenced (Foote et al. 2013), our first goal was to compare the mitogenome sequences of the Chilean samples with the New Zealand sequence to better understand matrilineal diversity.

The matrilineal relationships among the 3 biopsy sampled type D killer whales, the New Zealand type D sample, and a global dataset of 452 individuals were estimated from the mitogenome alignment using maximum-likelihood methods performed using PHYML v. 3.0 (Guindon et al. 2010), using the HKY + Inv + gamma model selected using jModelTest v. 1.1 (Posada 2008). 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.

Inferring kinship

Given that killer whales typically associate in kin-based social groups, it was important to establish the kinship among the 3 Chilean type D killer whales, as close kin would reflect sampling of the same genealogies. The relationships among the 3 biopsied type D killer whales were inferred using IBSrelate (Waples et al. 2019), which can estimate pairwise relatedness when data are limited, i.e. when population allele frequencies are not available. Kinship estimators included R0 (Rosenberg 2006), R1 (Waples et al. 2019), and KING-robust kinship estimator (Manichaikul et al. 2010) derived relatedness measures, each of which compares the ratios of different identity-by-state (IBS) genotypes. Ancestral state was inferred for each site of the killer whale genome by mapping closely related outgroup species and calling the consensus allele. Outgroup species used were the long-finned pilot whale Globicephala melas (accession number: SRR8861593), bottlenose dolphin Tursiops truncatus (accession number: SRR2148845), and white-beaked dolphin Lagenorhynchus albirostris (accession number: SRR11679513).

Estimating covariance and distance measures

To further understand the relationship between the New Zealand and Chilean samples, which offered the best opportunity to sample independent genealogies within type D killer whales, we compared the sharing of genome-wide alleles. Covariance and distance (IBS based) were estimated using ANGSD v0.933 (Korneliussen et al. 2014) from pseudo-haploid genotypes, i.e. randomly sampling an allele at each polymorphic site, to reduce bias from variation in coverage between the museum and modern (biopsy) samples. Only sites covered by uniquely mapping reads, with a minimum mapping quality of 30, adjusted to 50 for excessive mismatches, and minimum base quality of 20, adjusting q-scores around indels, a minimum allele frequency of 0.01, and a probability of being a variable site of P < 0.000001 were included. Transitions were excluded to reduce bias from postmortem DNA damage due to deamination of cytosine in the sequences obtained from the museum specimen (Lindahl 1996; Hofreiter et al. 2001). These filtering steps resulted in a pseudo-haploid genotype call at 5,960 transversions covered by at least 1 read in all samples including the museum sample. The majority (>98%) of single nucleotide polymorphisms (SNPs) were spaced at a physical distance of >10 kb apart (Supplementary Fig. 1). Pairwise distance between individuals was averaged across all included sites, where distance per site is 1 if different bases were called, and 0 if the same base was called, as per Korneliussen et al. (2014). Covariances were estimated as per Korneliussen et al. (2014).

To provide a frame of reference of within-population covariance and IBS, these analyses were repeated, comparing just 1 modern type D and the museum type D sample with published genomes (European Nucleotide Archive: ERS554443–ERS554448; Foote et al. 2016) from 5 killer whales from the small (<100 individuals) insular Southern Resident population of Washington State, USA (Parsons et al. 2009), which has a high incidence of inbreeding (Ford et al. 2011, 2018; Kardos et al. in press). Due to the low depth of coverage (~2×) data available for the Southern Residents, pseudo-haploid genotype calls were made at 2,233 transversions covered by at least 1 read in all samples including the museum sample.

All additional analyses were conducted on full depth of coverage autosomal chromosomes from the 3 biopsied Chilean killer whales. No further analyses incorporated the New Zealand sample.

Estimating theta (θ) from individual genomes

To understand the demographic history of type D killer whales and how this may have impacted genomic and potentially phenotypic variation, we first estimated a summary statistic of genome-wide variation. We estimated heterozygosity of autosomal regions from each of the 3 high coverage type D diploid genomes and made direct comparison to the globally sampled dataset of Foote et al. (2021). Under the infinite sites model (Kimura 1969), individual heterozygosity is a good, unbiased estimator of the population mutation rate, theta (θ; Watterson 1975). Theta was estimated directly from the filtered BAM files using the maximum-likelihood estimator (Lynch 2008) implemented in mlRho v.2.9. (Haubold et al. 2010).

Inferring changes in effective population size (N_e) from individual genomes

To better understand temporal changes in demography, we inferred N_e through time using the pairwise sequential Markovian coalescent (PSMC) (Li and Durbin 2011), a method suitable for studying ancient demographic history from a single unphased genome per population. The PSMC model estimates the Time to Most Recent Common Ancestor (TMRCA) of segmental blocks of the genome and uses information from the rates of the coalescent events to infer N_e at a given historical time, thereby providing a direct estimate of changes in past population size and structure (Li and Durbin 2011). A consensus sequence of the BAM file was then generated in fastq format using: first, SAMtools mpileup command with the −C50 option to reduce the effect of reads with excessive mismatches; second, 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, 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), with 100 bootstraps. The plot was scaled to time in years assuming a generation time of 25.7 yr and a mutation rate of 2.34 × 10⁻⁸ substitutions/nucleotide/generation (Foote et al. 2021). To perform bootstrapping, chromosomes were first split into shorter segments using the splitfa option. PSMC then uses a standard random sampling with replacement approach.

Calling ROH from individual genomes

ROH can further provide an indication of temporal changes in coalescent rates due to population size (Foote et al. 2021). ROH were identified using the window-based approach implemented in PLINK v1.07 (Purcell et al. 2007), as per Foote et al. (2021), incorporating the 3 type D killer whales with the Foote et al. (2021) dataset. Sliding window size was set to 300 kb, with a minimum of 50 SNPs at a minimum density of 1 SNP per 50 kb (across the dataset of 29 genomes) required to call an ROH. To account for genotyping errors, we allowed up to 3 heterozygote sites per 300 kb window within called ROHs, as per Ceballos et al. (2018). A length of 1,000 kb between 2 SNPs was required in order for them to be considered in 2 different ROHs.

Estimating genome-wide linkage disequilibrium profiles from individual genomes

A key impact of demographic bottlenecks is an increase in linkage disequilibrium (LD). To estimate LD, we calculated the correlation of zygosity (Δ) from the distribution of heterozygotes across all unmasked sites within the highest coverage type D genome (sample 202035) using mlRho v.2.9 (Lynch 2008; Haubold et al. 2010; Lynch et al. 2014). The 2-site zygosity classes comprise of H₀ (both sites are homozygous), H₁ (1 site is heterozygous), and H₂ (both sites are heterozygous). Estimation of Δ is based on the expectations of heterozygosity at the 2 sites due to recombination, drift and mutation (Strobeck and Morgan 1978; Lynch 2008; Haubold et al. 2010; Lynch et al. 2014), under the assumptions of the standard neutral model of evolution (Kimura 1968, 1983). Scaling Δ by the level of heterozygosity (θ) estimated for each genome by mlRho we can approximate LD (Δ/θ ≅ r²; Lynch et al. 2014), where r² is the measure of LD derived from the square of the correlation coefficient between the presence or absence of particular alleles at 2 different loci (Hill and Robertson 1968). We estimated Δ in steps of 1 bp between loci along the first 10 bps, in steps of 10 for the first 100 bps, in steps of 100 for the first 1,000 bp, in steps of 1,000 for the first 10,000 bp, and in steps of 10,000 for the first 100,000 bps of each chromosome.

Testing for introgression

The first genome sequences of type D killer whales provide an opportunity to test hypotheses regarding the distinct morphology, e.g. the “pilot whale-like” bulbous head shape (Pitman et al. 2011). Type D killer whales and long-finned pilot whales share an overlapping distribution (Pitman et al. 2011; Jefferson et al. 2015). Interspecific and intergeneric hybridization between a wide range of cetacean species is well documented and intrinsic barriers to interspecific introgression appear to be permeable to gene flow (Crossman et al. 2016). Given these observations, a simple tree-based test based upon excess sharing of derived alleles between 1 in-group population and an outgroup (Green et al. 2010) was applied to test whether similarity in morphological features of type D killer whales and long-finned pilot whales reflects ancestral introgression. The long-finned pilot whale sample was from an individual stranded on the sub-Antarctic Kerguelen archipelago, which is proximate (1,200 km from island to island) to the Crozet archipelago. Both the typical killer whale morphotype found in Crozet waters and the type D morphotype are observed directly interacting with the fisheries by depredating fish catches from long-line fisheries, both in Crozet as well as the Kerguelen EEZ, with the same type D killer whales individuals photo-identified in both localities (Tixier et al. 2016; Amelot et al. 2022). D-Statistics in which a “regular” Crozet killer whale and a type D killer whale were the ingroup, and a long-finned pilot whale the outgroup, with a haploid bottlenose genome representing the ancestral state at each site in the genome, were implemented in ANGSD (Korneliussen et al. 2014) and performed by sampling a single base at each position of the genome to remove bias caused by differences in sequencing depth at any genomic position. The definition used here is from Durand et al. (2011):

$${\displaystyle \begin{array}{l}{\displaystyle D=\frac{\mathrm{n}\mathrm{A}\mathrm{B}\mathrm{B}\mathrm{A}-\mathrm{n}\mathrm{B}\mathrm{A}\mathrm{B}\mathrm{A}}{(\mathrm{n}\mathrm{A}\mathrm{B}\mathrm{B}\mathrm{A}+\mathrm{n}\mathrm{B}\mathrm{A}\mathrm{B}\mathrm{A})}}\end{array}}$$

where nABBA is the number of sites where the Crozet killer whale (first ingroup) shares the ancestral allele with the dolphin (the ancestral state), and type D (second ingroup) and pilot whale (outgroup) share a derived allele; and, nBABA is the number of sites where the type D killer whale shares the ancestral allele with the dolphin, and the Crozet killer whale and pilot whale share a derived allele. Under the null hypothesis of no gene flow between the pilot whale and type D killer whale lineages, we expect an approximately equal proportion of ABBA and BABA sites and thus D = 0. The statistical significance of the deviation from 0 was assessed using a Z-score based on the assumption that the D-statistic (under the null hypothesis) is normally distributed with mean of 0 and a standard error estimated using a jackknife procedure.

Commands and code are available at https://github.com/AndrewDFoote/Orcanomics.

Results

Mapping and filtering of sequence data

Distribution of effective depth of coverage for each January 2019 Cape Horn sample is shown in Supplementary Fig. 2. The mean effective (i.e. postfiltering) coverage per sample was 21× for 202035, 17× for 202036, and 8× for 202037. All raw sequencing data used in this study are publicly archived under the Orcanomics BioProject (NCBI accession: PRJNA531206; Biosamples: SAMN20330599; SAMN20330592; SAMN20330591).

Sequencing of the New Zealand museum specimen resulted in sequence coverage of ≥1× at 25,811,700 bp of the 1.3 Gb repeat masked assembly. Genomic data from this specimen have been archived on the Aotearoa Genomic Data Repository (https://data.agdr.org.nz/), accession via https://repo.data.nesi.org.nz/discovery/TAONGA-AGDR00018/.

Analyses of potential misincorporations using mapDamage2 (Jónsson et al. 2013) revealed that sequencing reads exhibited characteristic postmortem damage patterns (Dabney et al. 2013), specifically an excess of C>T transitions at the 5ʹ termini as expected from deamination, and the complementary G>A transitions at the 3ʹ termini (Supplementary Fig. 3). Therefore, only transversions were considered in all downstream analyses that included the historical sample.

Matrilineal relationships

Mitochondrial genomes from the 3 type D killer whales sampled off Chile were all identical, indicating a shared matrilineal genealogy, consistent with the individuals being part of the same matrilineal kin-based social group (typical of killer whales).

Phylogenetic relationships based on the sequence data are presented in Supplementary Fig. 4. The mitogenomes of the 3 modern type D samples differed from the mitogenome from the type D New Zealand individual (previously characterized in Foote et al. 2013) by 4 nucleotide differences (Supplementary Fig. 4), equating to a similarity of 99.98% in a BLAST search. Three of the 4 differences between the modern and museum type D haplotypes occur in the hypervariable region, 2 of which are at potential mutational hotspots, and none of the alternative nucleotides at any of these sites are private to either type D haplotype (Supplementary Fig. 4). These details, together with the identical hypervariable region sequences obtained by both Sanger and high throughput sequencing (Foote et al. 2013), give confidence in the museum haplotype sequence, just as the identical sequences obtained from the 3 modern samples give confidence in the modern type D haplotype sequence. In summary, despite the small sample size, the modern and museum haplotypes indicate some variation in the mitochondrial DNA of type D, but which share a recent common matrilineal ancestor. An alignment of 159 mitogenome haplotype sequences from 455 individuals sampled globally (see Morin et al. 2015) and the new type D sequences is archived in the Dryad repository associated with this study (https://doi.org/10.5061/dryad.wpzgmsbr7).

Kinship

The 3 biopsied type D killer whales were all inferred to be female, based upon approximately equal mean depth of coverage of the longest autosome and the X-chromosome. Estimates of kinship statistics from genotype likelihoods in IBSrelate (see Waples et al. 2019) indicated close familial relationships among the 3 individuals, albeit with the caveat that IBSrelate assumes individuals are not inbred. IBS heterozygous genotypes increase R0 and R1 estimates, so full-siblings are expected to have higher R0 and R1 values than parent–offspring pairs. Therefore, we tentatively infer the 3 individuals represent a mother (202037) and her 2 full-sibling daughters (202035 and 202036; Table 1), but caution that lower R0 and R1 estimates due to reduced heterozygosity may also reflect the lower coverage of the genome of 202037. This potential close kinship means 2 of our sampled genomes (202035 and 202036) are not independent even for estimation of inbreeding from the previous generation, and all 3 will share approximately ≥50% of their genomic ancestry. Therefore, our results will reflect many shared genealogies. For this reason, the additional comparisons with the potentially unrelated museum specimen from the New Zealand stranding are particularly valuable.

Table 1.

R0, R1, and KING-robust kinship estimates from genotype likelihoods using IBSrelate (Waples et al. 2019) on pairwise comparisons of the 3 biopsied type D killer whale genomes.

Individual 1	Individual 2	Number of sites	Kinship estimator
			KING	R0	R1
202035	202036	9,679,128	0.263	0.115	0.957
202035	202037	9,469,799	0.288	0.056	0.878
202036	202037	9,552,730	0.276	0.071	0.852

Higher R0 and R1, and lower KING-robust values indicate higher relatedness. The number of sites included in each pairwise comparison is given in column 3..

Covariance and distance estimates

Consistent with the relatedness findings, the 3 type D killer whales sampled from the same social group off of Chile have highly correlated genotypes (Fig. 2). However, comparable high covariance and low genetic distance was also found between the New Zealand museum sample and modern type D genomes (Fig. 2). To fully put this into context, genetic covariance between the museum and modern type D samples is markedly higher, and genetic distance markedly lower, than estimates between contemporary Southern Resident killer whales: an insular population of <100 individuals (Supplementary Table 1; Ford et al. 2018). These findings strongly suggest limited genetic variation within the type D morph. We consider the alternative explanation that the 1955 New Zealand stranding was close kin of the group biopsied in 2019 as unlikely due to the time between sampling events, and that the modern samples originate from a different matriline to the museum specimen.

Fig. 2.

a) Distance and b) covariance matrices estimated from pseudo-haploid genotype calls at 5,960 sites with transversions (to avoid spurious C→T or G→A misincorporations resulting from postmortem DNA damage). Each matrix comprises a globally sampled dataset of modern killer whale genomes, 3 modern type D killer whales sampled off of Cape Horn, and the museum sample from the 1955 mass stranding of 17 individuals on Paraparaumu Beach, New Zealand (marked *). In panel (a), IBS increases from dark to light shading; in panel (b), covariance increases from light to dark shading. Samples are ordered based upon hierarchical clustering.

There is no clear strong affinity between type D and any other of the globally sampled populations included in this comparison (Fig. 2). A Principal Component Analysis (PCA) on the pruned covariance matrix constructed using 1 sample per population (Supplementary Fig. 5) largely recapitulated the previous inference of the relationships among populations (Foote et al. 2019). The type D genome apparently shares some of the variation in the Antarctic type B1 eco/morphotype segregating from North Pacific and North Atlantic killer whales along PC1. However, the differences between types D and B1 result in segregation along PC2.

Theta and effective population size

Estimates of theta (θ) generated by mlRho for each type D genome were 0.000178 to 0.000179 for 202035, 0.000174 to 0.000176 for 202036, and 0.000174 to 0.000175 for 202037. These values are between 0.1× and 0.5× the values of individual θ estimates from a global sampling of killer whale genomes (Foote et al. 2021). The closest previously published estimate was θ = 0.000370 for a killer whale from the critically endangered Scottish population, of which just 2 adult males have been sighted in the past 5 yr (Foote et al. 2021; Allendorf et al. 2022). As an additional comparison, we estimated θ = 0.000213 for the critically endangered vaquita (Phocoena sinus), from 20× mean depth of coverage short-read data (SRR15435925) mapped to the longest chromosome (CM018156.1; 185,845,356 bp) of the vaquita reference genome assembly (GCA_008692025.1, mPhoSin1; Morin et al. 2021). Genetic diversity in the type D genomes is thus more comparable to that found in the vaquita genome than to genomes of other killer whales. It is however important to remember that this statistic is based on the assumption that individual genomes are sampled from populations at demographic equilibrium. In reality, the action of linked selection is likely to directly impact estimates of demography and N_e.

The estimates of θ were to some extent supported by the estimated changes in N_e through time. PSMC samples local genealogies along the genome, estimating the TMRCA between alleles on each chromosome. The rate of coalescence at different time points reflects the population size at that time: many lineages will coalesce during periods of small population size, and conversely there will be less coalescence when population size is large. The trajectory of the PSMC plots for modern type D samples 202035 and 202036 (Fig. 3; Supplementary Fig. 6) suggest increased coalescence and small effective population size during the last glacial period. From the Last Glacial Maximum (LGM) to the start of the Holocene the plot levels off at a harmonic mean of N_e ≈ 150, roughly an order of magnitude less than the effective population size derived from our individual θ estimates using mlRho. This discrepancy may indicate an expansion in effective population size in type D during the Holocene, but this requires further investigation when more type D genetic samples are available.

Fig. 3.

PSMC plot of effective population size (N_e) through time, estimated from the high coverage genome sequence of a type D killer whale (202035). The thick red line represents the median estimate for this individual, and thin red lines represent 100 rounds of bootstrapping. The vertical dashed line indicates the approximate mean coalescent time (= 4N_e generations, assuming a generation time of 25.7 yr) for neutral mutations present at 10⁴ yr BP. Inferences further back in time from this point become increasingly stochastic, as illustrated by the variance among bootstrap values, due to the declining number of coalescent events from which to infer N_e. The plots of median estimates of N_e through time for individual genomes of Antarctic type B2, Chatham Islands and Norwegian killer whales are shown in shades of blue (see key) for comparison. These 3 individuals were selected as they reflect the 3 main trajectories of N_e through time found in a global dataset of 26 individuals (Foote et al. 2021).

The number of recombination events per each time interval was >20, which should avoid overfitting of the model (Vijay et al. 2018). However, the number of inferred recombination events was low in the modern type D sample (15,103 recombination events for sample 202035); which e.g. is an order of magnitude lower than in the Antarctic Peninsula type A killer whale genome (171,374 recombination events). Potentially, recombination events that result in no change in the genealogy, i.e. between recombination of identical haplotypes (see e.g. Deng et al. 2021) or those resulting in gene conversion rather than crossing-over (CO) (Wiuf and Hein 2000), will go undetected. This underestimation of recombination events due to CO of identical haplotypes would be expected to be higher in genomes from more inbred populations such as type D, and lower in admixed genomes such as the Antarctic Peninsula sample. These findings are consistent with low incidence of outbreeding in type D killer whales.

There are further caveats to consider when interpreting demographic inference from methods such as PSMC (see Beichman et al. 2017, 2018). The model assumes neutrality, when in reality purifying selection and its impact on linked sites (background selection or BGS; Charlesworth et al. 1993) can affect more than 95% of the genome and bias demographic inference (Ewing and Jensen 2016; Pouyet et al. 2018; Johri et al. 2021). Additionally, inference of population size at deep timescales depends upon there being sufficient loci within the genome that coalesce deeper back in time. Coalescent theory predicts such alleles will be relatively rare (Takahata and Nei 1985) in cases where the population has experienced a bottleneck (dependent upon the magnitude of the bottleneck), and only very few alleles within the genome will coalesce further back in time. The mean coalescent time for neutral mutations at a constant population size is expected to be 4N_e generations (approximately 15,000 yr in killer whales). Given the estimated N_e for type D at 10,000 yr BP (Fig. 3), the neutral mutations present at that time would have a geometric distribution of coalescent times with the mean at approximately 25,000 yr BP. Populations with declining N_e will retain some older variation reflective of the larger ancestral N_e. Consistent with these expectations, tracing the PSMC plot backwards in time from the LGM we see that inference of ancestral effective population size further back in time than 4N_e generations becomes increasingly noisy (indicated by the variation in bootstraps around the median), and then stops at 0.5 MYA (Fig. 3). We therefore caution that there are too few deep-time coalescent events to accurately infer effective population size and the deeper demographic history (>25 KYA) of this morphotype from coalescent-based approaches.

Lastly, all time point references are dependent upon the accuracy of the assumed mutation rate and generation time. Thus, we caution against a literal interpretation of the plot as an indication of precise changes in N_e through time, but rather as an indication that much of the variation in the genome of the type D killer whale coalesces in a recent (<25 KYA) ancestor.

Inbreeding

The 3 modern type D killer whale genomes had the highest median length of ROH (1.05, 1.23, and 1.19 Mb) among the 29 genomes analyzed (Fig. 4; Supplementary Fig. 7). The sum and number of ROH (SROH and NROH) are highly correlated among killer whale genomes (Foote et al. 2021). The 3 modern type D genomes are the greatest outliers of this correlation (Supplementary Fig. 8), due to harboring an excess of long ROH relative to the other killer whale genomes. ROH density along chromosome 1 highlights the extent and overlap in ROH among the 3 type D genomes (Fig. 4b). Regions in which no ROH are found in any of the 3 modern type D genomes are sparse and typically overlap with the absence of ROH in the global dataset (Fig. 4b), suggestive of such regions representing recombination hotspots (Kardos et al. 2017; Foote et al. 2021).

Fig. 4.

a) Kernel density (violin) plots of the length of ROH in a selection of killer whale genomes. White rectangle shows the interquartile range, and the black bar the median of the data. In addition to the genome of a type D killer whale (202035), a subset of genomes included in Foote et al. (2021) are shown for comparison, selected as they represented a range of demographic histories and genetic diversity and include the genomes with the lowest (Antarctic Peninsula) and highest (Scotland) genomic inbreeding. Antarctic type B2 had higher counts of ROH than type D, but these comprised mainly of short (<1 Mb) ROH, thereby differing from the pattern of autozygosity found in type D. A comparison with the full global dataset is shown in Supplementary Fig. 7. b) Physical mapping of ROH density along chromosome 1 of all ROH > 0.3 Mb in the global dataset used in Foote et al. (2021) and the 3 type D genomes. The black line in the top panel shows counts of overlapping ROH in nonoverlapping 100-kb windows for all genomes; the red line shows counts of overlapping ROH in the 3 type D killer whale genomes. The bottom panel presents ROH density per individual genome. Each 100-kb window containing an ROH is shaded red and blue in alternating rows.

Individual genomic inbreeding coefficients (F_ROH) were estimated from ROH as the proportion of the autosomes found in ROH >1.5 Mb, to allow direct comparison with published F_ROH estimates from the global dataset (Foote et al. 2021). Shorter ROH (≤1.5 Mb) are expected to coalesce deeper in the pedigree, e.g. during ancestral bottlenecks, and therefore subjected to more generations of purifying selection to remove deleterious recessive alleles than longer ROH. Consequently, short ROH are less likely to contribute to inbreeding depression (Szpiech et al. 2013; Stoffel et al. 2021), and hence it is common practice to remove short ROH from estimates of inbreeding coefficients. Genomes for the 3 modern type D samples 202035, 202036, and 202037 had estimated inbreeding coefficients of F_ROH = 0.65, 0.73, and 0.70, respectively. In comparison, the highest inbreeding coefficient previously estimated in the global dataset was F_ROH = 0.38, for the genome from a small Scottish population on the brink of extinction (Foote et al. 2021).

To draw comparison with the results of Kardos et al. (in press), we further estimated F_ROH including only ROH >1 Mb, and then ROH >10 Mb, and using the same mapping reference as per that study. Inbreeding coefficients of F_{ROH(>1 Mb)} = 0.75, 0.81, and 0.81 for 202035, 202036, and 202037, respectively, are higher than found in any of the 100 genomes from the endangered Southern Resident killer whale population (Kardos et al. in press). We estimated F_{ROH(>10 Mb)} of 0.04 for 202035 and 202037, and 0.06 for 202036, respectively. Whilst these values are indicative of ongoing high levels of inbreeding, they are comparable to that found in the Scottish genome (Foote et al. 2021) and Southern Resident killer whales (Kardos et al. in press). However, it should be noted that the methods for calling ROH differ between this study and Kardos et al. (in press), which may impact direct comparison. The length distributions of ROH in type D genomes suggest that the inbreeding coefficient in type D killer whales is a reflection of a long-term small population size as supported by the PSMC analysis, rather than an increase in inbreeding due to a more recent population decline during the Anthropocene.

Linkage disequilibrium

Demographic history, recombination, and admixture are key factors that influence LD along the genome (Pritchard and Przeworski 2001; Reich et al. 2001; Gray et al. 2009). Estimates of correlation of zygosity (Δ) scaled by θ are expected to represent LD (Δ/θ ≅ r²) and therefore reflect recombination rate (Lynch et al. 2014). The relationship among our estimates of scaled Δ at different distances between loci were as expected for the Antarctic Peninsula, Chatham Islands, Norway, type B2 and Scotland genomes (Fig. 5). These genomes were selected as they represented a range of demographic histories and genetic diversity (Foote et al. 2021). The results are qualitatively similar to those found in a comparison of Denisovan archaic human, and Chinese and African modern humans, and nonhuman primates in Figs. 4 and 5, respectively of Lynch et al. (2014). All the above genomes show parallel nonlinear declines in scaled Δ with increasing physical distance. The nonlinearity of this relationship is expected due to recombination events that result in gene conversion rather than crossing over (Lynch et al. 2014); gene conversion tracts are relatively short (Wiuf and Hein 2000) and so disproportionately (compared with CO recombination events) impact closely linked sites.

Fig. 5.

The relationship between the correlation of zygosity (Δ) scaled by estimates of theta (θ) with physical distance between sites for 6 killer whales from populations with different demographic histories, shown for Chromosome 1. Ordering in the legend reflects effective population size N_e, from smallest (top) to largest (bottom). Plots (a) and (b) show the same data plotted on a logarithmic and linear y axis plotted, respectively, the latter to include negative values.

The estimates of scaled Δ at any given physical distance between 2 loci reflects both the demographic histories (based on long-term N_e, see Foote et al. 2021) of the populations represented by each genome, and the inbreeding of the individual. For example, the Scottish killer whale, which has the lowest long-term N_e (omitting the type D), has the highest scaled Δ at any given distance, whereas the outbred Antarctic Peninsula killer whale has the lowest scaled Δ. This meets the expectation that LD will increase with lower long-term N_e (Slatkin 2008). All the killer whales included in the analysis presented in Fig. 5 are arguably sampled from nonequilibrium populations (Foote et al. 2021). However, despite violations of model assumptions, which include changes in effective population size and probable genome-wide selection (BGS), the method appears to be robust and able to provide plots that meet the expectations of how linkage between 2 loci would vary both with physical distance along the genome, and differences in long-term N_e. This is also evident in the Denisovan and western lowland gorilla (Gorilla gorilla gorilla) genomes included in Lynch et al. (2014), which are sampled from nonequilibrium populations.

The observed correlation in zygosity in the type D genome apparently deviates sufficiently from expectations under the standard neutral model of evolution to produce negative maximum-likelihood estimates of Δ. Scaled Δ starts 2 orders of magnitude higher than in the Scottish killer whale genome at distances between sites up to 10 bp; at distances between sites greater than 10 bp, estimates of Δ become negative (and are therefore not shown on a log-axis; Fig. 5). Δ is a measure of the degree to which the distribution of heterozygous sites deviates from expectations under random segregation. Demographic changes from a large ancestral effective population size to a severe population bottleneck, combined with extensive BGS would be expected to increase the zygosity class H₀, whilst reducing H₁ and H₂ (see methods) below the expectations of a population in approximate drift–mutation–recombination equilibrium. Additionally, extent of homozygous tracts in type D genomes will reduce detection of CO recombination events and is likely exacerbating the deviation from model expectations.

Thus, while this analysis does not provide a quantitative measure of LD in the type D killer whale genome, it does provide a qualitative indication of the extent to which the demographic history of type D killer whales deviates from the standard neutral model of evolution (relative to the other killer whales analyzed with this method).

Introgression

The D-statistic (Green et al. 2010) considers a tree-like relationship among 4 populations, in this case (((Crozet, type D), pilot whale), dolphin), and estimates whether either killer whale shares an excess of derived alleles with the pilot whale. Significant sharing of derived alleles between an ingroup sample and the pilot whale could indicate introgression from the pilot whale into that population. Estimation of D(Crozet, type D; pilot whale, dolphin) ± 1 SE = −0.0935 ± 0.0058 with a Z-score of −16.206 was inconsistent with our hypothesis of significant introgression from the pilot whale into type D, but not into the Crozet killer whale, with counts of ABBA sites = 14,976 (sites where type D and pilot whale shared a derived allele) and counts of BABA sites = 18,068 (sites where the Crozet killer whale and pilot whale shared a derived allele). We account for the excess sharing of derived alleles between the Crozet killer whale and pilot whale, relative to the number of derived alleles shared between type D and the pilot whale, as being due to derived alleles drifting to extinction in type D because of the lower long-term effective population size (see below). Whilst bottlenecks in the outgroup population and ancestral population of the ingroup have been shown not to influence the D-statistic (Durand et al. 2011), that study did not test the effect of bottlenecks in one of the ingroup populations.

Discussion

The first genome sequences of type D killer whales reveal the most extreme example of long-term inbreeding within this species to date (see for comparison Foote et al. 2021). The extent of the type D killer whale genome comprised by homozygous tracts exceeds that found in some of the most prominent examples of inbreeding in wild populations, e.g. eastern lowland (Gorilla beringei graueri) and mountain gorillas (G. b. beringei) (Xue et al. 2015) and Isle Royale wolves (Canis lupis) (Robinson et al. 2019). Our demographic reconstruction suggests that a long-term low effective population size, rather than a recent population collapse during the Anthropocene, underlies these very high inbreeding estimates. A consequence of this long-term low effective population size is that most variation coalesces in a relatively recent common ancestor. In other words, much of the genome represents a relatively shallow pedigree (<25,000 yr or ~1,000 generations). Ancient coalescent events reflecting the deeper pedigree are our window further back into the past (Wakeley 2009). The reduction of older coalescence within the genome limits our ability to infer the deeper demographic and evolutionary history of type D killer whales to the same extent that we have with other killer whale populations (e.g. Foote et al. 2019, 2021). Increasing the sample size (and thereby cross-coalescence among genomes) of unrelated type D genomes, and additional sequencing of the New Zealand stranded specimen, may provide greater resolution of the demographic history through time. A larger sample size would also provide a better understanding of the extent to which patterns described here reflect recent inbreeding in a few individuals, or long-term low effective population size within type D.

The timing of the inferred decline in N_e corresponds with a significant reduction of flow speed of the Antarctic Circumpolar Current (ACC) (Wu et al. 2021). Although type D has a circumpolar range, sightings suggest this is constrained to a narrow latitudinal sub-Antarctic range between 46°S and 62°S (Fig. 1). The impact of changes to the ACC flow rate during the last Glacial, and the narrow range width, may be driving factors that have impacted the demography of type D killer whales. However, it is important to caution that our inferences of demographic and evolutionary history of type D killer whales are obtained from patterns of coalescence within their genomes. Whilst these analyses give us some indication of the timing and scale of the changes in coalescent rates and effective population size, they do not explicitly inform us of the extrinsic factors responsible. Any interpretation of those causes is speculative.

We found no evidence of introgression into type D from the partially sympatric long-finned pilot whale, with which it shares some morphological traits. However, our characterization of the genomes of type D killer whales does provide insights into the potential genetic underpinning of the distinctive morphological traits by which type D was first defined (Pitman et al. 2011). Notably, the most morphologically distinct populations of killer whales, e.g. the Antarctic types, the Scottish West Coast population, and type D (Pitman and Ensor 2009; Pitman et al. 2011; Mäkeläinen et al. 2014), are also the populations with the lowest heterozygosity (Foote et al. 2021). Recessive mutations, which are expressed in homozygous genotypes, can contribute to the heritability of complex traits such as morphology (Lohmueller 2014). A family group from the North Pacific transient ecotype, known as the T2s, provide an illustrative example. The T2s included several individuals with varied head and jaw abnormalities, and at least 2 individuals with white pigmentation (Hoyt 1981; Ford et al. 1999). The white pigmentation in 1 individual was identified as a symptom of an autosomal recessive disorder, Chédiak-Higashi syndrome (Taylor and Farrell 1973). The extent of fixation of homozygous tracts in the genomes of type D killer whales through inbreeding, and the consequential expression of linked recessive alleles, are potential candidates for explaining the fixation of the distinct morphology of type D.

The long-term, low effective population size in type D is commensurate to drift, which is reflected in the low covariance of allele frequencies and IBS with other killer whale populations. Drift drives the random establishment of LD among loci. This may in turn play a role in strengthening endogenous barriers to gene flow in type D. Many studies have proposed that speciation is promoted through increasing LD between alleles at loci associated with endogenous reproductive isolation, and loci associated with exogenous selection pressures, e.g. local adaptation (Felsenstein 1974, 1981; Barton and Bengtsson 1986; Butlin 2005; Barton and de Cara 2009; Bierne et al. 2011). Both endogenous (those associated with intrinsic reproductive isolation) and exogenous alleles (those associated with local adaptation) will be targets of selection and will thus reduce gene flow of neutral alleles through the action of linked selection (Bierne et al. 2011). As LD increases among targets of selection, its action can become cumulative across the many loci (Turner 1967; Kruuk et al. 1999; Bierne et al. 2011; Flaxman et al. 2014; Oomen et al. 2020). Thus, increasing LD and selection can drive progress toward speciation, but this is in turn opposed by recombination, which will act to break up the association between alleles at 2 or more loci (Felsenstein 1981). The permeability to gene flow will therefore depend upon the local recombination rate (Barton 1979; Bierne et al. 2011). Based on our results, the genomes of type D killer whales have extensive tracts of homozygosity covering the majority of the autosomes, and the rate of recombination of different haplotypes is low relative to other killer whale genomes. In theory these characteristics should promote the reproductive isolation of type D from other killer whales.

Speciation could be further enhanced through assortative mating. Tixier et al. (2016) describe and quantify the complete social segregation of type D and the local Crozet Archipelago population of killer whales from over a decade of observations (Tixier et al. 2014). Type D and Crozet killer whales have been observed in sympatry when feeding on the same long-line fishing set (Tixier et al. 2016), so there is the opportunity for the 2 types to interact and interbreed. However, type D and Crozet killer whales actively segregate when they encounter one another (Tixier et al. 2016).

In summary, the first genome sequences of type D killer whales have proved to be as distinctive as the morphology of this type. Despite the challenges of working with a small sample size from a nonequilibrium population, our genomic analyses provide insight into how type D can retain this distinctiveness, despite co-occurrence with other killer whale populations. The killer whales found around the Antarctic continent: types B1, B2, and C, also have low genetic diversity. However, their genomes are characterized by an excess of short ROH (<1.5 Mb), relative to type D. A key question then is whether the genomes of the ancestors of types B1, B2, and C were once as homogenous as type D genomes, and whether over time, rare episodic gene flow and recombination will gradually break up the long haplotypes and erode LD in the genomes of type D killer whales. Alternatively, have type D killer whales progressed beyond the threshold where the barriers to gene flow with other killer whales can be broken down? Given our findings of extreme inbreeding and the absence of evidence of any outbreeding, the taxonomic status of type D killer whales deserves careful consideration.

Postscript: At the time of writing the final revision of this study and 67 yr after the last known type D stranding, a probable type D (based on photographs) has stranded in Magallanes, near Punta Delgada, Chile. We intend to rerun the analyses included in this paper incorporating this individual. Upon completion of these analyses the additional results will be deposited in the associated Dryad repository associated with this study (https://doi.org/10.5061/dryad.wpzgmsbr7).

Funding

ADF was funded by European Research Council (ERC) Consolidator Grant ERC-COG-101045346 “EXPLOAD.” Mridula Srinivasan and the Protected Species Science Branch of NOAA Fisheries Office of Science and Technology funded genome sequencing of DNA extracted from the type D biopsy samples. Sequencing of the Paraparaumu Beach sample was funded by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 663830 awarded to ADF.

Acknowledgments

We would like to thank Martin Stoffel and 1 anonymous reviewer for their comments which greatly improved the final manuscript, and Marty Kardos for providing valuable feedback on an earlier draft as part of the NOAA internal ­review process. Additionally, we are appreciative of the editorial team, Anjanette Baker, Bill Murphy, and Scott Baker, for accommodating the inclusion of the postscript right at the late stage of manuscript revisions. We thank Songhai Li and Yaolei Zhang generously provided their draft chromosomal reference assembly ahead of their own publication, and Ryan Waples for advice on running IBSrelate. Ngā mihi to Ngāti Toa Rangatira for allowing us to work with their taonga, and Colin Miskelly and Anton van Helden—curators at Museum of New Zealand—Te Papa Tongarewa for providing the historical sample. Ngā mihi to Te Rerekohu Tuterangiwhiu for translating the abstract to te reo Māori. Ben Wallis and the vessel Australis; the toothfish fishers and fishing companies of South Chile that helped with finding type Ds during the dedicated trip of 2019 (especially E. Infante from AOBAC). All other contributors to the opportunistic sightings of type Ds (fishery observers and crews onboard toothfish longliners off Prince Edward/Marion, Crozet, and Kerguelen islands with support from CAPFISH, SARPC, the MNHN, and the DPQM and Réserve Naturelle in the TAAF administration; the crew of the R/V Astrolabe, with support from the IPEV program no. 109 [Ornithoeco], and of the Bob Barker, etc.). Thanks also to Rio Seco Museum for collecting the sample of stranded animal referred to in the postscript. The CITES export permit #19CL000006WS was issued by the National Fisheries and Aquaculture Service of Chile. The research permit Res. 1811/2017 and its modification Res. 4402/2018 were issued by the Undersecretariat of Fisheries and Aquaculture of Chile. Remaining biopsy tissues are stored in the SWFSC Marine Mammal and Sea Turtle Research (MMASTR) collection.