RAD Sequencing and a Hybrid Antarctic Fur Seal Genome Assembly Reveal Rapidly Decaying Linkage Disequilibrium, Global Population Structure and Evidence for Inbreeding.

Recent advances in high throughput sequencing have transformed the study of wild organisms by facilitating the generation of high quality genome assemblies and dense genetic marker datasets. These resources have the potential to significantly advance our understanding of diverse phenomena at the level of species, populations and individuals, ranging from patterns of synteny through rates of linkage disequilibrium (LD) decay and population structure to individual inbreeding. Consequently, we used PacBio sequencing to refine an existing Antarctic fur seal (Arctocephalus gazella) genome assembly and genotyped 83 individuals from six populations using restriction site associated DNA (RAD) sequencing. The resulting hybrid genome comprised 6,169 scaffolds with an N50 of 6.21 Mb and provided clear evidence for the conservation of large chromosomal segments between the fur seal and dog (Canis lupus familiaris). Focusing on the most extensively sampled population of South Georgia, we found that LD decayed rapidly, reaching the background level by around 400 kb, consistent with other vertebrates but at odds with the notion that fur seals experienced a strong historical bottleneck. We also found evidence for population structuring, with four main Antarctic island groups being resolved. Finally, appreciable variance in individual inbreeding could be detected, reflecting the strong polygyny and site fidelity of the species. Overall, our study contributes important resources for future genomic studies of fur seals and other pinnipeds while also providing a clear example of how high throughput sequencing can generate diverse biological insights at multiple levels of organization.

Advances in short read sequencing technologies, in particular Illumina sequencing, have made it possible to generate genome assemblies as well as dense genetic marker datasets for practically any organism (Ekblom and Galindo 2011; Ellegren 2014). However, assemblies based solely on short read data tend to be highly fragmented, even with assembly strategies that incorporate medium length insert libraries (Gnerre ). Consequently, although such assemblies can be generated rapidly and cheaply, there has been growing interest in technologies that incorporate longer range information to improve scaffold length and contiguity. For example, Pacific Biosciences (PacBio) single molecule real-time (SMRT) sequencing generates read lengths in the order of several kilobases (kb) that have proven effective in gap filling, resolving complex repeats and increasing contig lengths across diverse taxa (English ; Conte and Kocher 2015; Pootakham ).

In parallel to these and related developments in genome sequencing technologies, reduced representation sequencing approaches such as restriction site associated DNA (RAD) sequencing (Baird ; Peterson ) are providing unprecedented levels of genetic resolution for population genetic and genomic studies (Morin ; Stapley ; Seeb ). By sequencing and assembling short stretches of DNA adjacent to restriction cut sites and interrogating the resulting tags for sequence polymorphisms, RAD sequencing can facilitate the acquisition of large genome-wide distributed single nucleotide polymorphism (SNP) datasets incorporating multiple individuals.

The above approaches show great promise for studying wild populations where genomic resources are typically absent. For example, information from model organisms with well-characterized genomes can facilitate studies of their wild relatives as long as patterns of synteny between the two can be established. Knowledge of synteny can facilitate gene annotation, assist in gene mapping and help to elucidate the genetic basis of fitness variation by identifying genes closely linked to loci responsible for inbreeding depression (Johnston ; Ekblom and Wolf 2014; Kardos ).

High density SNP markers mapped to a reference genome can furthermore provide insights into processes that shape levels of variation within genomes. For example, the positional information of genomic loci can be used to characterize patterns of linkage disequilibrium (LD). LD is a central concept in population genetics because it is closely associated with factors such as effective population size (N), genetic drift, historical fluctuations in population size, population structure, inbreeding and recombination (Slatkin 2008). Understanding the strength and extent of LD can aid in the inference of demographic history and has important implications for identifying genetic variants underlying key fitness traits through genome-wide association analyses or quantitative trait locus mapping (Carlson ; Miller ; Kardos et al. 2016). Nevertheless, the genomic pattern of LD has only been described in a handful of wild organisms. Typically, LD decays within a few tens to hundreds of kilobases (kb) in large and unstructured populations (Poelstra ; Kawakami ; Vijay ), but can extend for several megabases (Mb) in smaller, isolated, heavily bottlenecked and/or inbred populations of species such as wolves and sheep (Hagenblad ; Miller ).

In addition to facilitating the characterization of genome-wide patterns of variation, dense genomic markers can also be used to describe variation at the population and individual level, even without positional information. For example, studies are increasingly employing approaches such as RAD sequencing to obtain large datasets in order to reliably characterize population genetic structure (Malenfant ; Benestan ; Younger ) and many are uncovering patterns that had previously gone undetected (Reitzel ; Ogden ; Vendrami ). A precise understanding of population structure is critical for the delineation of management units for conservation (Bowen ) as well as for avoiding false positives in genome-wide association studies (Johnston ). Knowledge of population structure can also help to elucidate contemporary and historical barriers to gene flow (McRae ; Hendricks ) as well as patterns of extinction and recolonization (McCauley 1991).

A major topic of interest at the level of the individual is the extent to which inbreeding occurs in natural populations (Kardos et al. 2016) and its consequences for fitness variation and population demography (Keller and Waller 2002). Pedigree-based studies, typically of isolated island populations and often involving polygynous species, have uncovered widespread evidence of inbreeding in the wild (Marshall ; Townsend and Jamieson 2013; Nietlisbach ). However, the extent of inbreeding in large, continuous and free-ranging populations remains open to question. On the one hand, simulations have suggested that inbreeding will be absent from the vast majority of wild populations with the possible exception of highly polygynous and/or structured populations (Balloux ). On the other hand, associations between microsatellite heterozygosity and fitness (heterozygosity fitness correlations, HFCs) have been described in hundreds of species (Chapman ) and it has been argued that these are highly unlikely to arise in the absence of inbreeding (Szulkin ). Due to the high sampling variance of microsatellites, there has been growing interest in the use of high density SNP data to reliably quantify inbreeding, and recent empirical and simulation studies suggest that this can be achieved with as few as 10,000 SNPs (Kardos ; 2018). Consequently, with approaches like RAD sequencing, it should be possible to quantify the variation in inbreeding in arguably more representative wild populations.

