Population genomics of eastern oysters, Crassostrea virginica, in a well-mixed estuarine system: advancement and implications for restoration strategies

Abstract

Eastern oysters, Crassostrea virginica, are historically a keystone species in many of the estuaries in which they reside, providing critical ecosystem services. Because oyster populations have been on the decline for the past century, restoration initiatives currently are underway in many estuarine systems, including Great Bay Estuary (GBE), New Hampshire. Results of prior studies of eastern oyster population genomics cannot be applied directly to GBE, as it is a well-mixed estuarine system that is relatively contained, and the sources of recruits are split among cultivated and native. This study aimed to identify the population genomic structure of eastern oysters in GBE to facilitate determination of effective population size and estimation of genetic differentiation among subpopulations. Results showed moderate genomic differentiation among native, cultivated, and restoration C. virginica subpopulations in the Bay. A small number of breeders (N_e=163–276) was found in all subpopulations except the Lamprey River site (N_e=995). This research provides a contemporary snapshot of eastern oyster subpopulation structure at the genomic level in GBE that will facilitate restoration and enhanced management.

Background

Eastern oysters, Crassostrea virginica, are historically a keystone species in many of the estuaries in which they reside, with an estimated economic value of $5,500–$99,000 ha^–1 y^–1 [1]. They provide numerous functions that are critical to their ecosystems [2] including water filtration [2, 3], nitrogen conversion [4,5,6], act as a sanctuary for juvenile fishes in the form of oyster reefs [7], and reduce wave energy from storm surge to protect shorelines [8, 9]. A cause for concern is that eastern oysters have been on the decline throughout the entirety of their range for several decades, decreasing in oyster numbers by 64% and oyster biomass by 88% between the early 1900s and the early 2000s [2]. This could prove to be catastrophic for some estuarine ecosystems [2, 10, 11], particularly estuaries in New England where climate change is occurring at an alarming rate. The Gulf of Maine, where Great Bay Estuary (GBE), New Hampshire is located, currently is one of the fastest warming regions on Earth, where sea surface temperatures increased at a rate of 0.26 °C y^− 1 between 2004 and 2012 [12]. Eastern oyster populations in GBE presently are only 10% of what they were in the 1980s, with the most common causes of decline being attributed to disease (specifically MSX (multinucleated sphere with unknown affinity X) and Dermo), over harvesting, lack of substrate, and sedimentation leading to shell burial [13]. In GBE, MSX and Dermo are consistently present in eastern oysters; in 2018, MSX infections were found in roughly 10% of the eastern oyster population and Dermo in roughly 80% [14]. Studies that increase our knowledge of eastern oysters and that could lead to increased success of restoration strategies are of great importance.

Eastern oyster restoration projects in GBE currently are underway through collaborative efforts of multiple agencies, such as the University of New Hampshire, The Nature Conservancy, New Hampshire Sea Grant, and New Hampshire Fish and Game. Restoration efforts from 2000 to 2018 included deployment of shell, live spat-on-shell, and transplantation of large, reproductively active cultivated oysters (deemed “uglies”) onto pre-selected restoration sites [15]. Following a large spat-on-shell restoration initiative by The Nature Conservancy in 2011, it was observed that approximately 5.8 × 10⁴ oyster spat recruited to a 1-ha reef located at the Lamprey River [16], which subsequently showed indication of natural recruitment [13]. In 2013, a survey of multiple GBE oyster restoration sites indicated some levels of success (i.e., recruitment), but found numerous reef areas where few, if any, live oysters were found post-restoration [15,16,17]. In a 2015 project, only 3 of 9 sites had numbers of live oysters adequate to constitute a healthy, living oyster reef [13, 15]. A 2019 assessment provided similar results [15]. A 2021 project by The Nature Conservancy, titled “The Purchase Program” under their Supporting Oyster Aquaculture and Restoration initiative, deployed large aquacultured “uglies” at restoration sites (the success of this project is not yet known).

