Long-term aggregation of larval fish siblings during dispersal along an open coast

Study System.

Oceanographic conditions off the Oregon Coast in the US Pacific Northwest are characterized by an offshore, southward-flowing current (California Current) and an inshore, northward-flowing current (California Undercurrent) that has seasonal flow inversions in the top 50 m of water (39). Seasonal wind-driven upwelling occurs throughout spring and summer, with intermittent periods of relaxation that are often accompanied by an inshore flow of surface currents, and large recruitment events of intertidal and nearshore organisms (29, 40, 41). We obtained Bakun index values of daily upwelling at 45 °N, 125°W from the Pacific Fisheries Environmental Laboratory division of the National Oceanic and Atmospheric Administration (NOAA-PFEL). This index represents daily averages of wind-driven cross-shore transports computed from Fleet Numerical Meteorology and Oceanography Center 6-hourly surface pressure analysis. Positive values indicate offshore transport in units of cubic meters per second along each 100 m of coastline.

Study Species.

Splitnose rockfish (Sebastes diploproa) occur from Alaska to Baja California, and adults are most abundant at depths of 100–350 m. They live up to 103 y and have a relatively late (6–10 y) age of maturity, and females can produce up to 255,000 eggs per brood (42). Similar to other rockfish species, female splitnose are live-bearing: eggs hatch in the maternal ovary several days before extrusion, and larvae are capable of swimming on release (43). Females can store sperm from multiple males for several months before fertilizing the eggs, thereby releasing offspring with multiple paternity (44). Larvae and early juveniles can remain in the pelagic environment for up to 1 y before settling to a benthic habitat at a size of 30–50 mm (42). Pelagic juveniles frequently aggregate to drifting kelp mats (45), and the transition to benthic habitats occurs later in northern than southern latitudes, peaking in May–June in Southern California and in August–September in Oregon and Washington. The long pelagic duration of benthic species similar to splitnose rockfish is thought to enable larvae and pelagic juveniles to disperse large distances from the parental source (7). The reproductive output of the splitnose population is estimated to be ∼10 billion larvae per year, and there is no evidence of spatial population structure, likely because of their long pelagic larval duration and longevity (46).

Recruits.

We collected newly settled splitnose rockfish recruits near Depoe Bay, central Oregon, using standardized monitoring units for the recruitment of fishes (47). These collectors are made of black polyvinyl chloride mesh folded inside garden fencing that is shaped in a long cylinder (100 × 30 cm). This creates a 3D structure that simulates natural recruitment substrates such as a kelp canopy. We deployed seven replicate collectors at 1 m below the surface at sites where the depth was ∼15 m. Collectors were located 390–1,170 m offshore and were 425–1,315 m apart (Fig. 1). From April 20 to September 11, 2013, we collected newly recruited fish approximately every 2 weeks, using hinged butterfly-style nets to enclose each standardized monitoring unit and collect the recruited fishes. Collected fish were killed with MS-222, measured with calipers to the nearest millimeter, and stored at −80 °C. An exceptionally large (n = 538) pulse of splitnose recruits was collected on September 11 (Fig. 2). A fin-clip of each recruit was stored in 95% (vol/vol) ethanol for the kinship analysis. According to a size-at-age regression of juvenile splitnose rockfish (48), these recruits ranged in age from 120 to 180 d.

Reference Population.

Unbiased measures of genetic relatedness require a reference population of unrelated and noninbred individuals (49, 50). This assumption is violated when the focus population (i.e., fish recruits of the same cohort) is used as a reference population (50). Therefore, we used adult splitnose from nine different locations off the Oregon Coast as a reference population (Fig. S1). Tissue samples from 144 adult splitnose rockfish were collected during the 2015 NOAA West Coast Groundfish Bottom Trawl Survey in Aberdeen net trawls (15 min at ∼2.2 nm⋅hr⁻¹) at depths of 117–225 m and stored in 95% (vol/vol) ethanol.

Microsatellite Genotyping.

