Phylogeography of the Pacific Sardine, Sardinops sagax, across its Northeastern Pacific Range

Analysis of population structure is important to the study of species and to the field of evolutionary biology. Molecular analysis of populations in relation to hypothesized barriers to gene flow can elucidate questions of population structure. Biogeographic boundaries occur between areas with unique assemblages of fauna, and at times, they may align with intra-species phylogeographic breaks, indicating concordance of a biogeographic boundary with a barrier to gene flow (Avise 2000). A question that remains is how frequently this concordance occurs in oceans where barriers to dispersal are less obvious than in terrestrial systems.

The marine ecosystems of the northeastern Pacific Ocean, ranging from the warm waters of the Gulf of California north to the colder waters of Alaska are home to a diverse assemblage of marine organisms. This region is represented by several distinct assemblages of marine fishes representing multiple biogeographic or faunal provinces. Traditionally, these provinces have been described as the Aleutian province, extending from Nunivak Island in Alaska through the Dixon Entrance (U.S.-Canada border), the Oregonian province from the Dixon Entrance southward to Point Conception (Briggs and Bowen 2012 consider the southernmost boundary to be the Los Angeles region), the Californian Province extending from Point Conception (Briggs and Bowen 2012 consider the northernmost boundary to be the Los Angeles region) to Bahía Magdalena in Mexico, and the Cortez Province from Bahía Magdalena around the southern end of the Baja California Peninsula into the Gulf of California (Briggs 1974; Horn and Allen 1978; Boschi 2000; Hastings 2000; Briggs and Bowen 2012).

Arguably the most well-studied of the traditional biogeographic breaks within the range of the Pacific sardine is Point Conception, which has long been recognized as a biogeographic break separating the cold Oregonian and temperate Californian fish assemblages (Briggs 1974; Horn and Allen 1978; Boschi 2000). The waters surrounding this prominent headland are characterized by complex oceanographic patterns resulting from the convergence of the warmer surface waters of the Southern California Bight and the colder surface waters north of the point (Wang 1997; Hohenlohe 2004). Historically, these oceanographic differences have been considered to act as a barrier to dispersal of adult and larval organisms. Briggs and Bowen (2012), however, argue that this area is more of a transition zone than a boundary. With this in mind, it is interesting to ask if the same processes acting on species distributions in this region also act on an intraspecific level. That is, is Pt. Conception a barrier to gene flow (i.e., a phylogeographic break driving population separation), or does it act more as a diffusive membrane (i.e., allowing sufficient gene flow to prevent evolutionary divergence)? While some studies have shown the former, others have shown the latter (Burton 1998; Bernardi 2000; Dawson et al. 2001; Craig et al. 2011).

The geographic range of the Pacific sardine also includes other hypothesized barriers to dispersal. On the western coast of the Baja California Peninsula, the northern end of the Cortez province is marked by a strong temperature gradient at Bahía Magdalena (Palacios-Salgado et al. 2012). The biogeographic boundary at Punta Eugenia has also been implicated as a phylogeographic boundary zone for marine organisms, with many examples of nearshore species in the Eastern Pacific showing reciprocal monophyly and strong genetic divergence among disjunct populations in the Gulf of California (Terry et al. 2000; Bernardi et al. 2003). The Gulf of California is considered as part of the Cortez faunal province (Hastings 2000; Briggs and Bowen 2012; Palacios-Salgado et al. 2012), and, as with the boundary of the Oregonian and Californian provinces, some discrepancy exists on the exact delineation of its endpoints (Briggs and Bowen 2012).

In marine fishes, the primary driver of gene flow is dispersal, either as adult movement or egg and larval drift between regions or groups. Hypothesized barriers to dispersal such as Point Conception and Punta Eugenia may be more relevant to benthic, estuarine, or live-bearing fish with lower dispersal ability by acting as more effective barriers to gene flow over long time scales (Bernardi 2000; Bernardi and Talley 2000; Buonaccorsi et al. 2005; Johansson et al. 2008; Craig et al. 2011). The phylogeography of demersal fishes on the west coast of North America is broadly represented in the literature (e.g., Bernardi 2000; Bernardi and Talley 2000; Terry et al. 2000; Dawson et al. 2001; Cope 2004; Craig et al. 2011), yet studies on pelagic species are less so. Additional investigation of population genetic data for pelagic species in the northeastern Pacific Ocean may facilitate an understanding of biogeographical patterns and gene flow in marine fishes.