Most studies have emphasized the potential for larval recruitment of oysters in GBE, suggesting that recruitment may be efficient at sites near existing natural reefs with an established adult oyster population [18]. The historic decrease in the eastern oyster adult population in GBE may limit the potential for recruitment from native reefs [13]. Although the approximate densities and locations of oysters in GBE are known [19], questions have arisen relating to the number of breeders, gene flow, and the potential for recruitment to both natural and artificial oyster reefs [20]. Because restoration has been challenging in GBE, detailed information on the population structure of eastern oysters in this estuary is needed to enhance restoration efforts to allow for discovery of which oysters should be focused on for restoration (i.e., the populations with high effective breeding populations) and to determine the influence of native and cultivated oysters on the restoration populations.

A contemporary approach to examining population structure, population genomics, interrogates numerous regions of the genome to better understand genetic variation among populations [21]. This approach allows estimation of migration, population differentiation, effective number of breeders, random mating, etc. Historically, management decisions often have been made without consulting important genetic information [22], particularly with respect to marine invertebrates such as the eastern oyster. More recently, genomic studies are being introduced into marine management decisions [23,24,25,26]. Such studies (including this one) encompass several important parameters: inbreeding coefficients (F_IS), population subdivision (F_ST), effective population size (N_e), heterozygosity, and population clustering. More often than not, in many restoration schemes the issue of preserving genetic diversity arises [27] as hatchery-reared individuals that perform better are inadvertently selected without considering the inbreeding consequences of non-random mating. While this type of inadvertent selection can be addressed by strip-spawning in the hatchery setting, reproductive variance in the hatchery can lead to loss of genetic diversity and lower effective population sizes [28]. Recent studies of marine oyster populations indicate strong population structure on a fine-scale [29, 30], alluding to the value of understanding these fine-scale population characteristics in estuarine regions to develop more appropriate restoration and management strategies. Bivalve restoration practices, including eastern oysters, commonly include supplementing native populations with hatchery-reared individuals in hopes of increasing total reproductive output, higher recruitment, and ultimately successful reef restoration [13, 16, 31]. This may not be the best strategy if cultivated transplants have lower N_e, in which case alternative best practices for restoration may be necessary. This study aimed to address such scenarios.

This study used single nucleotide polymorphisms (SNPs) to answer the necessary questions regarding eastern oyster population genomics in a well-mixed estuarine system. Great Bay Estuary, NH has a unique set of characteristics that separate it from other estuarine systems. Generally small freshwater fluxes (2% of tidal prism) [32] and strong tidal mixing result in negligible stratification (except very close to the river mouths) within the system [33]. Therefore, it is difficult to extend the knowledge from previous eastern oyster population genomics studies in other geographic areas [29, 34, 35] to GBE, as well as other studies from different oyster species, such as Crassostrea gigas [56] and Ostrea lurida [30]. To provide population genomic information to inform GBE restoration practices, low-coverage, whole-genome sequencing (lcWGS) was performed on 210 individuals from 7 different sites in GBE with approximately 5× coverage. The lcWGS approach is a powerful low cost tool that has been shown to be effective for describing population structure [36, 37]. Despite previous concerns about the reduction in depth of coverage and confidence in individual genotype assignment, recent software allows for lcWGS-specific analyses for population genomics. Additionally, researchers obtain greater breadth of coverage and larger sample sizes when conducting lcWGS for the same cost [36]. This study aimed to (1) assess the population genomic structure of cultivated, native, and restoration subpopulations of eastern oysters in GBE, (2) determine levels of genetic differentiation among subpopulations, and (3) estimate the effective breeding size of each subpopulation sampled. Due to the unique characteristics of GBE, especially its well-mixed nature and strong tidal currents [33], it was expected that genomic structure of native oyster groups would be homogenous, that the restoration specimens would show signs of either native or cultivated heritage, whichever group was most successful at that site, and that major differentiation between native and cultivated populations would be evident due to unique heritage and possible reproductive variance effects.

Methods

Sample collection

