Eastern oysters, Crassostrea virginica, are facing rapid environmental changes in the northern Gulf of Mexico and can respond to these changes via plasticity or evolution. Plastic responses can immediately buffer against environmental changes, although this buffering may impact the organism’s ability to evolve in subsequent generations. While plasticity and evolution are not mutually exclusive, the relative contribution and interaction between them remain unclear. In this study, we investigate the roles of plastic and evolved responses to low salinity in C. virginica using a common garden experiment with four populations exposed to two salinities. We use three transcriptomic analyses (edgeR, PERMANOVA, and WGCNA) combined with physiology data to identify the effect of genotype (population), environment (salinity), and genotype-by-environment interaction on both whole organism and molecular phenotypes. We demonstrate that variation in gene expression is mainly driven by population, with relatively small changes in response to salinity. In contrast, the morphology and physiology data reveal that salinity has a larger influence on oyster performance than the population of origin. All analyses lacked signatures of genotype-by-environment interaction, and in contrast to previous studies, we find no evidence for population-specific responses to low salinity. However, individuals from the highest salinity estuary displayed highly divergent gene expression from other populations, which could potentially drive population-specific responses to other stressors. Our findings suggest that C. virginica largely rely on plasticity in physiology to buffer the effects of low salinity, but that these changes in physiology do not rely on large persistent changes in gene expression.

Oyster collection and conditioning

In December 2017 and January 2018, subtidal adult C. virginica oysters were collected from two estuaries in Louisiana, USA; Calcasieu Lake (CL; 29° 50′ 58′′ N, 93° 17′ 1′′ W) and Vermilion Bay (VB; 29° 34′ 47′′ N, 92° 2′ 4′′ W). Additionally, in August 2018 oysters were collected from two estuaries in Texas, USA; Packery Channel (PC; 27° 37′ 38′′ N, 97° 13′ 59′′ W) and Aransas Bay (AB; 28° 7′ 38′′ N, 96° 59′ 8′′ W). In August 2018, all four broodstock populations were naturally induced to spawn at the Auburn University Shellfish Laboratory (AUSL) in Dauphin Island, Alabama, as described in Marshall et al. (2021). Oyster spat were maintained in upwelling nursery systems until ~6mm in shell height, and subsequently deployed in bags at an AUSL-permitted grow-out site at Bayou Sullivan, Alabama (30° 21′ 52′′ N, 88° 12′ 57′′ W) before being moved in March 2019 to the Grand Bay Oyster Park (GBOP), Alabama (30˚ 22′ 15′′ N, 88˚ 19’ 0” W) for further growth (see Marshall et al., 2021 for environmental data). In February 2020, we transferred the oysters (~18 months old) from each stock from GBOP to the static systems at Louisiana State University Agricultural Center Animal and Food Sciences laboratory building (AFL) in Baton Rouge, Louisiana.

Oysters were scrubbed, tagged, and placed into one of two 400-L tanks with recirculating artificial seawater (Crystal Sea Marinemix, Marine Enterprises International, Baltimore, MD, USA) adjusted to a salinity of 20 ppt and temperature of 20°C. Each population had ~10 oysters in each of the two tanks. The temperature was gradually adjusted at a rate of 1.0°C every 2 days until the experimental temperature of 25°C was reached. Salinity in each tank was gradually adjusted at a rate of 3 units every 2–3 days until the target salinities were reached (6ppt or 12ppt). We used artificial seawater or carbon-filtered and aerated freshwater to adjust the salinity. Oysters were fed daily at ~5% of their meat dry weight with Shellfish Diet 1800® (Reed Mariculture Inc, Campbell, CA, USA).

After 3 weeks of acclimation at each salinity, we began measuring clearance rates and oxygen consumption rates. We had limited access to the laboratory due to COVID restrictions, which increased the time oysters spent in these tanks and increased overall mortality. This increased mortality led to an uneven sampling between the populations and salinities (Table S1). In July 2020 we measured the final shell height, shell width, shell weight, and dry meat weight of the surviving oysters for each population at each salinity (Table S1). Shell height was measured from shell umbo to distal edge and shell width was measured as the maximum distance between the two valves when closed (Galtsoff, 1964). Dry meat weight was measured after drying at 70 °C for 48 h. For gene expression analysis, we sampled 5 oysters per population per salinity (n=40). Approximately 0.5 cm² piece of gill tissue was dissected and immediately frozen in liquid nitrogen and placed in a -80˚C freezer until later extraction.

