DNA was prepared for long-read sequencing with Pacific Biosciences (PacBio; Menlo Park, California, USA) SMRTbell Library kits, following a manufacturer’s protocol [38]. Continuous long-read (CLR) sequencing was performed on 13 SMRT cells for a 10-hour movie on the PacBio Sequel at the BYU DNA Sequencing Center (DNASC). Short-read sequencing was performed in rapid run mode for 250 cycles in one lane on the Illumina (San Diego, California, USA) Hi-Seq 2500 at the DNASC after sonication with Covaris (Woburn, Massachusetts, USA) Adaptive Focus Acoustics technology and preparation with New England Biolabs (Ipswich, Massachusetts, USA) NEBNext Ultra II End Repair and Ligation kits with adapters from Integrated DNA Technologies (Coralville, Iowa, USA).

RNA was prepared with Roche (Basel, Switzerland) KAPA Stranded RNA-seq kit, following manufacturer recommendations. Paired-end sequencing was performed in high output mode for 125 cycles on the three samples together in one lane on the Illumina Hi-Seq 2500 at the DNASC.

DNA was prepared with Phase Genomics (Seattle, Washington, USA) Proximo Hi-C kit (Animal) using the Sau3AI restriction enzyme (cut site: GATC) following recommended protocols. Paired-end sequencing was performed in rapid run mode for 250 cycles in one lane on the Illumina Hi-Seq 2500 at the DNASC.

Double digest restricted site-associated (ddRAD) sequencing was used to measure intraspecific genetic variation across six sampling localities in the SWIO. We extracted total DNA using Qiagen DNeasy Tissue kits according to the manufacturer’s protocol (Qiagen, Inc., Valencia, California, USA). We visually examined the quality of DNA extractions using gel electrophoresis and by quantifying isolated DNA using a Qubit fluorometer (Life Technologies, Carlsbad, California, USA).

We modified a protocol developed by Peterson et al. [39] to prepare samples for ddRAD sequencing. We used the rare cutter PstI (5′-CTGCAG-3′ recognition site) and common cutter MspI (5′-CCGG-3′ recognition site). We carried out double digests of 150–200 ng total DNA per sample using the two enzymes in the manufacturer’s supplied buffer (New England Biolabs) for 8 hours at 37 °C. We randomly distributed samples from different localities across the sequencing plate to minimize bias during library preparation. We visually examined samples using gel electrophoresis to determine digestion success, then ligated barcoded Illumina adapters to DNA fragments [39]. After ligation, we pooled samples into 12 libraries and performed a clean-up using the QIAquick PCR Purification kit. We then performed the polymerase chain reaction (PCR) using Phusion Taq (New England Biolabs) and Illumina-indexed primers [39]. Library DNA concentration was checked using a Qubit fluorometer, followed by normalization, a second round of pooling into four libraries, and an additional QIAquick cleanup step. We then re-measured DNA concentration using a Qubit and combined equal amounts from each of the four pools into one. We analyzed this final pool using a BioAnalyzer (Agilent, Santa Clara, California, USA) and performed size-selection using a Pippin Prep (Sage Science, Beverly, Massachusetts, USA), selecting for fragments between 300 and 500 bp. This was followed by a final measure of concentration using a BioAnalyzer. We sent the library to the University of Oregon Genomics and Cell Characterization Core Facility where concentrations were verified via quantitative PCR (qPCR) before 100-bp single-end sequencing on an Illumina Hi-Seq 4000.

Illumina whole-genome sequencing (WGS) reads were corrected using Quake v0.3.5 (RRID:SCR_011839) [40], which depended upon old versions of R v3.4.0 (RRID:SCR_001905) [41] and the R package VGAM v0.7-8 [42, 43]. Quake attempts to automatically choose a k-mer cutoff, traditionally based on k-mer counts provided by Jellyfish (RRID:SCR_005491) [44]. To generate q-mer counts instead of k-mer counts, BFCounter v0.2 (RRID:SCR_001248) [45] was used. Quake suggested a q-mer cutoff of 2.33, which was subsequently used by the correction phase of Quake. Unlike the WGS reads, the Illumina DNA reads created with the Hi-C library preparation were not corrected.