The Antarctic fur seal (Arctocephalus gazella) is an important marine top predator that has been extensively studied for several decades, yet many fundamental aspects of its biology remain poorly understood. This highly sexually dimorphic pinniped has a circumpolar distribution and breeds on islands across the sub-Antarctic, with 97% of the population concentrated on South Georgia in the South Atlantic (Figure 1). The species was heavily exploited by 18th and 19th Century sealers and was thought to have gone extinct at virtually all of its contemporary breeding sites (Weddell 1825). However, in the 1930s a small breeding population was found at South Georgia (Bonner 1964; Payne 1977), which in the following decades increased to number several million individuals (Boyd 1993). While it is generally believed that the species former range was recolonized by emigrants from this large and rapidly expanding population (Boyd 1993; Hucke-Gaete ), one would expect to find little or no population structure under such a scenario. However, a global study using mitochondrial DNA resolved two main island groups (Wynen ) while microsatellites uncovered significant differences between South Georgia and the nearby South Shetland Islands (Bonin ), implying that at least two relict populations must have survived sealing.

Figure 1

Individual assignment to genetic clusters based on STRUCTURE analysis for K = 4 using 27,592 SNPs. Each horizontal bar represents a different individual and the relative proportions of the different colors indicate the probabilities of belonging to each cluster. Individuals are separated by sampling locations as indicated on the map.

Antarctic fur seals have been intensively studied for several decades at a small breeding colony on Bird Island, South Georgia, where a scaffold walkway provides access to the animals for the collection of detailed life history and genetic data. Genetic studies have confirmed behavioral observations of strong polygyny (Bonner 1964) by showing that a handful of top males father the majority of offspring (Hoffman ). Furthermore, females exhibit strong natal site fidelity, returning to within a body length of where they were born to breed (Hoffman and Forcada 2012), while adults of both sexes are highly faithful to previously held breeding locations (Hoffman ). Together these behavioral traits may increase the risk of incestuous matings. In line with this, heterozygosity measured at nine microsatellites has been found to correlate with multiple fitness traits including early survival, body size and reproductive success (Hoffman ; 2010; Forcada and Hoffman 2014). However, such a small panel of microsatellites cannot provide a very precise estimate of inbreeding (Slate ; Balloux ; Miller ) and therefore high density SNP data are required to generate more detailed insights into the variance in inbreeding in the population.

Here, we used PacBio sequencing to improve an existing Antarctic fur seal genome assembly comprising 8,126 scaffolds with an N50 of 3.1 Mb (Humble ). We additionally RAD sequenced 83 individuals, mainly from South Georgia but also from an additional five populations, to generate a large dataset of mapped genetic markers. The resulting data were then used to investigate synteny with the dog (Canis lupus familaris) and, within the focal South Georgia population, to characterize the pattern of LD decay as well as the variance in inbreeding. Finally, using data from both RAD sequencing and 27 microsatellites, we investigated the strength and pattern of population structure across the species range and compared the ability of the two marker types to resolve genetic differences between island groups. Our hypotheses were as follows: (i) We expected to find strong synteny between the fur seal and dog, the closest relative with an annotated, chromosome-level genome assembly; (ii) Given that Antarctic fur seals are free-ranging with large population sizes, we predicted that LD would decay very rapidly, although it is also possible that the historical bottleneck could have resulted in elevated levels of LD; (iii) We hypothesized that nuclear markers would detect the same two island groups as previously found with mitochondrial DNA as well as possibly resolve finer scale structuring. We furthermore predicted that RAD sequencing would provide greater power to capture genetic differences than microsatellites; (iv) Finally, we anticipated finding variation in inbreeding consistent with knowledge of the species mating system as well as previous studies documenting HFCs.

Materials and Methods

Hybrid genome assembly and PacBio DNA library preparation

We first used the program GapCloser v1.12 to fill gaps in the existing fur seal genome v1.02 (Humble ) (NCBI SRA: BioProject PRJNA298406) based on the paired end information of the original Illumina reads. This approach closed 45,852 gaps and reduced the amount of N space in the assembly from 115,235,953 bp to 78,393,057 bp (v1.1, Table 1). Following this, we generated SMRT sequencing data from the DNA used for the original genome assembly (NCBI SRA: BioSample SAMN04159679) following the protocol described in Pendleton . First, 10 μg of pure genomic DNA was fragmented to 20 kb using the Hydroshear DNA shearing device (Digilab, Marlborough, MA) and size-selected to 9–50 kb using a Blue Pippin according to the standard Pacific Biosciences SMRT bell construction protocol. The library was then sequenced on 64 PacBio RSII SMRT cells using the P6–C4 chemistry. This yielded a total of 58 Gb (∼19x) of sequencing data contained within 8,101,335 subread bases with a mean read length of 7,177 bp (median = 6,705 bp; range = 50–54,622 bp).

Table 1

Genome assembly statistics for successive improvements to the original Antarctic fur seal genome assembly

	v1.0.2 ALLPATHS3	v1.1 GapCloser	v1.2 PBJelly2	v1.3 Quiver	v1.4 Pilon