Eastern oysters (n = 210) were collected from seven sites throughout Great Bay Estuary, NH, USA (30 individuals per site) (Fig. 1). Native reefs in Lamprey River (LR), Squamscott River (SQ), Oyster River (OR), and Adam’s Point (AP) were sampled. Restoration oysters were collected from an established restoration site at Nannie Island (NI). Cultivated oysters were collected from two oyster farms near Fox Point (FP) and Cedar Point (CP) (Table 1). The MSX-resistant strain was derived from the Haskin Northeast High Survival (NEH) Diploid Oyster (obtained from Muscongus Bay Aquaculture, Inc.). Mantle tissue was collected from each individual and cleaned by gentle agitation for 1 min in 5 vol of 99% ethanol, followed by a 3 min soak in 5 vol of 3% sodium hypochlorite, and rinsed by soaking 1 min in 5 vol of 99% ethanol. Cleaned mantle tissue samples were preserved with 70% ethanol and stored at -20 °C until DNA extraction.

Fig. 1

Map of Great Bay Estuary in New Hampshire, USA. Native oyster reefs at Lamprey River (LR), Squamscott River (SQ), Oyster River (OR), and Adam’s Point (AP) are marked in shades of blue, the restoration site at Nannie Island (NI) is marked in yellow, and farms Fox Point (FP) and Cedar Point (CP) are marked in shades of red

Table 1 Sampling locations of Crassostrea virginica populations in Great Bay Estuary, NH, USA

DNA extraction and preparation

Genomic DNA of C. virginica was extracted from mantle tissue using the DNEasy PowerSoil Pro Kit following the manufacturer’s protocol (Qiagen, Hilden, Germany) and quantified using the Qubit 1× dsDNA High Sensitivity Assay Kit (Life Technologies, Foster City, CA, USA). Library construction and sequencing were conducted at University of New Hampshire’s Hubbard Center for Genome Studies using the Kapa BioSystems HyperPlus Kit (KR1145 -v3.16) and NovaSeq 6000 with an SP flow cell (paired-end 250 bp reads). Data were demultiplexed using bcl2fastq v2.20.0.422.

Sequence assembly, filtering, & SNP identification

The NovaSeq SP PE produced 250 base-pair paired-end (2 × 250 bp) reads in FASTQ format. Illumina adapters and low-quality bases (Q ≤ 20) were trimmed using Trimmomatic version 0.40 [38]. Quality trimmed paired-end reads were mapped to the C. virginica genome (NCBI assembly GCA_002022765.4) using Burrows-Wheeler Aligner (BWA) [39]. Mapped reads were sorted and indexed using SAMtools [40] and duplicate reads were marked using the Genome Analysis ToolKit (GATK) [41]. Following these steps, variant calling was performed wherein single nucleotide polymorphisms (SNPs) and insertions and deletions (indels) were called using the FreeBayes haplotype-based variant caller version 1.2.0 [42]. The output created by FreeBayes was a single file with all variants for all samples in VCF file format.

To enable downstream analysis of SNPs, the VCF file containing all variants was filtered using a series of steps in VCFtools version 4.2 [43]. First, all variants except SNPs were removed from the VCF file. Once only SNPs remained, the SNPs were filtered based on the following parameters: present in at least 50% of individuals, minimum quality ≥ 30, and minor allele count (MAC) ≥ 3. Following filtering, individuals that were missing more than 50% of their data were removed from further analysis (n = 42). This resulted in 18 individuals from site AP, 28 from CP, 30 from FP, 27 from LR, 24 from NI, 24 from OR, and 18 from SQ for downstream analyses. Thereafter, SNPs were filtered so that each was present in ≥ 85% of individuals, and each had a minor allele frequency of ≥ 5%. The SNP distributions across the C. virginica genome were plotted using a Manhattan plot and associated p-values were generated using PLINK [44].

Population genomic analyses

Analyses followed the dDocent pipeline [45]. Population genomic analyses and plots, excluding ADMIXTURE and N_e, were conducted in RStudio version 2023.9.1.494 [46]. Pairwise F_ST values were estimated to discover the amount of genetic variance that could be explained by population structure [47], calculated using the package hierfstat [48]. Expected (H_e) and observed (H_o) heterozygosity and inbreeding coefficients (F_IS) were calculated using the packages vcfR [49] and heirfstat, respectively. The genetically effective size of each subpopulation (N_e) was calculated using currentNe, which is based on the linkage disequilibrium method, is beneficial for small sample sizes, and provides more consistent confidence intervals [50]. Due to not meeting the assumptions required for estimating effective migration rate, N_em, this parameter was excluded from analysis [51], and F_ST was used as the sole parameter to estimate subpopulation differentiation. A discriminant analysis of principal components (DAPC) was conducted for all populations to determine population clustering based on SNP profiles with the package adegenet [52]. ADMIXTURE analyses [53], which disregarded known classification of each individual and used instead their SNP profiles to assign individuals into clusters, were completed for each individual and K-values were cross-validated to obtain the optimal number of clusters (K). Once completed, individual assignments based on the SNP profile were plotted in a bar chart as well as overlaid on a map of GBE [54].