Measuring clearance rates

The clearance rate, defined as the volume of water completely cleared of suspended particles per unit of time, was quantified using a static system (Widdows, 1985). Oysters were individually placed in 1 L beakers filled with 1 L of 0.5-μm filtered seawater and left undisturbed for ~ 60 min. Shellfish Diet 1800® was added to each beaker to give an initial concentration of 3 × 10⁴ cells mL⁻¹. Particle concentrations (>5 μm) were collected every 30 seconds until particle count declined to 50% of original values and was measured using a PAMAS Model S4031GO particle counter (PAMAS Partikelmess-und Analysesysteme GMBH, Rutesheim, Germany). We used air stones to reduce algal sedimentation. Beakers with algal cells and oyster shells inside were used as controls.

Individual oyster clearance rate (CRi) was calculated according to the equation CRi = [(b – b') × Vol(L) × 60 min h⁻¹], where b is the slope of the linear regression between the natural logarithm of cell concentration (cells mL⁻¹) and time (min) in a beaker, b′ is the slope for the control beaker, and Vol is the volume of seawater in the beaker (1 L). This method follows Cranford et al. (2016), in which a generalization of Coughlan’s (1969) method is developed to integrate the values of intermediate particle concentration samplings, which aims to reduce the impact of outliers. Accordingly, only linear regressions with an r² > 0.90 were included in these calculations. CRi measurements are in units of L h⁻¹.

Clearance rates were also standardized by shell height, specifically to a standard oyster of 80 mm (i.e., average shell height for the oysters) using the equation CRh =(?_std/?_exp)^?×CRi, where CRh is the clearance rate standardized by shell height, Hstd is the shell height of the standard oyster (80 mm), Hexp is the shell height of the experimental oyster, CRi is the individual oyster clearance rate, and b is the allometric exponent for shell height set to 1.78 (Cranford et al. 2011). CRh measurements are in units of L h⁻¹ 80 mm⁻¹.

Measuring oxygen consumption rates

To measure the oxygen consumption rate of fed oysters, known as the routine oxygen consumption rate (Bayne, 2017; Gosling, 2015), individual oysters were placed in 915-mL acrylic chambers filled with 0.5-μm filtered seawater from their respective tank. We used a ProOBOD probe (YSI Incorporated, Yellow Springs, OH, USA) to measure dissolved oxygen and temperature recorded every 30 s by a YSI Multilab 4010-3 m (YSI Inc.) connected to the probe. The readings started once the oyster’s valves opened and ran continuously for up to 90 min or until dissolved oxygen fell below 70% saturation.

Individual oyster oxygen consumption rate (OCRi) was calculated as OCRi = [(b – b') × Vol (L) × 60 min h⁻¹], where b is the slope of the linear regression of oxygen concentration (mg L−1) versus time, b′ is the slope for the control, and Vol is the volume of water (in L) in the chamber (chamber minus oyster volume). OCRi measurements are in units of mg O₂ h⁻¹. Only linear regressions with an r² > 0.95 were included in the calculations. Oxygen consumption rates were also standardized by dry meat weight (OCRw), specifically to a standard oyster of 1 g dry meat weight. We used the equation OCRw =(?_std/?_exp)? ×OCRi, where W_std is the dry meat weight of the standard oyster (1 g), W_exp is the dry meat weight of the experimental oyster, and b is the allometric exponent, 0.58 (Casas et al., 2018). OCRw measurements are in units of mg O₂ h⁻¹ g⁻¹.

Statistical analyses

We examined the effect of population, salinity, and population x salinity interaction on clearance rates (CRh), oxygen consumption rates (OCRw), and dry meat weight. We used a Shapiro-Wilk normality test (Shapiro & Wilk, 1965; Royston, 1982) to determine if the data followed a normal distribution and used either a log or square root transformation if needed. A two-way ANOVA was performed using the R function “aov” in the stats package (version 3.6.2) (p-value < 0.05). To examine pairwise multiple comparisons we used the post-hoc Dunn’s test with the Benjamini-Hochberg correction (p-value < 0.05) (Benjamini & Hochberg, 1995).

Gene expression analysis; sampling and initial processing

Total RNA was extracted using an E.Z.N.A.® Total RNA Kit I (Omega BIO-TEK Inc., Norcross, GA, USA; VWR) following the manufacturer's instructions. The yield and quantity were initially assessed using a NanoDrop 2000 spectrophotometer. Total RNA extracted from the 40 individuals was sent to the University of Texas at Austin’s Genomic Sequencing and Analysis Facility where RNA quality was confirmed using a 2100 Agilent Bioanalyzer on a Eukaryote Total RNA Nano chip, and libraries were produced using the Tag-Sequencing approach (Meyer et al., 2011). The resulting 40 libraries were pooled at equimolar ratios and sequenced across two lanes of an Illumina NovaSeq platform, with 100 base pair single-end reads.

Sequencing reads were trimmed of adapter sequences and base pairs with quality scores below 20 were removed using Trimmomatic (version 0.39) (Bolger et al., 2014). The trimmed reads were mapped to the C. virginica reference genome (Gómez-Chiarri et al., 2015) with known haplotigs removed (https://github.com/jpuritz/OysterGenomeProject/tree/master/Haplotig_Masked_Genome) using the single pass option for STAR RNA-seq aligner (version 2.6.0a) (Dobin et al., 2013). Three samples displayed poor sequencing quality (i.e., the lowest number of reads, poorest mapping rates, and isolated in PCA) and were subsequently removed for all downstream analyses. A count matrix was generated from the ReadsPerGene.out.tab output from STAR. To ensure statistical robustness, we used three different approaches to examine gene expression: edgeR, a PERMANOVA test, and a Weighted Gene Co-Expression Network Analysis (WGCNA).

Gene expression analysis; edgeR

We first assessed pairwise changes in gene expression associated with population and salinity using the package edgeR (version 3.5.2) (Robinson et al., 2010). We filtered all genes using the command (filterByExpr) and the remaining read counts were normalized using a trimmed mean of M-values (TMM) normalization method (Robinson & Oshlack, 2010). Broad changes in gene expression were first assessed using a principal component analysis (PCA) conducted using the R program vegan with euclidean distances calculated from log2 +1 transformed normalized counts obtained from the cpm() function in edgeR. Pairwise differential expression was measured using the genewise negative binomial generalized linear model implemented in the edgeR function glmQLFit, and significantly differentially expressed genes (DEGs) were identified based on False Discovery Rates (FDR) calculated using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995). Functional enrichment of differentially expressed genes was tested using a Fisher’s Exact Test (p-value < 0.05). All Fisher’s Exact Tests were conducted using the scripts originally developed by Wright et al., 2015 and available at: (https://github.com/z0on/GO_MWU/blob/master/GO_MWU.R). This method uses a binary input (DEGs = 1, non-DEGs = 0) to calculate if there is enrichment in the DEGs across three gene ontology (GO) categories: Molecular Function (MF), Biological Processes (BP), and Cellular Component (CC).

Gene expression analysis; PERMANOVA

We used a permutational multivariate analysis of variance test (PERMANOVA) to identify genes associated with population, salinity, or population x salinity interaction. We performed two PERMANOVA analyses to examine the effect on total gene expression and on a per gene basis. For both PERMANOVA analyses, log-transformed counts were evaluated using the R function ‘adonis2’ within the ‘vegan’ package (version 2.5-6). The PERMANOVA examined the effect of gene expression~ population + salinity + population*salinity with 10⁶ permutations (total expression) or 99,999 permutations (per-gene expression) (p-value < 0.05). The difference in the number of permutations between the total expression and per-gene analysis is the result of the higher computational requirements for the per-gene analysis. For the per-gene PERMANOVA, the p-values were corrected for multiple comparisons using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995). Functional enrichment of significant genes was tested using a Fisher’s Exact Test (p-value < 0.05).

Gene expression analysis; WGCNA

We used a Weighted Gene Co-Expression Network Analysis (WGCNA) (Langfelder & Horvath, 2008) to find clusters of genes with highly correlated expression patterns associated with population, salinity, oxygen consumption rate, final shell height, shell width, and dry meat weight (n = 37). We also ran a separate WGCNA analysis for individuals with clearance rate measured (n=23). For all WGCNA analyses, we restricted the number of genes to remove lowly expressed features by only using samples with greater than 5 counts per million in 80% of all samples. WGCNA was run using the 7,575 genes (7,533 genes for clearance rate) that passed this filter, a soft-threshold of 16, and a minimum module size of 30. Functional enrichment of significant modules was tested using a Fisher’s Exact Test (p-value < 0.05).