Number of scaffolds	8,126	8,126	6,170	6,170	6,169†
N90a	890,836 (768)	888,912 (768)	1,624,547 (387)	1,511,352 (387)	1,542,705 (387)
N50a	3,169,165 (233)	3,165,747 (233)	6,454,664 (108)	6,076,522 (108)	6,207,322 (108)
N10a	8,459,351 (25)	8,458,289 (25)	17,733,103 (11)	16,529,571 (11)	16,861,656 (11)
Longest scaffold (bp)	13,012,173	12,999,316	34,690,325	32,399,786	33,062,611
Total size (bp)	2,405,038,055	2,403,626,805	2,426,014,533	2,268,217,244	2,313,485,084
Gaps present (%)	4.79	3.26	0.62	0.57	0.55
Number of gaps	136,284	90,432	45,102	22,783	20,611
Average gap size (bp)	845.56	866.87	331.16	570.37	613.39

Size in bp (number of scaffolds)

Excluding the mitochondrial genome, which was filtered out by Pilon.

Next, we used PBJelly v15.8.24 and blasr (https://github.com/PacificBiosciences/blasr) with default parameters to align the PacBio sequencing reads to the gap-closed assembly to generate a hybrid genome assembly (v1.2). Lastly, we followed a two-step strategy to remove any indels introduced by single molecule real-time sequencing (Ross ). We first used Quiver (contained in the SMRT/2.3.0 suite: GenomicConsensus v0.9.2) with the refineDinucleotideRepeats option to perform initial assembly error correction. Due to this step being computationally demanding, we ran it separately for each scaffold. Next, we mapped the original Illumina reads (Humble ) to the quiver assembly (v1.3) using BWA MEM v0.7.15 (Li 2013) and used Picard tools to sort and mark duplicates. We then used PILON v1.22 (Walker ) to perform the final error correction step to generate assembly v1.4.

Genome alignment

We aligned the fur seal scaffolds from assembly v1.4 to the dog genome (Canis lupus familiaris assembly version CanFam3.1, GenBank accession number GCA_000002285.2) using LAST v746 (Kiełbasa ). First, the dog genome was prepared for alignment using the command lastdb. We then used lastal and last-split in combination with parallel-fastq to align the fur seal scaffolds against the dog genome. Using the program MafFilter v1.0.0, we then processed the resulting multiple alignment format (maf) file and estimated pairwise sequence divergence between the two species (Dutheil ). Finally, we extracted alignment coordinates from the maf file using bash commands to allow subsequent visualization with the R package RCircos v1.2.0 (Zhang ).

Sampling and DNA extraction

Tissue samples were collected from 57 Antarctic fur seal individuals from Bird Island, South Georgia. These comprised 24 triads consisting of 24 pups, 16 mothers and 17 fathers of which five mothers and two fathers sired multiple offspring (Figure S1). Additional samples were obtained from the main breeding colonies across the species range (Figure 1): Cape Shirreff in the South Shetlands (n = 6), Bouvetøya (n = 5), Kerguelen (n = 5), Heard Island (n = 5) and Macquarie Island (n = 5). Skin samples were collected from the inter-digital margin of the fore-flipper using piglet ear notching pliers and stored in 20% dimethyl sulphoxide saturated with NaCl at –20°. Skin samples from the South Shetlands were collected using a sterile 2mm biopsy punch and stored in 95% ethanol. Total genomic DNA was extracted using a standard phenol-chloroform protocol (Sambrook ).

Microsatellite genotyping

All samples were genotyped at 27 polymorphic microsatellite loci (see Table S1), previously found to be in Hardy-Weinberg equilibrium (HWE) in the study population at South Georgia and unlinked (Stoffel ; Peters ). These loci were PCR amplified in three separate multiplex reactions (see Table S1) using a Type It Kit (Qiagen). The following PCR profile was used for all multiplex reactions except for multiplex one: initial denaturation of 5 min at 94°; 28 cycles of 30 sec at 94°, 90 sec at 60°, and 30 sec at 72°, followed by a final extension of 30 min at 60°. The PCR profile of multiplex one only differed from this protocol in the annealing temperature used, which was 53°. Fluorescently labeled PCR products were then resolved by electrophoresis on an ABI 3730xl capillary sequencer and allele sizes were scored using GeneMarker v1.95. To ensure high genotype quality, all traces were manually inspected and any incorrect calls were adjusted accordingly.

RAD library preparation and sequencing

RAD libraries were prepared using a modified protocol from Etter with minor modifications as described in Hoffman . Briefly, 400 ng of genomic DNA from each individual was separately digested with SbfI followed by the ligation of P1 adaptors with a unique 6 bp barcode for each individual in a RAD library, allowing the pooling of 16 individuals per library. Libraries were sheared with a Covaris S220 and agarose gel size-selected to 300–700 bp. Following 15–17 cycles of PCR amplification, libraries were further pooled using eight different i5 indices prior to 250 bp paired-end sequencing on two Illumina HiSeq 1500 lanes.

SNP genotyping

Read quality was assessed using FastQC v0.112 and sequences were trimmed to 225 bp and demultiplexed using process_radtags in STACKS v1.41 (Catchen ). We then followed GATK’s best practices workflow for variant discovery (Poplin ). Briefly, individual reads were mapped to the Antarctic fur seal reference genome v1.4 using BWA MEM v0.7.10 (Li 2013) with the default parameters. Any unmapped reads were removed from the SAM alignment files usingSAMtools v1.1 (Li 2011). We then used Picard Tools to sort each SAM file, add read groups and remove PCR duplicates. Prior to SNP calling, we performed indel realignment to minimize the number of mismatching bases using the RealignerTargetCreator and IndelRealigner functions in GATK v3.6. Finally, HaplotypeCaller was used to call variants separately for each individual. Genomic VCF files were then used as input to GenotypeGVCFs for joint genotyping. The resulting SNP dataset was then filtered to include only biallelic SNPs using BCFtools v1.2 (Li 2011) to obtain a dataset of 677,607 SNPs genotyped in 83 individuals. Subsequently, we applied a variety of filtering steps according to the analysis being performed as shown in Figure S2 and described below.

SNP validation

To provide an indication of the quality of our SNP dataset, we attempted to validate a representative subset of loci using Sanger sequencing. First, we randomly selected 50 loci whose 70 bp flanking sequence contained no secondary SNPs and mapped uniquely to the fur seal reference genome and with initial depth of coverage and minor allele frequency (MAF) filters of 5 and 0.05 respectively. We then designed oligonucleotide primers using Primer 3 (Untergasser ) to PCR amplify each putative SNP together with 100–200 bp of flanking sequence. Each locus was PCR amplified in one fur seal individual that had been genotyped as homozygous at that locus and one that had been genotyped as heterozygous. PCRs were carried out using 1.5 µL of template DNA, 20 mM Tris–HCl (pH 8.3), 100 mM KCl, 2 mM MgCl2, 10x Reaction Buffer Y (Peqlab), 0.25 mM dNTPs, 0.25 mol/L of each primer and 0.5U of Taq DNA polymerase (VWR). The following PCR profile was used: one cycle of 1.5 min at between 59° and 62° depending on the primers used (Table S2), 60 sec at 72°; and one final cycle of 7 min at 72°. 5 µL of the resulting PCR product was then purified using shrimp alkaline phosphatase and exonuclease I (NEB) following the manufacturer’s recommended protocol. All fragments were then sequenced in both directions using the Applied Biosystems BigDye Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Scientific) and analyzed on an ABI 3730xl capillary sequencer. Forward and reverse reads were aligned using Geneious v10.2.3 (Kearse ). Heterozygous sites were identified as those with two peaks of roughly equal intensity but with around half the intensity of a homozygote.