Results

A total of 26,275,355 variants were identified using FreeBayes from the Crassostrea virginica low-coverage whole-genome sequencing. Following removal of non-SNP variants, filtering for quality, and removal of individuals with missing data, a complete SNP profile of 6,657 SNPs for 168 individuals remained for downstream analyses (Supplementary Table 1). SNPs were uniformly distributed throughout the C. virginica genome (Fig. 2). Excessive missing data (> 50%) in some individuals was caused by too low amounts of template DNA. Pairwise F_ST values (Table 2) were low for all comparisons among all sites (range 0.002–0.042). Among the native populations, the restoration site (NI) accounted for the highest F_ST values (0.030 < F_ST<0.042). No clear trends among native, restoration, and cultivated sites were observed for heterozygosity; the expected values for heterozygosity were similar among all sites (0.207<H_e<0.230), as were the observed values for heterozygosity (0.204<H_o<0.226) (Table 3). The lowest observed heterozygosity (H_o = 0.204) was in FP, and NI had the highest (H_o = 0.226).

A majority of SNPs (> 95%) conformed to the expectations of Hardy-Weinberg Equilibrium and all SNPs were retained in this study, even those that failed HWE expectations, as it has been shown in other species that removing them can cause a loss of relevant information during analysis [55]. Additionally, when SNPs were considered at the site level [35, 45], all SNPs conformed to HWE. All sites had negative F_IS values, ranging from − 0.201 at NI to -0.157 at SQ (Table 3). Effective population sizes (N_e) were variable across all sites ranging from 163 at SQ to 995 at LR (Table 3). Oyster subpopulations at cultivated sites had N_e of 216 at FP and 261 at CP. The restoration site NI had N_e of 224.

Fig. 2

Manhattan plot showing the distribution of SNPs throughout the genome of Crassostrea virginica and their associated p-values

Table 2 Pairwise F_ST values estimated using SNP frequency data for native, restoration, and cultivated Crassostrea virginica subpopulations in Great Bay Estuary. LR: Lamprey River, SQ: Squamscott River, OR: Oyster River, AP: Adam’s point, NI: Nannie Island, FP: Fox Point, CP: Cedar Point

Table 3 Descriptive statistics for Crassostrea virginica sample sites. Effective population size (N_e) and 90% confidence intervals (CIs), observed and expected heterozygosity (H_o and H_e, respectively), and inbreeding coefficient (F_IS) were calculated. LR: Lamprey River, SQ: Squamscott River, OR: Oyster River, AP: Adam’s Point, NI: Nannie Island, FP: Fox Point, CP: Cedar Point

Four distinct clusters of individuals were present in the DAPC: one with all native populations except the intertidal site (site AP) (shades of blue), one containing cultivated individuals (shades of red), the restoration population (yellow), and the intertidal population (light blue) (Fig. 3). Ellipses at 95% confidence encompassed a majority of individuals within each cluster. Cross-validation via ADMIXTURE analysis demonstrated that k = 3 was the optimal value for clustering this genetic dataset. Individual assignments (Fig. 4A) showed 3 distinct clusters, with influence from all clusters on each other. Native populations (LR, SQ, OR, AP) fell within cluster 2 (blue), cultivated populations (CP and FP) fell within cluster 1 (red), and the restoration population (NI) clustered on its own in cluster 3 (yellow). When overlaid on a map (Fig. 4B), location did not appear to affect an individual’s genetic profile, but rather their heritage (native, cultivated, restoration) was the salient feature.

Fig. 3

Discriminant analysis of principal components (DAPC) showing differentiation of SNP genotypes among seven eastern oyster subpopulations in Great Bay Estuary, NH. Ellipses represent 95% confidence intervals, and each dot represents an individual’s SNP profile. Individuals collected from native subpopulations are shown in shades of blue, cultivated subpopulations are shown in shades of red, and the restoration site is shown in yellow