An estimate of the number of k-mers present in the reads is required to run BFCounter. This number is based on the number of reads, the length of the reads, and k-mer size according to this equation: where n is the number of reads, l is the read length, k is the k-mer size, and T is the total number of k-mers (not necessarily unique or distinct) present in the reads. This assumes a uniform read length. If the reads are paired-end, n is still the number of reads, not the number of pairs of reads. Since ntCard v1.0.1 (RRID:SCR_022010) [46] was used to quickly get a picture for the k-mer coverage histogram, its reported value F1 was used instead of the equation as it is an estimate for T.

Illumina RNA-seq reads were corrected using Rcorrector v1.0.2 (RRID:SCR_022011) [47]. Rcorrector automatically chooses a k-mer cutoff based on k-mer counts provided by Jellyfish [44], which Rcorrector runs automatically for the user. Alternately, Jellyfish can be run externally or bypassed with an alternate k-mer counting program, and counts can subsequently be provided to Rcorrector, which may be started at “stage 3”. We bypassed Jellyfish by using BFCounter v0.2 [45] to count k-mers. Note that Rcorrector made no changes to the reads.

Several methods were attempted to correct the PacBio CLRs. Corrected reads from each method that did not fail were assembled, and the assembly results were used to choose the correction strategy. Ultimately, a hybrid correction strategy was employed. Typically, a “hybrid” correction strategy is one in which more than one data type (i.e., PacBio CLRs and Illumina short reads) are employed. This differs from a “self” correction strategy in which only the PacBio CLRs are used to correct themselves. Our hybrid strategy is not fully described by the word “hybrid” because both a “self” and “hybrid” strategy were serially employed; we referred to this strategy as a “dual” strategy (Figure 2). First, the reads were self-corrected using Canu v1.6 (RRID:SCR_015880) [48]. Second, the self-corrected reads were further corrected using Illumina short-reads (previously corrected with Quake) using CoLoRMap (commit #baa680, RRID:SCR_022012) [49]. Note that Illumina short reads had to be interleaved into a single file for CoLoRMap. Similarly, all PacBio reads were concatenated into a single file, and the headers were required to be unique up to the first space. CoLoRMap is a pipeline with a basic wrapper script. In practice, it makes more sense to run each step in the wrapper script as separate jobs to avoid re-computing if a failure (e.g., too much random access memory (RAM) or time) occurs in a downstream step; this is how we ran the pipeline.

Before conducting population genomic analyses, we performed outlier tests to identify loci putatively under selection. These are usually identified by a significant difference in allele frequencies between populations [97]. Specifically, we implemented two outlier detection methods that accommodate missing data: pcadapt v4.1.0 (RRID:SCR_022019) [97] and BayeScan v2.1 (RRID:SCR_022018) [98]. The assumption behind pcadapt is that loci associated with population structure, ascertained via principal component analysis (PCA), are under selection [97]. pcadapt is advantageously fast and able to handle many loci. The number of principal components (K) was chosen based on visualization of a scree plot of the eigenvalues of a covariance matrix. Once the K-value was chosen, the Mahalanobis distance (D test statistic) was calculated using multiple linear regression of the number of SNPs versus K [97, 99]. To account for false discovery rates, the p-values generated using the Mahalanobis distance D were transformed to q-values using the R v3.6.3 [41] q-value package v2.15.0 (RRID:SCR_001073) [100] with the cut-off point (𝛼) set to 10% (0.1).

BayeScan measures allele frequencies between different populations and identifies loci perceived to be undergoing natural selection based on their FST values [101, 102]. The method applies linear regression to generate population- and locus-specific FST estimates and calculates subpopulation FST coefficients by taking the difference in allele frequency between each population and the common gene pool. BayeScan incorporates uncertainties in allele frequencies owing to small sample sizes, as well as imbalances in the number of samples between populations [98]. We assigned each of the six sampling localities as a population. Our analyses were based on 1:50 prior odds and included 100,000 iterations and a false discovery rate of 10%. We used the default values for the remaining parameters and visualized results in R v3.6.3 following the developer’s manual [103]. After running both pcadapt and BayeScan, we used R to assess the number of outliers identified by both programs and subsequently removed outlier loci to generate a neutral dataset for downstream analyses. The outlier dataset results are presented in Figure 4. These datasets were also annotated with respect to the MAKER-based annotations of the genome assembly using SnpEff v5.0e (RRID:SCR_005191) [104].

To examine population structure, we used a model-based clustering method to reconstruct the genetic ancestry of individuals using sparse non-negative matrix factorization (sNMF) and least-squares optimization. Model-based analyses were performed using the package LEA v2.6.0 (RRID:SCR_022020) [105] in R. The sNMF function in LEA estimates the number of ancestral populations and the probability of the number of gene pools from which each individual originated by calculating an ancestry coefficient and investigating the model’s fit through cross-entropy criterion [106]. We calculated and visualized cross-entropy scores of K population clusters ranging from 1–10 with 10 replicates. To complement sNMF, we also used PCA, a distance-based approach based on variation in allele distributions, implemented in VCFtools v0.1.16 [69]. For sNMF and PCA analyses, no populations were assigned a priori. We assigned each of the six sampling localities as populations for subsequent visualization, grouped into four “island groups” based on the proximity of some of the atolls that comprised our sampling localities (Figure 5). The five Seychelles atolls we sampled were spread among three separate clusters of islands, which are commonly referred to as the “outer” island groups owing to the geographic locations of these outlying coralline islands relative to the densely populated, granitic “inner” islands of the Seychelles archipelago. The island groups consisted of (i) Amirantes (St. Joseph’s Atoll), (ii) Farquhar (Farquhar and Providence Atolls), (iii) Aldabra (Aldabra and Cosmoledo Atolls), as well as (iv) Mauritius (St. Brandon’s Atoll; Table 1). We computed summary statistics in R v3.6.3, including pairwise FST estimates (StAMPP v1.6.1 (RRID:SCR_022022) [107]), isolation by distance via the Mantel Rand test (adegenet v2.1.3 (RRID:SCR_000825) [108]), and expected and observed heterozygosity (hierfstat v0.5-7 (RRID:SCR_022021) [109]) to compare genetic diversity and differentiation between the four island groups.

Paired-end, short-read sequencing (Illumina) yielded 109.5 million pairs of reads comprising 53.86 Gbp. The mean and N50 read lengths were 245.981 and 250, respectively. Continuous long-read sequencing (PacBio) generated 9.5 million reads with a total of 69.85 Gbp. The mean and N50 read lengths were 7,352.726 and 13,831, respectively. The longest read was 103,889 bp. The read length distribution is plotted in Figure 5. Result summaries for both sequencing runs are available in Table 5.

RNA-seq from the three tissues (i.e., gill, heart, and liver) generated 270.7 million pairs of reads totaling 67.2 Gbp. The gill tissue yielded 107.7 million pairs of reads, with a total of 26.7 Gbp. The heart tissue generated 19.6 Gbp across 78.8 million pairs of reads. The 84.2 million pairs of reads from the liver tissue comprised 20.9 Gbp. Across all three tissues, the mean and N50 read lengths were 124.122 and 125, respectively. The combined results from all three tissues are summarized in Table 5.

Sequencing yielded 88.7 million pairs of reads comprising 44.28 Gbp. The mean and N50 read lengths were 249.493 and 250, respectively. A summary of these results is presented in Table 5.

After data processing using ipyrad, we recovered a mean of 114,324 reads per individual and an average of 107,105 loci per individual. Following filtering for missing data, minor allele frequency, and linkage disequilibrium, the dataset contained 66 individuals and 38,355 SNPs. BayeScan, being a more conservative outlier detection method than pcadapt, did not identify any outliers; we thus used only outlier detection results from pcadapt. Subsequent removal of pcadapt outliers (N = 155) resulted in a neutral dataset containing 38,200 SNPs with 9% missing data.

When Quake corrects paired-end reads, three outcomes are possible for each pair of reads: (i) both reads are either already correct or correctable, (ii) one read is either correct or correctable and the other is low-quality, or (iii) both reads are low-quality. Of the original 218.96 million reads (109.5 million pairs of reads), Quake corrected 62.7 million of them and removed 51.6 million of them. 5.97 million pairs of reads were discarded because both reads were rated as erroneous. 39.6   million pairs of reads had one read that was correct or correctable and one read that was low-quality; these were also discarded. The remaining 63.88 million pairs of reads were either correct or correctable and were kept in the final read set containing 29.11 Gbp of sequence.

No corrections were made to the RNA-seq reads by Rcorrector [47].

The dual-correction strategy (self-correction followed by hybrid-correction) reduced the number of reads from 9.5 million to 2.79 million and the total number of bases from 69.85 Gbp to 36.79 Gbp. The mean and N50 read lengths were changed from 7,354 and 13,831 to 13,193 and 15,483, respectively. The longest read was 63,271 bases. The distribution of read lengths can be viewed in Figure 6.