We genetically analyzed 513 splitnose recruits and 95 individuals from the reference population. The remaining 25 recruit samples were suspected of cross-contamination in the laboratory processing, and therefore removed from the analysis. The 95 reference population samples represent all nine sampling locations (Fig. S1). Within location, samples were randomly selected for the analysis. DNA was extracted following a standard silica-based method for all samples (51).

We optimized PCR conditions of 35 microsatellite markers isolated from congeneric species, and screened the loci using 16 of our splitnose recruits. We successfully amplified 24 markers that had multiple alleles per locus. We then selected the 20 most polymorphic markers (Table S1) to genotype all of the remaining samples. PCR products were visualized using an ABI 3730xl DNA Analyzer, and alleles were sized using GENEMAPPER SOFTWARE 5 (Applied Biosystems).

We adopted a conservative approach and removed all of the recruit samples that had any missing loci, leaving 491 (96%) fully genotyped recruits (at 20 loci) to conduct the kinship analysis. One sample from the reference population that was missing >10 loci was also removed, leaving 94 total samples (99%). We regenotyped 95 (19%) randomly selected recruits, and calculated the genotyping error rate by dividing the number of discordant alleles by the total number of scored alleles. The genotyping error rate across all 20 loci was 2% ± 2% (mean ± SD; Table S1). HWP and LD were tested independently for the recruit and reference samples, using GENEPOP 4.2 (52). Significance (P < 0.05) of HWP and LD was tested after Bonferroni adjustments (53). We evaluated the effect of the sample size affecting the power to detect departures from HWP by applying a paired t test on the coefficient of inbreeding (F_IS). The mean F_IS was not significantly different between the reference population and the recruits (t₁₉ = −1.22; P = 0.24), indicating that the power to detect loci out of HWP is greater in the recruits. However, F_IS was positive in both cases, indicating that both recruits and reference population have deficient heterozygosity. Presence of null alleles is frequently associated with this condition. Using MICRO-CHECKER (54), we detected the presence of null alleles in nine loci of the recruits and three loci of the reference population (Table S1). Thus, presence of null alleles may contribute to a deficiency of heterozygosity in the recruits. However, null alleles alone are unlikely to account for similar deficiencies in the reference population. To test whether the deficiency of heterozygosity in the reference population was the result of a Wahlund effect, adult genotypes were grouped into northern (n = 46) and southern (n = 48) sample sites, respective to 44 °N (Fig. S2), and we calculated the index of fixation (F_ST) between the groups using FSTAT v2.9.3.2 (55). After 1,000 permutations, we found evidence of a small but significant amount of genetic differentiation (F_ST = 0.002; P = 0.042) that, along with the presence of null alleles, likely contributed to the heterozygote deficiency.

Data Analysis.

Pairwise relatedness (r) was measured with the Triadic IBD estimator of relatedness, using COANCESTRY 1.0.1.5 (18, 56), which applies a maximum likelihood method that estimates pairwise relatedness using the genotype of a third individual as a control. This estimator measures relatedness effectively (16, 18, 57, 58) and is robust to minor error rates, null alleles, deviation of HWP, and LD, as well as small amounts of genetic structure (18). Allele frequencies of the adult population were used as reference in COANCESTRY (50). The number of missing alleles in the reference population relative to the recruits ranged from 0 to 15 (Table S1). This was expected because the sample size of the reference population was smaller than the sample size of the recruits (59). Thus, we assigned a low (≤0.005) frequency value to all missing alleles in the reference population. Although this approach might introduce marginal bias for r estimates, it is likely a conservative bias, as most of the assigned frequencies were higher than those observed in the recruits, thus underestimating r (18).