Fig. 4

ADMIXTURE analysis based on SNP profiles of Crassostrea virginica individuals sampled in Great Bay Estuary. Panel A) Individual cluster assignments from ADMIXTURE analysis. B) Summary of individual assignments in each population obtained from ADMIXTURE analysis overlaid on a map of Great Bay Estuary, NH, USA

Discussion

Assessing population structure and genetic variance for declining population restoration programs is critical to maximize the success and benefits of these programs; however, this has not been done for many marine populations. Preserving high genetic variance in restorations is a major goal as this enhances long-term population viability. For the Crassostrea virginica population in the well-mixed GBE, this is the first study of genetic structure. The pipeline used as a guide [45] yielded a total number of SNPs similar to other studies of oysters [30, 35, 56]. Although the C. virginica genome has a high number of repetitive regions [57], the SNPs found in this study did not excessively cluster with known repeat regions (Fig. 2). There are several key takeaways from the SNP analyses that will enhance restoration activities. First, this study demonstrated that eastern oyster populations in this system cluster together based on heritage rather than geographic location. This knowledge will facilitate determination of the source of recruits to the native population as suggested by other researchers attempting to assess relative contribution of sources [58]. Second, very low population differentiation was present, indicating that despite the dramatic reduction in oyster numbers, geneflow remains substantial among the various sites. Third, the genetically effective population sizes (N_e) were very low for all samples except LR, which is consistent with other estimates of N_e from small local C. gigas populations [59]. The effective size was considerably less than the estimated census sizes (N_c) for most sites, as expected. For example, the effective number of breeders at SQ (N_e = 163, the lowest encountered in this study) was 0.1% of the estimated census size at SQ [15]. This difference between N_e and N_c is consistent with other marine studies and is in fact expected as a consequence of unequal sex ratios (small oysters are predominantly male), reproductive variance (oysters mass spawn, so the relative locations can have a large effect on the numbers of offspring attributed to each adult), population fluctuations, and reduced population size (certainly a factor for GBE oysters due to the extended period of harvest and disease-related mortality) [21, 60]. The LR oysters, with a high estimated N_e, presently is the healthiest group of oysters in GBE. The remaining native populations with N_e of approximately 200, will retain ~ 99.6% of genetic variation each generation, which should maintain population viability over the short-term allowing continued focus on conserving habitat. But if/when those populations continue to lose breeding adults and ultimately begin to approach N_e of 100, genetic variance decay will accelerate, inbreeding will increase, and the ability of the GBE oyster population to adapt to climate change, disease, etc. will be at risk [61].

Together, these results paint a picture of genomic exchange within a well-mixed estuarine system that differs from previous studies conducted in other locations such as Chesapeake and Delaware Bays. In those estuarine systems, there is strong evidence of genomic geographic differentiation and isolation by distance [35, 62, 63], whereas in the present study, geographic location had little influence on genomic profiles. Very low levels of genetic differentiation were observed among the native oyster subpopulations (SQ, LR, OR, AP) as was expected given the small total area of GBE and its well-mixed character. No relevant genetic differentiation was observed between the cultivated subpopulations (FP, CP). The values among native GBE populations were considerably lower than F_ST estimates for eastern oyster population genomic studies in Gulf of Mexico [34], Delaware Bay [64], and Chesapeake Bay [35, 62,63,64]. Low levels of differentiation also were observed between the two cultivated populations and the four native populations, which was not expected given the selective breeding heritage of the cultivated oysters in this study. These low values are not a consequence of repeat regions because the SNPs were not predominantly associated with repeats. Although the F_ST values indicate small differences among the groups, the direct comparison of the native versus cultivated samples (F_ST=0.025) was 4× the genetic divergences among samples within those two groups; an observation similar to the divergence measured between native and selectively bred C. gigas [59]. Even the restoration site (NI) showed low to very low levels of genetic differentiation from both the native and cultivated groups, respectively. Since the NI restoration site is adjacent to a known native reef, and cultivated oysters of the same strain studied here have on multiple occasions over the past two decades been deployed on the restoration site, it was predicted that the NI sample would show genetic similarity to either the native or cultivated populations, whichever group was the most represented in recruits to that population. The results of F_ST, DAPC, and ADMIXTURE combined imply that the cultivated strain has had slightly more of an effect on the current NI population than the native oysters. This is an important finding that informs restoration practices in this area by confirming that recruitment to the artificial reef is in part due to reproduction of cultivated strains, which may have better survival and success due to their heritage from selected strains.