Pacific sardine (Sardinops sagax) are clupeid fish that are at times the dominant forage fish in the California Current (Preti et al. 2001; Emmett et al. 2005; Lo et al. 2011; McClatchie et al. 2017). Pacific sardine spawn at various intensities year-round through a wide range of temperatures often with peaks in activity during the spring (Ahlstrom 1954; Ahlstrom 1959; Ahlstrom 1965; Hernandez-Vazquez 1994; McFarlane et al. 2005). Recent examination of CalCOFI and IMECOCAL data across Alta and Baja California from 1963-2015 showed that sardine larvae were present across all seasons but were typically at their lowest abundance during winter (Tran 2023). Sardine also spawn in the Gulf of California during the cooler winter months (Hammann et al. 1998; Aceves-Medina et al. 2009). Sardine eggs are planktonic (Zweifel and Lasker 1976), and after hatching larvae absorb the yolk as they develop (Lasker 1965) and may take between 5 and 8 d to deplete the yolk sack depending on the oceanic conditions (Politikos et al. 2018). In the California Current Ecosystem (southern Canada to southern Baja California), the range for suitable egg and larval survival is thought to be 13–16 °C (Lasker 1964; Agostini et al. 2007; Dorval et al. 2019). Juvenile sardines are strict filter feeders (Arthur 1976), but as adults they also exhibit particle feeding (Emmett et al. 2005). Tagging studies have shown that at least a portion of Pacific sardine biomass migrates extensively each year, moving northward in the spring and returning to warmer southern waters in the fall (Clark and Janssen 1945; Lo et al. 2011; Dorval et al. 2019). These studies have shown that older, larger fish tend to migrate farther than smaller, younger fish (Clark and Janssen 1945). Parasite data have also suggested that not all individuals make an annual migration (Baldwin 2010; Baldwin et al. 2012; Jacobson et al. 2019).

Previous genetic studies have given no indication of deep genetic structure in Pacific sardine (Hedgecock et al. 1989; Lecomte et al. 2004; García-Rodríguez et al. 2011). However, these studies analyzed comparatively few individuals and/or used methods that were common practice at the time of publication but are no longer widely used (e.g., allozymes). In an effort to further our understanding of the phylogeography of Pacific sardine and to understand the effect of hypothesized barriers to gene flow over evolutionary time scales in a coastal pelagic species, this study expands on previous research by including a larger sample size and a more fine-scale/broader geographic sampling range. Given the ecology and high motility of the Pacific sardine, it was considered unlikely that hypothesized barriers to gene flow would act as long-term barriers to gene flow in Pacific sardine along their northeastern Pacific range.

Tissue samples from the USA were sourced from the archival collection at the Southwest Fisheries Science Center in La Jolla, California, that were collected during surveys of coastal pelagic species aboard NOAA research vessels (Cutter and Demer 2008). Samples from Mexico were collected through collaborative port/fisheries sampling programs. All samples consisted of muscle tissue or fin clips stored in 100% ethanol at ambient temperature. Twelve sampling locations were selected. These included Vancouver Island (Canada), Washington State (USA), northern, central, and southern Oregon (USA), Crescent City (USA), Pescadero (USA), Monterey (USA), Morro Bay (USA), Point Conception (USA), San Diego (USA), Ensenada (Mexico), Isla Cedros (Mexico), Bahía Magdalena (Mexico), El Desemboque (Mexico), and Guaymas (Mexico) (Fig. 1). All samples were collected from 2002-2006. Total genomic DNA extraction was performed using the Qiagen DNeasy 96 Blood and Tissue Kit following manufacturer’s protocol with the exception of the final elution step which was performed using 50 µl of DI water.