Linkage disequilibrium decay

Prior to quantifying linkage disequilibrium, we filtered the SNP dataset as shown in Figure S2A. First, to minimize the occurrence of unreliable genotypes, we removed individual genotypes with a depth of coverage below eight or above 30 using VCFtools (Danecek ). Genotypes with very low depth of coverage have a greater likelihood of being called incorrectly as it can be difficult to distinguish between homozygotes and heterozygotes when very few reads are present. Similarly, genotypes with very high depth of coverage are more likely to be spurious as high coverage can result from misalignment due to the presence of paralogous loci or repeats (Fountain ). Second, as including SNPs from short scaffolds can downwardly bias LD values, we retained only SNPs located on the longest 100 scaffolds of the assembly (min length = 6.6 Mb, max length = 33.1 Mb). Third, as an additional quality filtering step, we used information on known parental relationships to identify loci with Mendelian incompatibilities using the mendel function in PLINK v1.9 and removed these from the dataset. Fourth, to minimize any possible confounding effects of population structure and relatedness, we focused on adult individuals sampled from the single largest population of South Georgia. Fifth, to provide an informative dataset while further minimizing genotyping error, we discarded SNPs with a MAF of less than 0.1 and/or called in less than 50% of individuals using PLINK. As a final quality control step, we also removed SNPs that did not conform to Hardy-Weinberg equilibrium (HWE) with a p-value threshold < 0.001 using the–hwe function in PLINK.

Using the final dataset of 25,068 SNPs genotyped in 33 individuals (Figure S2A), we used the –r2 function in PLINK to quantify pairwise LD between all pairs of SNPs located within 500 kb of each other. We visualized LD decay with distance by fitting a nonlinear regression curve using the nls package in R, where the expected value of r2 under drift-recombination equilibrium () was expressed according to the Hill and Weir function (Hill and Weir 1988):where n is the number of gametes scored, P = 4N, N is the effective population size and c is the recombination fraction between sites (Remington ).

Population structure

Prior to quantifying population structure, we filtered the full SNP dataset as shown in Figure S2B. We did not initially filter the dataset for SNPs with low depth of coverage, as for the analysis of population structure we wanted to retain as many SNPs as possible that were genotyped across all of the populations. We also did not remove individuals with large amounts of missing data in order to maximize the representation of all populations in the final dataset. Nevertheless, because closely related individuals can bias population genetic structure analysis by introducing both Hardy-Weinberg and linkage disequilibrium (Rodríguez-Ramilo and Wang 2012; Wang 2017), we used known parentage information to remove adults and related pups (full and half siblings) from the South Georgia dataset. This strategy also enabled us to attain a more even number of samples per population (Puechmaille 2016). Second, SNPs with a MAF of less than 0.05 and/or called in less than 99% of individuals were discarded using VCFtools. Third, we removed SNPs that were out of HWE at the population level with a p-value threshold < 0.001 using the –hwe function in PLINK. Fourth, SNPs were pruned for LD using the –indep function in PLINK. We used a sliding window of 50 SNPs, a step size of 5 SNPs and removed all variants in a window above a variance inflation factor threshold of 2, corresponding to r2 = 0.5.

Using the final dataset of 27,592 SNPs genotyped in 37 individuals (Figure S2B), we first visualized population structure by performing a principal components analysis (PCA) using the R package adegenet (Jombart 2008). We then used a Bayesian clustering algorithm implemented by the program STRUCTURE to identify the number of genetic clusters (K) present in the dataset. We performed STRUCTURE runs for values of K ranging from one to six, with five simulations for each K and a burn-in of 100,000 iterations followed by 1,000,000 Markov chain Monte Carlo iterations. We used the admixture and correlated allele frequency models without sampling location information. The R package pophelper (Francis 2017) was then used to analyze the STRUCTURE results, parse the output to CLUMPP for averaging across iterations and for visualizing individual assignment probabilities. The optimal K was selected based on the maximum value of the mean estimated ln probability of the data (Ln Pr(X | K)) as proposed by Pritchard and the ΔK method of Evanno . For comparison, we also implemented the above analyses using microsatellite data for the same individuals.