Genetically effective population size (N_e) is a critical tool to estimate evolutionary history and the potential for loss of genetic variability, especially in a species with declining populations [65, 66]. Because it is affected by iteroparity, reproductive success, age at maturity, lifespan, and other life-history traits, estimates of N_e generally are much smaller than the census size of a population. This is particularly important in marine populations [67, 68], implying that only a small portion of individuals in marine populations act as breeders. With the exception of the site at Lamprey River, N_e estimations in this study (Table 3) showed relatively low effective population sizes, as expected, reflecting the declining population size in GBE and the limited numbers of breeding individuals for the cultivated populations. These values of N_e for wild oysters, ranging from 163 to 995, are similar to those found in other estuaries along the eastern United States, such as Delaware Bay [69] and Chesapeake Bay [35]. These N_e values can inform both planning of conservation efforts and assessment of their success and impact.

Both clustering methods illustrated that location within the western rim of this estuary system did not play an appreciable effect on genetic profile, as native reefs spanning the entire GBE all were genetically similar (Fig. 4A). While results of DAPC showed divergence of the intertidal oysters at Adam’s Point from the subtidal oysters, the ADMIXTURE and F_ST results showed this divergence was minimal. Rather than isolation by distance playing a role in the population genomic structure of eastern oysters in GBE, heritage (native, cultivated, restoration) had a stronger influence on population clustering and assignment (Figs. 3 and 4), an observation that is contrary to studies examining population structure in Chesapeake Bay [35]. In GBE, population structure by heritage rather than distance is likely due to the well-mixed nature of this unique estuarine system. Studies in the future should continue to examine the population genomics of eastern oysters in this estuary system. As oyster farming has grown exponentially in Great Bay Estuary [70], the effects of a well-mixed estuary system might not yet be seen on the oysters’ genomic profiles. In the future, studies may show that there are no unique genomic profiles among native, cultivated, and restoration populations. Therefore, these studies should continue to occur at a regular pace.

The overall unique population stratification among native, cultivated, and restoration oyster subpopulations within GBE can inform managers and help them to make more precise decisions for eastern oyster restoration. Current restoration strategies deploy large, cultivated oysters to a central site within GBE, a native reef that once was highly productive but has since lost almost the entirety of its population. One major concern with this strategy has been ensuring the genomic diversity of each subpopulation is sustained to protect against threats such as disturbances associated with a changing climate. The current genomic results demonstrate a sustained level of genomic diversity with low levels of isolation by distance. Restoration efforts in GBE should keep in mind the importance of sustaining the current level of diversity to ensure best chances of survival of eastern oyster groups.

Conclusion

Discovering and understanding the population genomics of cultivated, native, and restoration subpopulations of Crassostrea virginica in GBE is critical to streamlining restoration practices and ensuring management strategies are using best practices. This research, in conjunction with knowledge of the local ecosystem, provides a more complete picture of eastern oyster population structure in GBE that is relevant to oyster restoration. This study showed low to very low genomic differentiation among native, cultivated, and restoration C. virginica subpopulations. As expected, native and cultivated subpopulations exhibited differences evident of their unique heritage. Because the restoration site is known to be influenced by both the native and cultivated subpopulations (formerly a natural reef, proximal to an existing productive native reef, and having been supplemented with cultivated oysters), it was expected that the oysters at the restoration site would exhibit genomic characteristics of native, cultivated, or a mixture of both subpopulations. Interestingly, the genomic analysis shows not a distinct signal of either ancestral subpopulation but instead a strikingly homogenous unique signal. Restoration strategies can utilize these genomic results to enhance selection of restoration sites (i.e., Which subpopulations are in most need of additional breeding individuals? How many individuals are needed per area?). The data also help to identify areas that have the greatest conservation value (e.g., Lamprey River has a large breeding population).