A 1500 bp portion of the Cytochrome b gene was amplified in two ∼750 bp fragments using the polymerase chain reaction (PCR). Primers were designed in house using reference sequences available from GenBank (Table 1). PCR reactions consisted of 1 μl 10X buffer (670 mM Tris, 166 mM (NH₄)₂SO₄, 20 mM MgCl₂, 100 mM 2-mercaptoethanol, pH = 8.8), 1 μl dNTPs at 2 mM, 0.25 μl BSA at 20 mg/ml, 0.5 μl forward and reverse primer at 10 μM, 5.7 μl MilliQ water, 0.05 μl Taq 5 U/ul, and 1 μl DNA at 10-100 ng/μl. Thermal cycling parameters were: 94 °C for 2:00 min, 35 cycles of 94 °C for 30 sec, 47 °C for 30 sec, 72 °C for 1 min 15 sec, with a final extension at 72 °C for 5 min. PCR products were visualized by gel electrophoresis using 2% agarose gel with ethidium bromide to check for successful amplification. PCR products were purified by adding 6 μl of ExoSap solution (Exonuclease I at 3.33 U/μl and Shrimp Alkaline Phosphatase at 0.66 U/μl) and heating to 37 °C for 15 min followed by 80 °C for 15 min. Sanger sequencing was performed using BigDye Terminator v3.1 using the following master mix: 2 μl of 5x Sequencing Buffer, 1 μl BigDye Terminator v3.1 Ready Reaction Mix, 0.3 μl primer at 10 μM, 5.2 μl of MilliQ water, and 3.5 ul PCR product. Sequence reactions were cleaned using ethanol/salt precipitation and centrifugation. 86 μl of a high salt precipitation solution (3 ml of 3 M NaOA_C pH 5.2 ± 0.1, 22.5 ml of MilliQ H₂O, 62.5 ml of 100% EtOH) was added and samples were centrifuged for 20 min at 3700 rpm. Samples were drained and washed with 150 μl of 70% ethanol twice and resuspended in 10 μl of HiDi formamide. Samples were then denatured at 95 °C for five min, and sequenced on an ABI3730 Genetic Analyzer.

Contiguous consensus sequences were concatenated from forward and reverse strands, edited, aligned, and trimmed to 1062 bp using Sequencher v5.0 (Sequencher sequence analysis software, Gene Codes Corporation, Ann Arbor, MI USA). Those not reaching this size threshold were eliminated from further analysis. After size selection, 434 individuals were retained. Aligned sequences were analyzed for population structure using Arlequin v3.5. (Excoffier and Lischer 2010), and the R package strataG (Archer et al. 2016). All analyses done in R were performed using R Statistical Software v4.2.2 (R Core Team 2022). Haplotype and nucleotide diversities, as well as neutrality tests for Tajima’s D and Fu’s F_S were calculated in Arlequin by sampling site, and all tests were run with 10,000 permutations. A comparison of pairwise differences was computed in Arlequin with 1000 permutations with the significance level set to 0.05. Fixation indices (Φ_ST and Φ_CT) were calculated using strataG v2.5.01. Φ statistics were calculated for samples grouped into a Pacific group and a Gulf of California (GOC) group, with Bahia Magdalena tested in each group. Φ_ST was also tested for all samples pooled into one single group. Overall haplotype and nucleotide diversities, Tajima’s D, Fu’s F_S using Nei’s genetic distance, and the population parameter Θ_S were calculated in the R package pegas v1.1 (Paradis 2010). Estimates of time parameters and population parameters Θ_o and Θ₁ (Table 2) were fitted to observed mismatch distributions to determine coalescence times, following Li (1977) and Rogers and Harpending (1992). Coalescence analysis requires an estimate of generation time and mutation rate. Generation time was estimated following Depczynski and Bellwood (2006) using age of maturity at 1 year and a maximum age of 10 years. Divergence estimates for cytochrome b at c. 2%/Myr between species and 1% within lineages have been obtained for a number of reef fishes (e.g., Muss et al. 2001), and this rate was provisionally applied here in the absence of a known rate for sardine. Demographic history and coalescence time were also analyzed by producing a Bayesian Skyline Plot (BSP) created in Beast v.2.7.6 (Bouckaert et al. 2019) using a GTR substitution model with 4 gamma categories and an in-run estimated shape parameter. 500 × 10^8 generations were run discarding the first 50 × 10^8 as burn in. In order to visualize the sequence data, a minimum-spanning haplotype network was created using PopART v.1.7 (Leigh and Bryant 2015; Bandelt and Röhl 1999). Tracer v.1.7 (Rambaut et al. 2018) was used to create the Bayesian Skyline Plot.