Inbreeding coefficients

Prior to quantifying inbreeding, we filtered the SNP dataset as shown in Figure S2C. First, we wanted a dataset with as few gaps as possible so we discarded one individual with more than 90% missing data. Second, we removed individual genotypes with a depth of coverage below eight or above 30 using vcftools. Third, we removed loci with Mendelian incompatibilities, and fourth, we again restricted the dataset to the focal population of South Georgia. Fifth, we discarded SNPs with a MAF of less than 0.05 and/or called in less than 70% of individuals using vcftools. Finally, we filtered the SNPs for HWE as described previously and pruned linked SNPs out of the dataset using the –indep function in PLINK with the parameters shown above.

Using the final dataset of 9,853 SNPs genotyped in 56 individuals (Figure S2C), we calculated four genomic estimates of individual inbreeding: standardized multi-locus heterozygosity (sMLH; Coltman ), an estimate based on the variance of additive genotype values (), an estimate based on excess homozygosity () and an estimate based on the correlation of uniting gametes, which gives more weight to homozygotes of the rare allele at each locus () (Yang ). sMLH was calculated using the sMLH function in the R package inbreedR (Stoffel ) and and were calculated in GCTA v1.24.3 (Yang ). When inbreeding varies among individuals, heterozygosity measured at one part of the genome will be correlated with heterozygosity measured at another, a phenomenon termed identity disequilibrium (ID) (Weir and Cockerham 1973; Szulkin ). We quantified the extent of ID using the measure in the R package inbreedR (Stoffel ) where significant values provide support for variance in inbreeding in the population. Finally, we compared the resulting value with the variance in our inbreeding coefficients to determine the expected correlation between estimated ( and realized levels of inbreeding (Szulkin ) given as:

Data availability

The PacBio sequencing reads and the final Antarctic fur seal genome assembly are available at the European Nucleotide Archive (ENA) under BioProject ID PRJEB26995. The raw RAD sequencing reads are available at the short-read archive (SRA) of the National Biotechonology Centre of Information (NCBI) under BioProject ID PRJNA473050, SRA accession SRP148937. Code for the analysis is available at https://github.com/elhumble/seal_rad_2018. Supplemental material and SNP genotypes are available at Figshare: https://doi.org/10.25387/g3.6713975.

Results

Hybrid genome assembly

We used PacBio sequencing to improve an existing Antarctic fur seal genome assembly. Using PBJelly, we were able to close a total of 45,394 gaps, resulting in a 40% reduction in overall gap space (assembly v1.2, Table 1). Subsequent assembly correction with Quiver resulted in a total of 11,319,546 modifications to the PBJelly assembly consisting of 291,179 insertions, 1,117,226 substitutions and 9,911,141 deletions. Finally, PILON corrected 653,246 homozygous insertions (885,794 bp), 87,818 deletions (127,024 bp) and 34,438 single-base substitutions and closed an additional 2,172 gaps in the Quiver assembly. Overall, gap closing and error correction resulted in a hybrid Antarctic fur seal assembly with a total length of 2.3 Gb (v1.4, Table 1). The number of scaffolds in the genome was reduced from 8,126 to 6,169 such that 50% of the final assembly is now contained within the longest 108 scaffolds (Table 1).

Genome synteny

To investigate synteny between the Antarctic fur seal and the dog, we aligned the fur seal scaffolds to the dog genome (CanFam3.1). We estimated overall sequence divergence between the two species to be 13.8%. Visualization of the full alignment revealed that all of the dog chromosomes are represented in the fur seal assembly (Figure S3). Alignment of the 40 longest fur seal scaffolds (min length = 10.7 Mb, max length = 33.1 Mb) revealed strong chromosomal synteny between the two genomes, with the vast majority of the fur seal scaffolds mapping exclusively or mainly to a given dog chromosome (Figure 2). Specifically, for 37 of the scaffolds, over 90% of the total alignment length was to a single dog chromosome, with 26 of those aligning exclusively to a single dog chromosome. Only one scaffold (S4 in Figure 2) aligned in roughly equal portions to two different dog chromosomes (62% to D5 and 38% to D26).

Figure 2

Synteny of the longest 40 Antarctic fur seal scaffolds (10.7–33.1 Mb; right, prefixed S) with dog chromosomes (left, prefixed D). Mapping each fur seal scaffold to the dog genome resulted in multiple alignment blocks (mean = 2.1 kb, range = 0.1–52.8 kb) and alignments over 5 kb are shown.

RAD sequencing and SNP discovery

RAD sequencing of 83 fur seal individuals generated an average of 5,689,065 250 bp paired-end reads per individual. After mapping these reads to the reference genome, a total of 677,607 biallelic SNPs were discovered using GATK’s best practices workflow for variant discovery (see Materials and methods for details), with an average SNP coverage of 727. We then filtered the dataset in three different ways (Figure S2) to generate datasets suitable for the analysis of LD decay, population structure and inbreeding respectively.

To provide an indication of the quality of our SNP dataset, we used Sanger sequencing to validate 50 randomly selected loci. For each locus, we sequenced a single heterozygote and a single homozygote individual based on the corresponding GATK genotypes. For 40 of these loci, we successfully obtained genotypes for both individuals (Table S2). Concordance between the GATK and Sanger genotypes was high, with 76/80 genotypes being called identically using both methods, equivalent to a validation rate of 95%. The four discordant genotypes were all initially called as homozygous with GATK but subsequently validated as heterozygous with Sanger sequencing.

LD decay

The pattern of LD decay within South Georgia was quantified based on 25,068 SNPs genotyped in 33 individuals and located on the 100 longest fur seal scaffolds. LD was found to decay rather rapidly, with r2 reaching the background level by around 400 kb (Figure 3). Strong LD (r2 >= 0.5) decayed by 5 kb and moderate LD (r2 >= 0.2) by around 15 kb.

Figure 3

Plot of linkage disequilibrium (r2) against physical distance between SNPs in the Antarctic fur seal. LD was calculated using 25,068 filtered SNPs from the 100 largest scaffolds of 33 South Georgia adults. Gray points indicate observed pairwise LD values. The dark gray curve shows the expected decay of LD in the data estimated by nonlinear regression.

Next, we used a dataset of 37 pups genotyped at 27 microsatellites and 27,592 SNPs to quantify the pattern and strength of population structure across the species circumpolar range. PCA of the microsatellite dataset uncovered weak clustering with South Georgia, the South Shetlands and Bouvetøya tending to separate apart from Kerguelen, Heard and Macquarie Islands along the first PC axis (Figure 4A). However, considerable scatter and no clear pattern of separation was found along either PC2 or PC3 (Figures 4A and 4C). By contrast, population structure was more clearly defined in the PCA of the SNP dataset. Specifically, the first PC axis clearly resolved two distinct island groups, the first comprising South Georgia, the South Shetlands and Bouvetøya and the second comprising Kerguelen, Heard Island and Macquarie Island (Figure 4B). Within the first island group, Bouvetøya clustered apart from South Georgia and the South Shetlands along PC2 (Figure 4B) while all three locations clustered apart from one another along PC3 (Figure 4D).

Figure 4

Scatterplots showing individual variation in principal components (PCs) one and two (panels A and B), and one and three (panels C and D) derived from a principal component analysis conducted using 27 microsatellites (panels A and C) and 27,592 SNPs (panels B and D). The amount of variance explained by each PC is shown in parentheses.

To test whether population structure could be detected without prior knowledge of the sampling locations of individuals, we used a Bayesian approach implemented within STRUCTURE (Pritchard ). This program works by partitioning the data set in such a way that departures from Hardy-Weinberg and linkage equilibrium within the resulting groups are minimized. Separately for the microsatellite and SNP datasets, five replicate runs were conducted for each possible number of clusters (K) ranging from one, implying no population differentiation, through to six, which would imply that all of the populations are genetically distinct. For the microsatellite dataset, Ln Pr(X | K) and ΔK both peaked at two, indicating support for the presence of two genetically distinct populations (Figure S4A and C). Membership coefficients for the inferred clusters are summarized in Figure S5A and indicate the presence of a Western population comprising individuals from South Georgia, the South Shetlands and Bouvetøya, and an Eastern population comprising individuals from Kerguelen, Heard Island and Macquarie Island.

For the RAD dataset, Ln Pr(X | K) also peaked at two but remained high for K = 3 and 4, while ΔK reached its maximum at K = 4 (Figure S4B and D). To explore this further, we plotted membership coefficients for K = 2–6 for both the microsatellite and SNP datasets. For the former, no evidence of population structure was found beyond K = 2, with successive increases in K merely introducing additional admixture (Figure S5A). By contrast, for the latter, plots corresponding to K values greater than two clearly resolved further hierarchical structure (Figure S5B). Results for K = 4 are shown in Figure 1, in which Kerguelen, Heard and Macquarie Islands are resolved as a single population, while South Georgia, the South Shetlands and Bouvetøya can be readily distinguished based on their corresponding group membership coefficients.

Inbreeding

Finally, we investigated inbreeding in the focal population at South Georgia using data from 9,853 SNPs genotyped in 56 individuals (Figure 5A). Identity disequilibrium differed significantly from zero (0.0052; bootstrap 95% confidence interval = 0.0008–0.0091, P = 0.023, Figure 5B) providing evidence for variance in inbreeding within the sample of individuals. Each individual’s level of inbreeding was quantified from the SNP dataset using four different genomic inbreeding coefficients (sMLH, , and ). All four of these measures were inter-correlated, with correlation coefficients (r) ranging from 0.69 to 0.84. (Figure 5C–E). Furthermore, the variances of all four measures fell within the 95% confidence interval of (Figure 5B), suggesting that the expected correlation between the estimated and realized level of inbreeding does not differ significantly from one.

Figure 5

(A) Distribution of genomic inbreeding coefficients () for 56 individual fur seals from South Georgia; (B) Distribution of identity disequilibrium () estimates from bootstrapping over individuals. The vertical dashed line represents the empirical estimate and the horizontal black line shows the corresponding 95% confidence interval based on 1000 bootstrap replicates. The vertical colored lines represent the variance in four different inbreeding coefficients. Panels C, D and E show pairwise Pearson’s correlation coefficients between and , and sMLH.

Discussion

Advances in high throughput sequencing technology have afforded researchers the opportunity to generate genome assemblies and genomic marker datasets for virtually any species from which high quality DNA can be collected. These resources allow a broad variety of questions in ecology and evolution to be addressed with greater power and precision than was possible with traditional methods. In this study, we utilized PacBio sequencing to refine an existing Antarctic fur seal genome assembly and combined this with RAD sequencing to characterize synteny with the dog genome, elucidate the rate of LD decay, resolve global population structure and quantify the variance in inbreeding. Our results provide new insights at multiple levels of organization that enrich our understanding of an important Antarctic marine top predator and indicate the general promise of these and related approaches for tackling broad-reaching questions in population and evolutionary genetics.

An important outcome of this study is a significantly improved Antarctic fur seal genome assembly. This was achieved through three iterative steps involving gap filling, inclusion of long PacBio reads and assembly error correction respectively. Overall, the number of scaffolds was reduced by around one quarter, while N50 almost doubled to over 6 Mb and the proportion of gaps was reduced by around an order of magnitude to around half a percent. This represents an improvement over existing pinniped assemblies such as the walrus (Odobenus rosmarus divergens, GenBank accession number GCA_000321225.1) and Weddell seal (Leptonychotes weddellii, GenBank accession number GCA_000349705.1), which both have lower N50 values (2.6 and 0.9 Mb respectively). The improved Antarctic fur seal genome will therefore serve as an important resource for the wider pinniped community. However, there is still considerable room for improvement as a handful of other marine mammal genome assemblies incorporating longer range information show higher levels of contiguity, such as the polar bear (Ursus maritimus, N50 = 15.9 Mb), killer whale (Orcinus orca, N50 = 12.7 Mb) and Hawaiian monk seal (Neomonachus schauinslandi, N50 = 22.2 Mb) (Liu ; Foote ; Mohr ).

To further quantify genome quality and to explore patterns of synteny, we mapped the scaffolds of our new assembly to the dog genome. The resulting alignment revealed almost complete coverage of the dog chromosomes. This is in line with the observation that the total length of the assembly has not changed appreciably between versions and suggests that the assembly is near-complete, with the exception of the Y-chromosome for which sequence data are currently lacking as the genome individual is a female. In general, carnivore genomes show high levels of synteny (Arnason 1974; Ferguson-Smith and Trifonov 2007), with pinnipeds in particular exhibiting highly conserved karyotypes indicative of slow rates of chromosomal evolution (Beklemisheva ). By contrast, the domestic dog has an extensively re-arranged karyotype differentiated from the ancestral carnivore karyotype by over 40 separate fission events (Nie ). To provide insights into the extent of conservation of chromosomal blocks between seals and dogs, we mapped the longest 40 fur seal scaffolds to the dog genome. We found a clear pattern whereby all but one of the scaffolds mapped exclusively or mainly to single chromosomes, indicating the conservation of large genomic tracts often several Mb in length. The remaining scaffold mapped to two dog chromosomes in roughly equal proportions, suggestive of either a fission event in the lineage leading to dogs or a fusion event in the lineage leading to seals. By focusing only on the largest scaffolds, we had little power to detect multiple chromosomal rearrangements, although these are to be expected given a substantial increase in the number of chromosomes in dogs (2n = 74) relative to the seal (2n = 36) (Gustavsson 1964; Arnason 1974). Nevertheless, the observed high degree of synteny is consistent with previous studies revealing both strong sequence homology and the conservation of polymorphic loci between seals and dogs (Osborne ; Hoffman ).

SNP discovery and validation

Our study identified a total of 677,607 SNPs in a discovery pool of 83 individuals. These markers will be useful for future studies including the planned development of a high-density SNP array. However, not all SNPs are suitable for every analysis due to differential sensitivity to missing data, low depth of sequencing coverage and the inclusion of low frequency alleles (Shafer ). We therefore carefully considered how best to filter our SNP dataset for each of our main analyses. For LD decay, we applied relatively strict filters as we sought a high-quality dataset with consistently high coverage across individuals. For population structure, it was important to have as many SNPs as possible represented in all of the sampling locations, so we did not remove genotypes with low coverage but instead filtered the data to retain SNPs genotyped in at least 99% of individuals. Conversely, for the estimation of inbreeding, we honed in on a reduced subset of higher quality SNPs with greater average depth of coverage, that were in Hardy-Weinberg and linkage equilibrium, and which showed no evidence of Mendelian incompatibilities.

Even with stringent filtering, it is possible to retain SNPs in a dataset that have been called incorrectly. We therefore attempted to validate 50 randomly selected loci by Sanger sequencing selected individuals with homozygous and heterozygous genotypes as determined from the RAD data. For the 40 loci that we were able to successfully sequence, around 95% of the Sanger genotypes were identical to the RAD genotypes. Although this validation step required additional experimental effort, our results compare favorably with other studies (Cruz ; Bourgeois ) and thus give us confidence in the overall quality of our data.

We used the genomic positions of SNPs mapping to the largest 100 scaffolds to quantify the pattern of LD decay in the focal population of South Georgia. We found that LD decays rapidly, with moderate LD extending less than 20 kb. This is despite the species having experienced a population bottleneck in the 19th century, which would be expected to increase LD. A direct comparison with other organisms is hindered both by a paucity of data for most species and by the use of different measures for quantifying LD. However, our results are broadly in line with other wild vertebrate populations such as Alaskan gray wolves, polar bears and flycatchers, where moderate LD also extends less than 10–30 kb (Gray ; Malenfant ; Kardos et al. 2016). Extended LD has been documented in a number of species but in most cases this is associated with extreme bottlenecks, such as those experienced during domestication (Harmegnies ; McKay ; Meadows ). Although Antarctic fur seals are generally believed to have also experienced a strong historical bottleneck, a recent Bayesian analysis suggested that this may have been less severe than thought, with the effective population size probably falling to several hundred (Hoffman ). Furthermore, the population recovered from the bottleneck within a few generations, which could have mitigated the increased genetic drift and inbreeding effects that elevate and maintain strong LD. Additionally, the population is currently estimated to number around 2–3 million individuals (Boyd 1993) and is one of the most genetically diverse pinnipeds (Stoffel ). Therefore, given that LD is a function of both recombination rate and population size (Hill 1981), the rapid decay of LD in this species might also be a reflection of high long-term effective population sizes.

To provide further insights into the recovery of Antarctic fur seals globally, we quantified population structure across the species geographic range. Microsatellite genotypes provided evidence for two major geographic clusters, the first corresponding to South Georgia, the South Shetlands and Bøuvetoya, and the second corresponding to Kerguelen, Heard and Macquarie Island. By contrast, the RAD data uncovered an additional level of hierarchical structure, resolving South Georgia, the South Shetlands and Bøuvetoya as distinct populations. This is consistent with simulation studies suggesting that thousands of SNPs should outperform small panels of microsatellites at resolving population structure (Haasl and Payseur 2011) as well as with more recent empirical studies that have directly compared microsatellites with SNPs (Rašić ; Vendrami ). Furthermore, many of our populations had sample sizes of around five individuals yet could still be clearly distinguished from one another. This is in line with a recent simulation study suggesting that sample sizes as small as four individuals may be adequate for resolving population structure when the number of markers is large (Willing ). Thus, our results have positive implications for studies of threatened species for which extensive sampling can be difficult but where understanding broad as well as fine-scale population structure is of critical importance.

It is generally believed that Antarctic fur seals were historically extirpated from virtually all of their contemporary breeding sites across the sub-Antarctic, with the possible exception of Bøuvetoya, where sealing expeditions were more sporadic (Christensen 1935) and around a thousand breeding individuals were sighted just a few decades after the cessation of hunting (Olstad 1928). South Georgia was the first population to stage a major recovery, probably because a number of individuals survived at isolated locations inaccessible to sealers around the South Georgia mainland (Bonner 1964). Consequently, several authors have speculated that emigrant individuals from the expanding South Georgia population may have recolonized the species former range (Boyd 1993; Hucke-Gaete ). However, Wynen resolved two main island groups with mtDNA, while Bonin found that significant differences between the South Shetland Islands and South Georgia with microsatellites. Our results build on these studies in two ways. First, the two major clusters we resolved using both microsatellites and RAD sequencing are identical to those identified by Wynen , suggesting that broad-scale population structure is not simply driven by female philopatry but is also present in the nuclear genome. Second, within the Western part of the species range, we not only found support for the South Shetlands being different from South Georgia, but also Bøuvetoya, suggesting that relict populations probably survived at all three of these locations. By contrast, no sub-structure could be detected within the Eastern part of the species range, which taken at face value might suggest that a single population survived sealing in this region. Consistent with this, historical records suggest that fur seals went locally extinct at Heard and Macquarie islands (Page ; Goldsworthy ) and these populations may therefore have been recolonized by surviving populations in the Kerguelen archipelago. Thus, our study highlights the importance of relict populations to species recovery while also providing some evidence in support of local extinctions having occurred.

Delving a level deeper, we investigated individual variation in the form of inbreeding. A recent meta-analysis has shown that small panels of microsatellites are almost always underpowered to detect variation in inbreeding (Miller ). By contrast, a handful of recent studies have shown that tens of thousands of SNPs are capable of accurately quantifying inbreeding (Hoffman ; Huisman ; Bérénos ; Chen ; Kardos ). While empirical studies to date have largely focused on small, isolated populations where inbreeding may be common, it is less clear how prevalent inbreeding could be in larger, free-ranging populations. We found several lines of evidence in support of inbreeding in fur seals. First, was significantly positive indicating identity disequilibrium within the sample of individuals. Second, the variance of all genomic inbreeding coefficients was found to lie within the 95% confidence intervals of and therefore we can expect our estimates to reflect the realized level of inbreeding in the population. Third, all inbreeding coefficients were strongly intercorrelated, suggesting that our markers are uncovering consistent information about variation in genome-wide homozygosity caused by inbreeding.

Our results are in some respects surprising given that Antarctic fur seals number in the millions and are free-ranging and highly vagile. However, the species is also highly polygynous, with a handful of top males fathering the majority of offspring (Hoffman ) and females exhibiting strong natal site fidelity (Hoffman and Forcada 2012), which could potentially lead to matings between individuals with shared ancestry. As demographic effects can also generate variance in inbreeding sensu lato, we also cannot discount the possibility that the historical bottleneck contributed toward the variation we see today. To test this, we would need to quantify the length distribution of runs of homozygosity, which would require denser SNP data.

Our work builds upon another recent study that used RAD sequencing to quantify inbreeding in wild harbor seals (Hoffman ) where a higher estimate of was found, indicative of a greater variance in inbreeding within the sample. However, the study focused on stranded seals, many of which died of lungworm infection and may therefore have been enriched for unusually inbred individuals. In the current study, pups were sampled at random from within a single breeding colony, together with their parents. Consequently, our sample should be more representative of the underlying distribution of inbreeding within the population. In line with this, our estimate of is more similar to those obtained in wild populations of other polygynous mammals such as Soay sheep and red deer (Bérénos ; Huisman ).

Our results are consistent with previous studies documenting HFCs for numerous traits in the South Georgia population (Hoffman ; 2010; Forcada and Hoffman 2014) and suggest that these may well reflect inbreeding depression. More generally, literally hundreds of studies have documented HFCs across the animal kingdom (Coltman and Slate 2003) and it has been strongly argued that these HFCs are highly unlikely to occur when there is no variance in inbreeding (Szulkin ). The fact that we found variation in inbreeding in a large, free-ranging population is consistent with this notion and therefore contributes toward a growing body of evidence suggesting collectively that inbreeding could be more common in wild populations than previously thought.

Conclusion

We have generated an improved genome assembly for an important Antarctic marine top predator and used RAD sequencing to provide diverse insights from the level of the species through the population to the individual. Focusing on the larger South Georgia population, we characterized rapid LD decay and uncovered significant variation in individual inbreeding, while population-level analyses resolved clear differences among island groups that emphasize the importance of relict populations to species recovery. RAD sequencing and related approaches might conceivably be applied to other wild species to characterize patterns of LD decay, elucidate fine scale population structure and uncover the broader prevalence of inbreeding and its importance to wild populations.