Testing the estimator performance on simulated data is a critical process to predict true relatedness based on measures of r (18, 58, 60). We therefore tested our analytical power by measuring r in simulated datasets composed of individuals with defined relationships. A dataset of 1,000 half-sibling pairs and 1,000 full-sibling pairs (sharing one and both parents, respectively) was used to calculate the expected r estimates among individuals and provide criteria that could be used to distinguish related individuals from unrelated individuals. To predict the number of expected false-positive siblings at a given r threshold, we calculated r in 100 datasets composed of n = 491 unrelated individuals (120,295 pairwise comparisons) and averaged the number of false pairs falling above alternate thresholds. The number of simulated unrelated fish was selected to match the sample size of analyzed recruits. We additionally simulated 1,000 first-order cousin pairs to estimate the number of false-positive siblings that might be attributed to cousin relationship. All these simulations were conducted with COANCESTRY, using the allele frequencies of the 491 completely genotyped splitnose recruits. We found that r = 0.35 was an optimal cutoff that minimized the number of false-positive pairs (unrelated individuals with an r ≥ 0.35) while still preserving power to identify true sibling pairs (Fig. S2). The 100 simulations of n = 491 unrelated individuals yielded ∼6 ± 2 (mean ± SD) false-positive pairs, which is a low proportion (4.6 × 10⁻⁵) considering the 120,295 possible pairwise comparisons of each dataset. Using the same threshold (r ≥ 0.35), at least 16% of the simulated half-siblings and 90% of the simulated full-siblings were detected, whereas less than 1% of the simulated cousin pairs fell above the threshold. Thus, we concluded that most of the pairs observed above a cutoff of r = 0.35 are likely to be siblings and not cousins. A lower r cutoff (r = 0.3) resulted in a higher proportion of false positives (36 ± 7; mean ± SD), hampering the discrimination of real siblings from false positives. A higher r cutoff (r = 0.4) reduced the detection of real siblings (5% half-siblings; 78% full-siblings), resulting in an increased number of false negatives. We performed a similar power analysis with six other commonly used r estimators implemented in COANCESTRY (56), and found that the Triadic IBD estimator was the most precise and conservative. We tested the performance of COLONY (61) and ML-Relate (62) on simulated data as alternative methods of sibling identification. COLONY produced ∼750 false-positive siblings and only identified 37% of the 1,000 simulated full-sibling pairs correctly; 58% were incorrectly identified as half-siblings, and the remaining 5% were classified as unrelated. ML-Relate performed slightly better, detecting 88% and 85% of 1,000 simulated full- and half-sibling pairs, but producing 15 false full-sibling pairs and ∼6,300 false half-sibling pairs. Consistent with these results, only 27 of the 695 sibling pairs identified by COLONY among the recruits coincided with the siblings identified with the Triadic IBD estimator, whereas all of the 74 sibling pairs identified by the Triadic IBD estimator were found within the 9,198 sibling pairs identified with ML-Relate (Fig. S3). On the basis of these simulations, using the Triadic IBD estimator appeared to be the most conservative approach.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S3.

Number of siblings identified among 491 corecruiting splitnose rockfish using Colony (light gray-filled circle; 695 pairs), ML-relate (white-filled circle; 9,198 pairs), and the Triadic IBD estimator of relatedness (dark gray-filled circle; 74 pairs). Circles are not drawn to scale, and values indicate number of sibling pairs. Coinciding sibling pairs identified by multiple approaches are indicated by overlapping circles. There is relatively little overlap among the sibling pairs identified using Colony and the Triadic IBD estimator. A power analysis based on simulated data indicated that the Triadic IBD estimator performed best.

As a consequence of siblings sharing common parents, a relatively small effective number of breeders would be expected in a recruitment cohort containing siblings. Using LDNE v2.0 (63), a program that estimates the effective population size based on linkage disequilibrium information, we estimated the effective number of breeders for the entire sample collection (n = 491), the putative siblings (n = 57), and the collection minus the putative siblings (n = 434). The effective number of breeders is expected to be larger in the 434-recruit group than in the entire collection because we removed most of the linkage disequilibrium signature induced by siblings sharing common parents.

All analyses were performed in R v3.2.1 (64), using packages ‘plyr’ v1.8.3 (65) and ‘reshape2’ v1.4.1 (66). Figures were created using ARCGIS 10.2 (ESRI 2014), and R packages ‘ggplot2’ v2.0.0 (67) and ‘venneuler’ v1.1–0 (68). The network visualization was created with R package ‘igraph’ v1.0.1 (69), which distributes each node (i.e., fish) maximizing the occupied space of the figure while minimizing crossing-over of connected nodes.