On the one hand, overall haplotype diversity was large. Of the 434 individuals examined, 305 unique haplotypes were found (h = 0.992; GenBank accession #OR480110-OR480414). On the other hand, nucleotide diversity was low (Θπ = 0.0070; Table 3). The most common haplotype was found in 29 individuals and was distributed across the entire sampling range, with the exception of Bahía Magdalena (Fig. 2). The largest haplotype diversity was observed at Guaymas, whereas the smallest haplotype diversity was seen at Vancouver Island (Table 3). Both Tajima’s D and Fu’s F_S were negative and significant when measured site by site, and was also negative and significant overall (Table 3). Θ_S was estimated at 37.74943.

Estimated coalescence time using the traditional mismatch method was approximately 235,400 years before present. The Bayesian Skyline Plot, however, identified two periods of rapid population growth, one beginning ∼280 kya and a second beginning ∼110 kya (Fig. 3). Present day effective population size (Ne) was estimated at nearly 100 million.

Φ_ST, the measurement of variance within sampling sites, showed small but significant values, and Φ_CT, which measures variance among groups, was small and non-significant when Bahía Magdalena was included in the Pacific group, but was significant when included in the GOC group (Table 4). Pairwise Φ_ST values ranged from 0.0000 to 0.04487, with several significant but small values (Table 5).

Throughout the studied range in the northeastern Pacific Ocean, biogeographic boundaries align with phylogeographic breaks in some demersal species, however these patterns differ among species (Bernardi 2000; Bernardi and Talley 2000; Dawson et al. 2001; Craig et al. 2011; Chabot et al. 2015). Many demersal species experience limits to gene flow related to habitat associations, although few appear constrained by Point Conception as phylogeographic breaks seem to align with other geographic boundaries, especially for fishes that are restricted to specific habitats or are live-bearing (Bernardi 2000; Bernardi and Talley 2000; Terry et al. 2000; Dawson et al. 2001; Craig et al. 2011). For example, the California halibut (Paralichthys californicus) does not show phylogeographic breaks across its studied range despite its relatively sedentary adult behavior. Strong connectivity likely results from the approximately 30-day pelagic egg and larval stage conferring sufficient dispersal ability to cross barriers (Craig et al. 2011). Given the mobility of the Pacific sardine as a pelagic fish and its broad habitat requirements, Point Conception is an unlikely barrier to gene flow for sardine which is borne out by the data herein.

The Punta Eugenia biogeographic boundary is a biogeographic break that coincides with a phylogeographic break for several fishes (Bernardi et al. 2003) and is hypothesized to have been connected to the current day Gulf of California by an ancient Neogene seaway (Riddle et al. 2000). The closure of this seaway is thought to have resulted in the present-day disjunct distributions of several marine fishes (Bernardi et al. 2003). These disjunct distributions may also have resulted from more recent dispersal of organisms around Cabo San Lucas during times of glacial cooling that is not possible under current climate conditions (Riddle et al. 2000; Bernardi et al. 2003). This is due to the presence of warm waters at the southern extreme of the Baja California Peninsula, which act as a barrier to temperate species (Bernardi 2000; Bernardi and Talley 2000; Bernardi et al. 2003). The possibility for ongoing connection by migration in deeper water or in cooler months, as well as the mobility and tolerance of a broad temperature range makes it unlikely that this biogeographic boundary limits gene flow for Pacific sardine which is supported by the data herein.

The small divergences between Pacific sardine populations support a lack of barriers to long-time scale gene flow in this motile pelagic species when compared to demersal fishes in the same geographic range. Both overall (Table 4) and pairwise genetic distances (Table 5) suggest substantial gene flow across the Point Conception and Punta Eugenia biogeographic boundaries as well as evolutionary time scale connectivity between the Gulf of California and the Pacific Ocean. The genetic analysis of the Pacific sardine shows no evidence of population structure and high genetic diversity across the sampled range which is consistent with previous studies (Hedgecock et al. 1989; Lecomte et al. 2004; García-Rodríguez et al. 2011). This is supported by large overall haplotype diversity (h = 0.992) and the large number of unique haplotypes. Among the 434 individuals, there were 305 unique haplotypes and a majority of shared haplotypes were only observed twice (Fig. 2). Nucleotide diversity was also low (Θπ = 0.0070). Such patterns are consistent with previous findings for Pacific sardine, and indicate an evolutionarily recent population expansion of this population (Bowen and Grant 1997; Grant and Bowen 1998; Avise 2000). The coalescence time estimate of ∼235,400 yr before present indicates an evolutionarily recent coalescence and agrees with previous estimates (Bowen and Grant 1997; Grant and Bowen 1998; García-Rodríguez et al. 2011). Other studies on pelagic populations of fishes which have shown large panmictic populations have similarly high haplotype diversities (Lecomte et al. 2004; Theisen et al. 2008; Zheng et al. 2015; Min et al. 2019).

Neutrality tests also support the notion that the population may have undergone a large and evolutionarily recent expansion. As with high haplotype diversity, negative Fu’s F_S and Tajima’s D have been consistently found in recent studies of panmictic pelagic fishes (Lecomte et al. 2004; Zheng et al. 2015; Min et al. 2019). A negative D indicates an excess of low-frequency polymorphisms relative to neutral expectation, which may reflect rapid population growth and can mean that the population has undergone a bottleneck in its evolutionary past. Nevertheless, it is important to remember that a selectively neutral site could be linked to a site under selection and this could cause negative values to appear in the absence of past bottlenecks (Tajima 1989). As cytochrome b is mitochondrial, it is inherited as an entire unit with the rest of the mitochondrion and therefore could be under any form of selective pressure that other parts of the mitochondrial genome would be experiencing. Fu’s F_S was consistently negative and significant. Older mutations that occurred closer to the time of the most recent common ancestor are expected to occur more frequently, while younger mutations indicative of a recent expansion would be infrequent. Negative values of Fu’s F_S imply an excess of rare alleles at low frequency and may suggest a recent expansion (Fu 1997). The high value of Θ_S, which is positively correlated with population size, suggests a large effective population size. The estimated coalescent time agrees with previous work (García-Rodríguez et al. 2011) and concurs with a recent population expansion as indicated by the neutrality tests.

Another hypothesis suggests that high h and low Θπ may reflect “selective sweepstakes reproduction,” in which natural selection and the low probability of survival to adulthood (i.e., reproductive skew) combine to create a genetic pattern that also includes an excess of small mutations (Niwa et al. 2016; Árnason et al. 2023; Eldon and Stephan 2023). Given the biological nature of sardines, sweepstakes recruitment is expected, thus it is possible that both processes (rapid population expansion over evolutionary time scales and selective sweepstakes reproduction) have worked together to produce the current genetic architecture of Pacific sardine.

The BSP of population size shows two population expansion events, one beginning ∼280 kya and a second beginning ∼110 kya. We speculate that the ∼280 kya event corresponds to the same event estimated at ∼235 kya by the traditional coalescence approach. Population expansions in marine fishes are often considered in the context of changes in sea surface temperatures (SST) as they are reasonable indicators of suitable habitat. Around 400,000 yr BP, an increase in the amplitude of climate cycles across the globe known as the Mid-Brunhes Transition occurred. Interglaciation periods were warmer than before and resulted in warmer SST (Morley and Dworetzky 1991; Barth et al. 2018; Kender et al. 2019). Similar trends were present in the California Current, as interglacial warming periods and associated warmer sea surface temperatures have been observed by sedimentary core sample analysis (Heusser 1998; Lyle et al. 2010). Both of the expansion events identified by the BSP correspond to rapid increases in SST in the California current system following periods of SST below ∼13 °C (Herbert et al. 2001). Given that 13-16 °C is the optimal range for Pacific sardine egg and larval survival (Lasker 1964; Agostini et al. 2007; Dorval et al. 2019), we hypothesize that reproductive success during these cooler periods was reduced and prevented population growth.

Overall, the low but significant fixation indices support the null hypothesis of no genetic structure within the Pacific sardine of the northeastern Pacific that reflects evolutionary time scale processes. However, the significant values observed in the fixation between groups are worth reviewing, as significant differences between sampling groups were observed when Bahía Magdalena was grouped with the GOC (Table 4). This could imply that there is more similarity in the genetic architecture of Pacific sardine in Bahia Magdalena and the GOC. That said, spurious significance of small values of Φ_ST are common and may not be biologically meaningful (Waples 1998; Hedrick 1999; Meirmans and Hedrick 2011). Regardless of whether or not subtle structure exists between a Bahia Magdalena/GOC group and the rest of the Pacific, there is no evidence of phylogenetic structure north of Bahia Magdalena.

The molecular data presented here and the previous molecular studies on this species support the hypothesis of a single panmictic population. The data herein agree with the previous research and align with the traditional notion of substantial gene flow between pelagic fish populations. In the absence of data to the contrary, it is reasonable to accept the hypothesis that there is no genetic structure in the Pacific sardine in the northeastern Pacific Ocean. Mitochondrial studies will not reflect a recent separation as there is not enough time for lineage sorting (Avise 2000), and future genetic studies using different markers or analyses may yield different results and would expand our understanding of the genetic architecture of Pacific sardine in their northeast Pacific range.

The study was funded in part by the NOAA NMFS Genomics Strategic Initiative. We would also like to acknowledge the archival collection at the Southwest Fisheries Science Center which supplied tissue samples, numerous colleagues who facilitated port sampling in Mexico, G. Longo and the Northwest Fisheries Science Center for their help with computational analysis, and the members of the SWFSC Genetics Physiology and Aquaculture group who provided valuable assistance, comments, and guidance throughout this study.