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 remains unclear. In this study, we investigated the roles of plastic and evolved responses of C. virginica acclimated to low salinity using a common garden experiment with four populations exposed to two salinities. We used three transcriptomic analyses (edgeR, PERMANOVA and WGCNA) combined with physiology data to identify the effect of genotype (population), environment (salinity) and the genotype–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 the genotype×environment interaction and, in contrast to previous studies, we found no evidence for population-specific responses to low salinity. However, individuals from the highest salinity estuary displayed highly divergent gene expression from that of 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.

Climate change is causing rapid alterations to many of the world's most productive coastal ecosystems. These changes include alterations in temperature, pH, salinity and dissolved oxygen (Doney et al., 2012; Breitburg et al., 2015; Gunderson et al., 2016; Schmidtko et al., 2017). There is a need to understand how keystone species will respond to these changes to better predict which species or populations are under threat from future anthropogenic change. Plasticity and evolution are two processes that enable organisms to respond to environmental change (Parmesan, 2006; Hoffmann and Sgró, 2011; Seebacher et al., 2015). Plasticity is a single-generation process that refers to the variation in phenotypes produced from a single genotype (West-Eberhard, 2003; DeWitt and Scheiner, 2004). Although plasticity is not expected to generate additional changes in successive generations beyond the scope attained in the initial response, plastic responses can act as an immediate buffer against environmental changes (Seebacher et al., 2015). In contrast, evolution occurs over multiple generations and describes the process in which allele frequencies change over time because of differential reproductive and survival success among genotypes. Although these two process are unique, plasticity and evolution are not mutually exclusive; plasticity itself can evolve (Via and Lande, 1985). However, the relative contributions of and interactions between plasticity and evolution in populations responding to environmental change remain unclear (Price et al., 2003; De Jong, 2005; Ghalambor et al., 2007; Wund, 2012; Forsman, 2015; Hendry, 2016).

Our study focused on the eastern oyster, Crassostrea virginica, which is a foundation species that spans the majority of the northern Gulf of Mexico (nGOM) estuaries (Nevins et al., 2014; Volety et al., 2014; La Peyre et al., 2019). In nGOM estuaries, salinity is one of the main environmental variables affecting the distribution and performance of this species (Parker, 1960; La Peyre et al., 2009; Rybovich et al., 2016). Crassostrea virginica can tolerate extreme salinity (0–45 parts per thousand, ppt) over short periods of time, but populations are best sustained at intermediate salinity (4–35 ppt) (Lavaud et al., 2021; Swam et al., 2022). The main concern is that local salinity regimes are expected to change rapidly as a result of anthropogenic alterations to coastal hydrology and from changes in rainfall patterns (Bishop et al., 2019; Cayan et al., 2010; Das et al., 2012). The broad distribution of C. virginica along this major salinity gradient makes this system ideal for using a space-for-time substitution (Pickett, 1989) to clarify the roles of plastic and evolved responses to environmental change.

Local adaptation describes a particular type of evolved response in which populations become optimized to their local conditions and differ in their stress tolerance, such that the response of any given population will not be representative of the species across its range (Sanford and Kelly, 2011). Previous research took advantage of C. virginica’s broad distribution along the nGOM salinity gradient and examined individuals from four populations originating from distinct salinity regimes (Marshall et al., 2021a,b). Marshall et al. (2021a,b) investigated the response to salinity in field and lab settings and found evidence for local adaptation, with the population originating from the highest salinity estuary (Packery Channel) having significantly higher mortality and smaller size under low salinity conditions.

Variation in stress tolerance among populations provides an opportunity to investigate mechanisms of low salinity tolerance in C. virginica. The present study focused on physiological and transcriptomic responses to low salinity across the four populations studied in Marshall et al. (2021a,b) to determine the relative contribution of plastic and evolved responses, as well as to test for signatures of local adaptation in different traits. We incorporated comparative transcriptomics because it is a particularly useful approach for quantifying the relative contributions of plasticity and evolution to phenotypic responses to environmental change (DeBiasse and Kelly, 2016). This is because the expression level of every transcript is a phenotype resulting from the combination of the genotype (G, genetic variation in expression), environment (E, plasticity, the environmental effect on expression) and the genotype×environment interaction (G×E, genetic variation in plasticity) (Rockman, 2008; Levine et al., 2011; Zhou et al., 2012; Koch and Guillaume, 2020; Sirovy et al., 2021).

To investigate the roles of plastic and evolved responses to low salinity in C. virginica, we used a common garden experiment with four populations exposed to two different salinity regimes. We aimed to quantify the contributions of plasticity, genetic variation and genetic variation in plasticity on C. virginica’s new steady state after acclimation to low salinity. To accomplish this, we combined measurements of morphology (i.e. final shell height, shell width, shell mass, dry meat mass), physiology (i.e. clearance rate, oxygen consumption rate) and 3′ RNA sequencing (TagSeq) to identify the effect of genotype, environment and the genotype×environment interaction on the oyster's response to low salinity. Based on the results from Marshall et al. (2021a,b), we expected that if local adaptation is occurring then the oysters from the highest salinity estuary (Packery Channel) should show a greater drop in performance (i.e. lower clearance rate, higher oxygen consumption rate) under low salinity. Additionally, if adaptation to low salinity is associated with adaptive changes in gene expression in response to low salinity events, we expected either less change in gene expression in oysters from Packery Channel under low salinity or a gene expression signature that is dominated by stress response genes. As one of the first studies to combine transcriptomics, morphology and physiology to characterize the reaction norms of C. virginica populations exposed to low salinity, our study provides insights into population-specific responses to environmental change.

Oyster collection and conditioning

In December 2017 and January 2018, 200 subtidal adult Crassostrea virginica (Gmelin 1791) oysters were collected from two estuaries in Louisiana, USA (Figs 1 and 2B); Calcasieu Lake (29°50′58′′N, 93°17′1′′W, mean monthly salinity 16.2±2.8 ppt) and Vermilion Bay (29°34′7′′N, 92°2′4′′W, mean monthly salinity 7.4±1.6 ppt). Additionally, in August 2018, 150 oysters were collected from two estuaries in Texas, USA (Figs 1 and 2B); Packery Channel (27°37′38′′N, 97°13′59′′W, mean monthly salinity 35.5±5.1 ppt) and Aransas Bay (28°7′38′′N, 96°59′8′′W, mean monthly salinity 23.0±6.9 ppt). In August 2018, all four broodstock populations were naturally induced to spawn at the Auburn University Shellfish Laboratory (AUSL) in Dauphin Island, AL, USA, as described in Marshall et al. (2021a,b). In brief, oysters were spawned in individual containers, and gametes were collected for controlled crosses. The crosses ranged from 7 males×11 females to 12 males×38 females depending on the population (Marshall et al., 2021a,b). Larvae were pooled by population, and settled oyster spat were maintained in upwelling nursery systems until ∼6 mm in shell height, and subsequently deployed in bags at an AUSL-permitted grow-out site at Bayou Sullivan, AL, USA (30°21′52′′N, 88°12′57′′W, mean daily salinity during grow-out 20.5±5.7 ppt) 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. Monthly mean salinity (±s.d.) for this site ranged from 9.0±2.4 ppt in May to 25.0±1.0 ppt in September (see Marshall et al., 2021a,b, for more 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, LA, USA.

Fig. 1.

Map of the broodstock collection sites. Eastern oysters, Crassostrea virginica, were collected from Packery Channel (PC), Aransas Bay (AB), Calcasieu Lake (CL) and Vermillion Bay (VB).

Fig. 1.

Map of the broodstock collection sites. Eastern oysters, Crassostrea virginica, were collected from Packery Channel (PC), Aransas Bay (AB), Calcasieu Lake (CL) and Vermillion Bay (VB).

Fig. 2.

Temperature and salinity at the collection sites. Average monthly temperature (A), salinity (B), and upper (solid lines) and lower (dashed lines) 5th percentile of temperature (C) and salinity (D) for the 10 years prior to oyster collection at the four sites.

Fig. 2.

Temperature and salinity at the collection sites. Average monthly temperature (A), salinity (B), and upper (solid lines) and lower (dashed lines) 5th percentile of temperature (C) and salinity (D) for the 10 years prior to oyster collection at the four sites.

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 (6 ppt or 12 ppt). We used artificial seawater or carbon-filtered and aerated freshwater to adjust the salinity. Oysters were fed daily at ∼5% of their meat dry mass 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 because of 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, and prevented us from collecting physiological data for all individuals (Table S1). In July 2020, we measured the final shell height, shell width, shell mass and dry meat mass of all oysters for each population at each salinity (Table 1). 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 mass was measured after drying at 70°C for 48 h. For gene expression analysis, we sampled 5 oysters per population per salinity (n=40). An approximately 0.5 cm2 piece of gill tissue was dissected and immediately frozen in liquid nitrogen and placed in a −80°C freezer until later extraction.

Table 1.

Average shell height, dry meat mass, clearance rate and oxygen consumption rate for oysters at 6 ppt and 12 ppt salinity

Average shell height, dry meat mass, clearance rate and oxygen consumption rate for oysters at 6 ppt and 12 ppt salinity
Average shell height, dry meat mass, clearance rate and oxygen consumption rate for oysters at 6 ppt and 12 ppt salinity

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 2 l beakers filled with 2 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×104 cells ml−1. Particle concentrations (>5 μm) were counted every 30 s 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, in l h−1) was calculated according to the equation CRi=[(bb′)×V×60], where b is the slope of the linear regression between the natural logarithm of cell concentration (cells ml−1) and time (min) in a beaker containing an oyster, b′ is the slope for the control beaker, and V is the volume of seawater (in l) in the beaker (2 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 r2>0.90 were included in these calculations.

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=(Hstd/Hexp)b×CRi, where CRh (l h−1 80 mm−1) 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).

Measuring oxygen consumption rate

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-3m (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 (V̇O2,i, in mg O2 h−1) was calculated as V̇O2,i=[(bb′)×V×60], where b is the slope of the linear regression of oxygen concentration (mg l−1) versus time for the chamber containing the oyster, b′ is the slope for the control, and V is the volume of water (in l) in the chamber (chamber minus oyster volume). Only linear regressions with an r2>0.95 were included in the calculations. Oxygen consumption rates were also standardized by dry meat mass (V̇O2,m, in mg O2 h−1 g−1), specifically to a standard oyster of 1 g dry meat mass. We used the equation V̇O2,m=(Mstd/Mexp)b×V̇O2,i, where Mstd is the dry meat mass of the standard oyster (1 g), Mexp is the dry meat mass of the experimental oyster, and b is the allometric exponent, 0.58 (Casas et al., 2018).

Statistical analyses

We examined the effect of population, salinity and the population×salinity interaction on clearance rates (CRh), oxygen consumption rates (V̇O2,m), and dry meat mass. We used a Shapiro–Wilk normality test (Shapiro and Wilk, 1965; Royston, 1982) to determine whether 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<0.05). To examine pairwise multiple comparisons, we used the post hoc Dunn's test with the Benjamini–Hochberg correction (P<0.05) (Benjamini and 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 (Puritz et al., 2022) 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 principal component analysis, 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 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 and Oshlack, 2010). Broad changes in gene expression were first assessed using a 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 and Hochberg, 1995). Functional enrichment of DEGs was tested using a Fisher's exact test (P<0.05). All Fisher's exact tests were conducted using the scripts originally developed by Wright et al. (2015) (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 whether 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 two permutational multivariate analysis of variance tests (PERMANOVA) to examine the effects of population, salinity or the population×salinity interaction on gene expression. The first PERMANOVA tested for an effect on total gene expression (averaging across all genes), while the second PERMANOVA tested on a per-gene basis so we could identify the specific genes associated with population, salinity or population×salinity interaction. 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 106 permutations (total expression) or 99,999 permutations (per-gene expression) (P<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 and Hochberg, 1995). Functional enrichment of significant genes was tested using a Fisher's exact test (P<0.05).

Gene expression analysis: WGCNA

We used a weighted gene co-expression network analysis (WGCNA) (Langfelder and 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 mass (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 7575 genes (7533 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<0.05).

Morphology and physiology

To better understand the effect of population and salinity on oyster performance, we measured dry meat mass (n=69), oxygen consumption rate (n=69) and clearance rate (n=44). We found that population of origin was correlated with clearance rate; however, a post hoc Dunn's test did now show any significant differences between pairwise comparisons. Salinity was negatively correlated with dry meat mass, indicating that oysters from all populations had more meat tissue in the lowest salinity treatment (6 ppt) (Fig. 3A, Tables 1 and 2). In contrast, clearance rates were in general positively correlated with salinity, which indicates reduced feeding in the lowest salinity treatment (Fig. 3B, Tables 1 and 2. Although clearance rates were affected by both population of origin and salinity, the strength of the correlation was greater for salinity as was the trend for all measured traits. We did not observe any interactive effects of population×salinity for any traits (Table 2). An effect of population or salinity on oxygen consumption rate was not found (Table 2).

Fig. 3.

Effects of salinity and population on dry meat mass and clearance rate of Crassostrea virginica. Box plots (median, 1.5× interquartile range and outliers) show the effect of salinity on oyster dry meat mass (n=69; A) and population+salinity on clearance rate (standardized by average shell height, 80 mm; n=44; B).

Fig. 3.

Effects of salinity and population on dry meat mass and clearance rate of Crassostrea virginica. Box plots (median, 1.5× interquartile range and outliers) show the effect of salinity on oyster dry meat mass (n=69; A) and population+salinity on clearance rate (standardized by average shell height, 80 mm; n=44; B).

Table 2.

Effect of salinity, population and their interaction on dry meat mass, clearance rate and oxygen consumption rate based on a two-way ANOVA test

Effect of salinity, population and their interaction on dry meat mass, clearance rate and oxygen consumption rate based on a two-way ANOVA test
Effect of salinity, population and their interaction on dry meat mass, clearance rate and oxygen consumption rate based on a two-way ANOVA test

Summary statistics of TagSeq

Transcriptome sequencing using TagSeq produced a total of 299 million reads, with an average of 8.1 million reads per sample (range: 1,542,904–61,628,689) after trimming. This depth is more than sufficient for the TagSeq method (Meyer et al., 2011). Star mapping resulted in 86.36% of reads mapping to the reference genome (uniquely mapped reads and multi-mapped reads averaged 63.98% and 22.37%, respectively). Filtering mapped reads based on expression resulted in a total of 20,138 genes kept for all downstream analyses.

Gene expression analysis: edgeR and PERMANOVA

PCA revealed a large effect of population on gene expression, which was driven largely by oysters from Packery Channel (Fig. 4). We examined the effect of population, salinity and the interactive effect of population×salinity on total gene expression using a PERMANOVA test, which identified a significant effect of population (P=1e−06), and no significant effect of salinity (P=0.09) or the population×salinity interaction (P=0.52) (Table S2).

Fig. 4.

Principal component analysis (PCA) plot showing distances between samples. The PCA was performed on gene counts that had been filtered and normalized using edgeR (n=37 samples). The shapes represent the two salinity treatments (circles, 6 ppt; squares, 12 ppt) and colors differentiate source populations.

Fig. 4.

Principal component analysis (PCA) plot showing distances between samples. The PCA was performed on gene counts that had been filtered and normalized using edgeR (n=37 samples). The shapes represent the two salinity treatments (circles, 6 ppt; squares, 12 ppt) and colors differentiate source populations.

Differences in gene expression between populations and salinities were identified using both edgeR and a PERMANOVA test on a per-gene basis. For edgeR, we first compared the gene expression profiles of individuals from one population versus those of individuals from the other three populations. We detected 264 genes differentially expressed in Vermilion Bay (195 up-regulated and 69 down-regulated), 259 genes in Calcasieu Lake (211 up-regulated and 48 down-regulated), 51 genes in Aransas Bay (43 up-regulated and 8 down-regulated) and 3158 genes in Packery Channel (1670 up-regulated and 1488 down-regulated). The PERMANOVA test identified 2473 genes associated with population, with 2138 of those genes overlapping with the edgeR results (Fig. S1). GO enrichment detected only one molecular function, galactosyltransferase activity, for the genes identified by the PERMANOVA (Table S3). No enrichments were detected for genes identified by edgeR.

Next, we tested for genes differing by salinity treatment. For edgeR, we compared the gene expression profiles of oysters from the 6 ppt treatment versus the 12 ppt treatment. We examined this for all oysters regardless of population, on a per-population basis, and also looked at Packery Channel separately plus the combination of the remaining three populations. edgeR detected 20 genes differing by salinity when looking at all oysters, 26 genes when looking at Packery Channel separately plus the combination of the remaining three populations, and when focusing on each population independently, each had only one gene differing between salinities. Of the 26 genes responding to salinity in Packery Channel and the remaining three populations, the majority had greater fold-change in expression in Aransas Bay, Calcasieu Lake and Vermilion Bay than in Packery Channel (Fig. 5). The PERMANOVA test identified 13 genes differentially expressed between salinities, with 9 of those genes overlapping with the edgeR results (Fig. S1). GO enrichment did not detect any enrichments for genes identified by edgeR or PERMANOVA. Additionally, the PERMANOVA test did not detect any genes with a significant population×salinity interaction.

Fig. 5.

log Fold-change in expression for transcripts differentially expressed in response to salinity (6 versus 12 ppt). The one-to-one line divides transcripts with greater fold-change in expression in Aransas Bay (AB), Calcasieu Lake (CL) and Vermilion Bay (VB) than in Packery Channel (PC) (n=37).

Fig. 5.

log Fold-change in expression for transcripts differentially expressed in response to salinity (6 versus 12 ppt). The one-to-one line divides transcripts with greater fold-change in expression in Aransas Bay (AB), Calcasieu Lake (CL) and Vermilion Bay (VB) than in Packery Channel (PC) (n=37).

Gene expression analysis: WGCNA

The WGCNA analysis was used to detect groups of genes responding similarly to salinity, population, final shell height, shell width, shell mass, dry meat mass, clearance rate and oxygen consumption rate (n=37). A gene co-expression matrix was generated with modules that represent similarly expressed genes. We found that 12 out of 19 modules were significantly correlated with population, 1 module with salinity, 2 modules with shell height, 2 modules with shell width, 1 module with shell mass, 3 modules with dry meat mass and 2 modules with oxygen consumption rate (Fig. S2). When limiting the analysis to individuals with clearance rate measured (n=23), we found 2 of the 17 modules were significantly correlated with clearance rate (Fig. S3). Functional enrichment revealed gene ontologies for population, shell height, shell width and oxygen consumption rate (Table 3).

Table 3.

Top enriched GO terms identified by WGCNA analysis

Top enriched GO terms identified by WGCNA analysis
Top enriched GO terms identified by WGCNA analysis

We present one of the first studies to combine transcriptomics, morphology and physiology to characterize the reaction norms of C. virginica populations reared in a common garden field site and subsequently exposed to low salinity. By using a comprehensive approach with three unique transcriptomic analyses combined with morphological and physiological data, we were able to identify the effect of genotype (population), environment (salinity) and the genotype×environment interaction on both whole-organism and molecular phenotypes. Our RNAseq analyses (edgeR, PERMANOVA and WGCNA) largely agree in finding that gene expression was mainly affected by the population of origin and minimally affected by differences in salinity. In contrast, the morphological and physiological data suggest that low salinity has a larger influence on oyster performance than the population of origin, but that the negative effects of low salinity did not differ by population. Our gene expression, morphological and physiological data all consistently lacked signatures of a genotype×environment interaction, and more specifically local adaptation, for the salinity conditions tested (6 ppt and 12 ppt). Overall, we found that the steady state of C. virginica shows large evolved differences in baseline gene expression but almost no expression response to sustained low salinity, and a plastic physiological response to salinity that is similar across populations.

Morphology and physiology: plastic response does not differ among populations

The ANOVA indicated an effect of population of origin on clearance rates (Table 2). A post hoc Dunn's test did now show any significant differences between pairwise comparisons, but there was a trend towards lower clearance rates at 12 ppt for the Packery Channel population, relative to the other three. We observed a highly plastic response for dry meat mass and clearance rate in response to low salinity, but that plasticity did not differ across populations (no G×E) (Table 2, Fig. 3), suggesting populations have not optimized their plastic response to their local salinities. (Fig. 3A). However, one caveat is that our limited sample size and narrow range of salinity exposure may be concealing a stronger population signal or signs of a genotype×environment interaction.

The lack of evidence for population-specific responses to low salinity is somewhat surprising, given that Marshall et al. (2021a,b) found that the population originating from the highest salinity estuary (Packery Channel) had significantly higher mortality and smaller morphology under low salinity conditions. However, this previous study only found signs of local adaptation occurring at the lowest and highest salinity treatments (44 ppt and 2 ppt), which is outside the range of our exposure and may suggest that population-specific responses are only observed at the most extreme salinities or when salinity differences are larger than in our study. Conversely, Heilmayer et al. (2008) support our results in finding no evidence for population-specific responses when examining the interactive effects of salinity×temperature on condition index or cellular stress when exposing oysters to equally low salinity (2 ppt) and medium-high salinity (25 ppt). Furthermore, in agreement with our findings, Hughes et al. (2017) did not observe increased survival or growth for local populations when conducting a reciprocal transplant with six populations spread across three field sites in the southeastern Atlantic Bight. Overall, our morphological and physiological data suggest that oysters in this region largely use plasticity to buffer the effects of salinity changes and are less likely to buffer these changes through evolved differences, or more specifically local adaptation.

We were surprised to find that all four oyster populations had more meat tissue in the lowest salinity treatment (6 ppt), especially as other studies have found decreases in body condition index in lower salinity treatments (Wilson et al., 2005; Heilmayer et al., 2008). However, one possible explanation for this is that oysters spawned in the 12 ppt tanks, while no spawning was observed for oysters in the 6 ppt tank. As spawning largely depletes oyster resources, it is plausible that oysters allocated more energy towards reproductive efforts when exposed to a more favorable salinity (>6 ppt) and that oysters in the lowest salinity treatment were more stressed and not directing as many resources towards spawning efforts. This seems likely when considering that feeding rates were lowest in the 6 ppt treatment, which suggests that these oysters were not performing as well as oysters in higher salinity. Lastly, spawning itself can largely explain the differences in dry meat mass as the gonad can make a large part of the total mass of the oysters and is significantly reduced when gametes are released.

Gene expression analyses: evolved differences but lacking plasticity

Interestingly, in contrast to the morphological and physiological data, most differential expression occurred between populations (evolved response) with minimal differences between salinity treatments (plastic response) (Fig. 4). The most prominent differences between populations came from comparisons between Packery Channel and the remaining three populations (Fig. 4). The large difference in gene expression between individuals from Packery Channel and the other three populations is somewhat expected based on the population structure of eastern oysters in this region. Previous studies using allozymes, microsatellites and single nucleotide polymorphism (SNP) arrays have observed a consistent genetic break between C. virginica in southern Texas versus the remaining nGOM, with a transition zone possibly near Corpus Christi (Buroker, 1983; King et al., 1994; Varney et al., 2009; Anderson et al., 2014; Thongda et al., 2018; Hajovsky et al., 2021). Furthermore, this genetic break is shared with many other marine organisms including Atlantic croakers, black drum and fiddler crabs (Barnwell and Thurman, 1984; Woltmann et al., 2014; Anderson et al., 2019; Williford et al., 2021). The reason for this genetic divergence remains unclear, but some possibilities include a hydrogeographic barrier occurring near Laguna Madre (King et al., 1994), or potentially natural selection for increased tolerance to high salinity in southern Texas (Hedgpeth, 1953; Breuer, 1962; Eierman and Hare, 2016).

The separation between Packery Channel and the remaining three populations is also reflected in log fold-changes in expression between 6 and 12 ppt, where gene expression was in a similar direction for all populations, but the magnitude of change was smaller for oysters from Packery Channel, which might have already been stressed at 12 ppt (Fig. 5). This is plausible as Marshall et al. (2021a,b) observed lower growth and higher mortality for oysters from Packery Channel compared with the other three populations when exposed to a low salinity site. However, although individuals from Packery Channel showed divergence in baseline gene expression, we identified similar levels of plasticity in expression under low salinity. This suggests that the difference in performance between populations exposed to low salinity is unlikely to be the result of adaptive changes in differential gene expression.

One potential caveat is that RNA sequencing only captures a snapshot in time, so our failure to detect a large effect of salinity on gene expression may be due to sampling mRNA at a time that was not relevant to these traits, as opposed to the other studies which observed these differences. Like many osmoconformers, oysters use multiple mechanisms to regulate cell volume under changing salinities, and these mechanisms are deployed over different time scales. Short-term responses involve transport of inorganic ions (an emergency response to prevent dangerous swelling or shrinking of cells) (Pierce, 1982). This emergency response is followed by the accumulation or release of organic osmolytes, including free amino acids and quaternary amines to establish osmotic balance (Hosoi et al., 2003). Once this balance is established, the animal will not need to perform an additional osmoregulation, assuming the salinity remains stable at its new level. This may explain why we observed relatively few DEGs across salinities: because we are measuring gene expression at the new steady state for each salinity, we likely failed to detect any adaptive changes in gene expression occurring in the immediate acclimation period. Studies on other oyster species generally found large transcriptomic responses occurring within the first 8 days of exposure, but these differences may diminish as individuals reach a new steady state (Zhao et al., 2012; Meng et al., 2013; Maynard et al., 2018). Additionally, differences in gene expression between salinities could have been a result of alternatively spliced isoforms which are not captured by the TagSeq method. Even with these caveats in mind, our data suggest that oyster populations in the nGOM broadly respond to salinity in similar ways, with little in the way of salinity×population effects, either for physiology or for gene expression.

Gene expression; the potential driver of differences in morphology, feeding and oxygen consumption rates

WGCNA analysis was able to identify clusters of genes with highly correlated expression patterns associated with final shell width, shell mass, dry meat mass, oxygen consumption rate and clearance rate. Subsequent GO analysis gave us the ability to discover biologically relevant categories for the gene clusters associated with some of these traits (Table S3). For example, we discovered aminoacyl-tRNA enzymes were among the GO terms associated with oxygen consumption rate. These aminoacyl-tRNA enzymes have been previously linked to shell repair in pearl oysters and have been shown to have a role in hard clams in response to air exposure (Xiong et al., 2021; Zhou et al., 2021). Additionally, 10 of the 25 GO enrichments detected were associated with Packery Channel, further reflecting the genetic split between oysters in southern Texas and the remaining nGOM. Our ability to detect genes and relevant gene ontologies associated with morphological and physiological differences between individuals suggests that gene expression may be an important driver of differences in morphology, feeding and oxygen consumption rates in C. virginica.

Conclusions

Our study focused on the role of plastic and evolved responses to low salinity in C. virginica using measurements of morphology, physiology and 3′ RNA sequencing. We demonstrate that C. virginica exhibits a highly evolved response to low salinity when examining gene expression and a relatively small response when focusing on morphology and physiology. In contrast to Marshall et al. (2021a,b), we found no evidence for population-specific responses among our four populations. Individuals from the highest salinity estuary (Packery Channel) displayed highly divergent gene expression from that of the other three populations though, which may indicate a potential for population-specific responses or evolved differences in other traits besides the response to salinity. Our findings suggest that oysters in this region largely rely on plasticity in morphology and physiology to buffer the effects of low salinity; the lack of observed differences in gene expression may be driven by the fact that oysters perform little active osmoregulation once they have conformed to a new stable salinity.

We would like to thank Dr Jennifer Pollack for her help in project design, grant writing and collecting oysters from the two Texas estuaries. We would also like to thank Scott Rikard at the Auburn University Shellfish Laboratory in Dauphin Island, Alabama, for help in spawning and rearing oysters, Sarah Spellman and Caitlin Robitaille (Auburn University) for assistance in hatchery production, and Glen Chaplin (Auburn University) for field grow-out of the oysters. Portions of this research were conducted with high-performance computing resources provided by Louisiana State University (http://www.hpc.lsu.edu). We would also like to thank Maggie Thomas, Phuong Duong and Kenny-Khai Vo for their assistance in RNA extractions.

Author contributions

Conceptualization: J.F.L.P., M.W.K.; Methodology: K.A.S., S.M.C., J.F.L.P., M.W.K.; Validation: K.A.S.; Formal analysis: K.A.S.; Investigation: K.A.S.; Resources: J.F.L.P., M.W.K.; Data curation: K.A.S.; Writing - original draft: K.A.S.; Writing - review & editing: K.A.S., S.M.C., J.F.L.P., M.W.K.; Supervision: M.W.K.; Funding acquisition: M.W.K.

Funding

This research was supported by National Science Foundation BioOCE 1731710 to M.W.K. and J.F.L.P., Louisiana Sea Grant Award NA14OAR4170099 to M.W.K. and J.F.L.P., and the Alfred P. Sloan Foundation research fellowship awarded to M.W.K. During the project, K.A.S. was supported by a National Science Foundation Graduate Research Fellowship under grant no. 00001414.

Data availability

The count data, oyster traits, list of GO terms and scripts used for analyses are available from Dryad (Sirovy et al., 2022): https://doi.org/10.5061/dryad.wh70rxwq9. The raw sequencing data can be found in the NCBI SRA database under BioProject submission SUB11217727.

Anderson
,
J. D.
,
Karel
,
W. J.
,
Mace
,
C. E.
,
Bartram
,
B. L.
and
Hare
,
M. P.
(
2014
).
Spatial genetic features of eastern oysters (Crassostrea virginica Gmelin) in the Gulf of Mexico: northward movement of a secondary contact zone
.
Ecol. Evol.
4
,
1671
-
1685
.
Anderson
,
J. D.
,
O'leary
,
S. J.
and
Cooper
,
P. T.
(
2019
).
Population Structure of Atlantic croakers from the Gulf of Mexico: evaluating a single–stock hypothesis using a genomic approach
.
Mar. Coastal Fish.
11
,
3
-
16
.
Barnwell
,
F. H.
and
Thurman
,
C. L.
(
1984
).
Taxonomy and biogeography of the fiddler crabs (Ocypodidae: Genus Uca) of the Atlantic and Gulf coasts of eastern North America
.
Zool. J. Linn. Soc.
81
,
23
-
87
.
Bayne
,
B. L.
(
2017
).
Biology of Oysters
.
Academic Press
.
Benjamini
,
Y.
and
Hochberg
,
Y.
(
1995
).
Controlling the false discovery rate: a practical and powerful approach to multiple testing
.
J. R. Stat. Soc. Ser. B (Methodological)
57
,
289
-
300
.
Bishop
,
D. A.
,
Williams
,
A. P.
and
Seager
,
R.
(
2019
).
Increased fall precipitation in the southeastern united states driven by higher-intensity, frontal precipitation
.
Geophys. Res. Lett.
46
,
8300
-
8309
.
Bolger
,
A. M.
,
Lohse
,
M.
and
Usadel
,
B.
(
2014
).
Trimmomatic: a flexible trimmer for Illumina sequence data
.
Bioinformatics
30
,
2114
-
2120
.
Breitburg
,
D. L.
,
Salisbury
,
J.
,
Bernhard
,
J. M.
,
Cai
,
W.-J.
,
Dupont
,
S.
,
Doney
,
S. C.
,
Kroeker
,
K.
,
Levin
,
L.
,
Long
,
W. C.
,
Milke
,
L.
et al. 
(
2015
).
And on top of all that… Coping with ocean acidification in the midst of many stressors
.
Oceanography
28
,
48
-
61
.
Breuer
,
J. P.
(
1962
).
An ecological survey of the lower Laguna Madre of Texas, 1953-1959
.
Publ. Inst. Mar. Sci. Univ. Texas
8
,
3
-
183
.
Buroker
,
N. E.
(
1983
).
Population genetics of the American oyster Crassostrea virginica along the Atlantic coast and the Gulf of Mexico
.
Mar. Biol.
75
,
99
-
112
.
Casas
,
S. M.
,
Lavaud
,
R.
,
La Peyre
,
M. K.
,
Comeau
,
L. A.
,
Filgueira
,
R.
and
La Peyre
,
J. F.
(
2018
).
Quantifying salinity and season effects on eastern oyster clearance and oxygen consumption rates
.
Mar. Biol.
165
,
90
.
Cayan
,
D. R.
,
Das
,
T.
,
Pierce
,
D. W.
,
Barnett
,
T. P.
,
Tyree
,
M.
and
Gershunov
,
A.
(
2010
).
Future dryness in the Southwest US and the hydrology of the early 21st century drought
.
Proc. Natl. Acad. Sci. USA
107
,
21271
-
21276
.
Coughlan
,
J.
(
1969
).
The estimation of filtering rate from the clearance of suspensions
.
Mar. Biol.
2
,
356
-
358
.
Cranford
,
P. J.
,
Ward
,
J. E.
and
Shumway
,
S. E.
(
2011
).
Bivalve filter feeding: variability and limits of the aquaculture biofilter
. In
Shellfish Aquaculture and the Environment
(ed. S. E. Shumway), pp.
81
-
124
.
Cranford
,
P. J.
,
Strohmeier
,
T.
,
Filgueira
,
R.
and
Strand
,
Ø.
(
2016
).
Potential methodological influences on the determination of particle retention efficiency by suspension feeders: Mytilus edulis and Ciona intestinalis
.
Aquat. Biol.
25
,
61
-
73
.
Das
,
A.
,
Justic
,
D.
,
Inoue
,
M.
,
Hoda
,
A.
,
Huang
,
H.
and
Park
,
D.
(
2012
).
Impacts of Mississippi River diversions on salinity gradients in a deltaic Louisiana estuary: ecological and management implications
.
Estuarine Coastal Shelf Sci.
111
,
17
-
26
.
De Jong
,
G.
(
2005
).
Evolution of phenotypic plasticity: patterns of plasticity and the emergence of ecotypes
.
New Phytol.
166
,
101
-
118
.
DeBiasse
,
M. B.
and
Kelly
,
M. W.
(
2016
).
Plastic and evolved responses to global change: what can we learn from comparative transcriptomics?
J. Hered.
107
,
71
-
81
.
Dewitt
,
T. J.
and
Scheiner
,
S. M.
(
2004
).
Phenotypic Plasticity: Functional and Conceptual Approaches
.
New York
,
USA
:
Oxford University Press
.
Dobin
,
A.
,
Davis
,
C. A.
,
Schlesinger
,
F.
,
Drenkow
,
J.
,
Zaleski
,
C.
,
Jha
,
S.
,
Batut
,
P.
,
Chaisson
,
M.
and
Gingeras
,
T. R.
(
2013
).
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics
29
,
15
-
21
.
Doney
,
S. C.
,
Ruckelshaus
,
M.
,
Emmett Duffy
,
J.
,
Barry
,
J. P.
,
Chan
,
F.
,
English
,
C. A.
,
Galindo
,
H. M.
,
Grebmeier
,
J. M.
,
Hollowed
,
A. B.
,
Knowlton
,
N.
et al. 
(
2012
).
Climate change impacts on marine ecosystems
.
Annu. Rev. Mar. Sci.
4
,
11
-
37
.
Eierman
,
L. E.
and
Hare
,
M. P.
(
2016
).
Reef-specific patterns of gene expression plasticity in eastern oysters (Crassostrea virginica)
.
J. Hered.
107
,
90
-
100
.
Forsman
,
A.
(
2015
).
Rethinking phenotypic plasticity and its consequences for individuals, populations and species
.
Heredity
115
,
276
-
284
.
Galtsoff
,
P. S.
(
1964
).
The American Oyster, Crassostrea virginica (Gmelin)
.
Fishery Bulletin of the Fish and Wildlife Service
64
,
219
-
238
.
Ghalambor
,
C. K.
,
Mckay
,
J. K.
,
Carroll
,
S. P.
and
Reznick
,
D. N.
(
2007
).
Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments
.
Funct. Ecol.
21
,
394
-
407
.
Gómez-Chiarri
,
M.
,
Warren
,
W. C.
,
Guo
,
X.
and
Proestou
,
D.
(
2015
).
Developing tools for the study of molluscan immunity: The sequencing of the genome of the eastern oyster, Crassostrea virginica
.
Fish Shellfish Immunol.
46
,
2
-
4
.
Gosling
,
E.
(
2015
).
Marine Bivalve Molluscs
.
John Wiley & Sons
.
Gunderson
,
A. R.
,
Armstrong
,
E. J.
and
Stillman
,
J. H.
(
2016
).
Multiple stressors in a changing world: the need for an improved perspective on physiological responses to the dynamic marine environment
.
Annu. Rev. Mar. Sci.
8
,
357
-
378
.
Hajovsky
,
P.
,
Beseres Pollack
,
J.
and
Anderson
,
J.
(
2021
).
Morphological Assessment of the Eastern Oyster Crassostrea virginica throughout the Gulf of Mexico
.
Mar. Coast. Fish.
13
,
309
-
319
.
Hedgpeth
,
J. W.
(
1953
).
An introduction to the zoogeography of the northwestern Gulf of Mexico with reference to the invertebrate fauna
.
Publications of the Institute of Marine Science (University of Texas)
3
,
107
-
224
.
Heilmayer
,
O.
,
Digialleonardo
,
J.
,
Qian
,
L.
and
Roesijadi
,
G.
(
2008
).
Stress tolerance of a subtropical Crassostrea virginica population to the combined effects of temperature and salinity
.
Estuarine Coast. Shelf Sci.
79
,
179
-
185
.
Hendry
,
A. P.
(
2016
).
Key questions on the role of phenotypic plasticity in eco-evolutionary dynamics
.
J. Hered.
107
,
25
-
41
.
Hoffmann
,
A. A.
and
Sgró
,
C. M.
(
2011
).
Climate change and evolutionary adaptation
.
Nature
470
,
479
-
485
.
Hosoi
,
M.
,
Kubota
,
S.
,
Toyohora
,
M.
,
Toyohora
,
H.
and
Hayashi
,
I.
(
2003
).
Effect of salinity change on free amino acid content in Pacific oyster
.
Fish Sci.
69
,
395
-
400
.
Hughes
,
A. R.
,
Hanley
,
T. C.
,
Byers
,
J. E.
,
Grabowski
,
J. H.
,
Malek
,
J. C.
,
Piehler
,
M. F.
and
Kimbro
,
D. L.
(
2017
).
Genetic by environmental variation but no local adaptation in oysters (Crassostrea virginica)
.
Ecol. Evol.
7
,
697
-
709
.
King
,
T. L.
,
Ward
,
R.
and
Zimmerman
,
E. G.
(
1994
).
Population structure of eastern oysters (Crassostrea virginica) inhabiting the Laguna Madre, Texas, and adjacent bay systems
.
Can. J. Fish. Aquat. Sci.
51
,
215
-
222
.
Koch
,
E. L.
and
Guillaume
,
F.
(
2020
).
Additive and mostly adaptive plastic responses of gene expression to multiple stress in Tribolium castaneum
.
PLoS Genet.
16
,
e1008768
.
La Peyre
,
M. K.
,
Gossman
,
B.
and
La Peyre
,
J. F.
(
2009
).
Defining optimal freshwater flow for oyster production: effects of freshet rate and magnitude of change and duration on eastern oysters and Perkinsus marinus infection
.
Estuaries Coasts
32
,
522
-
534
.
La Peyre
,
M. K.
,
Aguilar Marshall
,
D.
,
Miller
,
L. S.
and
Humphries
,
A. T.
(
2019
).
Oyster reefs in northern Gulf of Mexico estuaries harbor diverse fish and decapod crustacean assemblages: a meta-synthesis
.
Front. Mar. Sci.
6
,
666
.
Langfelder
,
P.
and
Horvath
,
S.
(
2008
).
WGCNA: An R package for weighted correlation network analysis
.
BMC Bioinformatics
9
,
1
-
13
.
Lavaud
,
R.
,
La Peyre
,
M. K.
,
Justic
,
D.
and
La Peyre
,
J. F.
(
2021
).
Dynamic energy budget modelling to predict eastern oyster growth, reproduction, and mortality under river management and climate change scenarios
.
Estuarine Coast. Shelf Sci.
251
,
107188
.
Levine
,
M. T.
,
Eckert
,
M. L.
and
Begun
,
D. J.
(
2011
).
Whole-genome expression plasticity across tropical and temperate Drosophila melanogaster populations from eastern Australia
.
Mol. Biol. Evol.
28
,
249
-
256
.
Marshall
,
D. A.
,
Coxe
,
N. C.
,
La Peyre
,
M. K.
,
Walton
,
W. C.
,
Rikard
,
F. S.
,
Pollack
,
J. B.
,
Kelly
,
M. W.
and
La Peyre
,
J. F.
(
2021a
).
Tolerance of northern Gulf of Mexico eastern oysters to chronic warming at extreme salinities
.
J. Therm. Biol.
100
,
103072
.
Marshall
,
D. A.
,
Casas
,
S. M.
,
Walton
,
W. C.
,
Rikard
,
F. S.
,
Palmer
,
T. A.
,
Breaux
,
N.
,
La Peyre
,
M. K.
,
Pollack
,
J. B.
,
Kelly
,
M.
and
La Peyre
,
J. F.
(
2021b
).
Divergence in salinity tolerance of northern Gulf of Mexico eastern oysters under field and laboratory exposure
.
Conserv. Physio.
9
,
coab065
.
Maynard
,
A.
,
Bible
,
J. M.
,
Pespeni
,
M. H.
,
Sanford
,
E.
and
Evans
,
T. G.
(
2018
).
Transcriptomic responses to extreme low salinity among locally adapted populations of Olympia oyster (Ostrea lurida)
.
Mol. Ecol.
27
,
4225
-
4240
.
Meng
,
J.
,
Zhu
,
Q.
,
Zhang
,
L.
,
Li
,
C.
,
Li
,
L.
,
She
,
Z.
,
Huang
,
B.
and
Zhang
,
G.
(
2013
).
Genome and transcriptome analyses provide insight into the euryhaline adaptation mechanism of Crassostrea gigas
.
PLoS ONE
8
,
e58563
.
Meyer
,
E.
,
Aglyamova
,
G. V.
and
Matz
,
M. V.
(
2011
).
Profiling gene expression responses of coral larvae (Acropora millepora) to elevated temperature and settlement inducers using a novel RNA-Seq procedure
.
Mol. Ecol.
20
,
3599
-
3616
.
Nevins
,
J. A.
,
Pollack
,
J. B.
and
Stunz
,
G. W.
(
2014
).
Characterizing nekton use of the largest unfished oyster reef in the United States compared with adjacent estuarine habitats
.
J. Shellfish Res.
33
,
227
-
238
.
Parker
,
R. H.
(
1960
).
Ecology and distributional patterns of marine macro-invertebrates, Northern Gulf of Mexico
. In:
Recent Sediments, Northwest Gulf of Mexico
(ed.
F. B.
Phleger
and
T. H.
van Andel
), pp.
302
-
337
.
Tulsa, Okla
:
Tulsa, American Association of Petroleum Geologists
.
Parmesan
,
C.
(
2006
).
Ecological and evolutionary responses to recent climate change
.
Annu. Rev. Ecol. Evol. Syst.
37
,
637
-
669
.
Pickett
,
S. T. A.
(
1989
).
Space-for-Time Substitution as an Alternative to Long-Term Studies
. In:
Long-Term Studies in Ecology
(ed. G. E. Likens), pp.
110
-
135
.
New York
:
Springer-Verlag
.
Pierce
,
S. K.
(
1982
).
Invertebrate cell volume control mechanisms: a coordinated use of intracellular amino acids and inorganic ions as osmotic solute
.
Biol. Bull.
163
,
405
-
419
.
Price
,
T. D.
,
Qvarnström
,
A.
and
Irwin
,
D. E.
(
2003
).
The role of phenotypic plasticity in driving genetic evolution
.
Proc. R. Soc. Lond. Ser. B Biol. Sci.
270
,
1433
-
1440
.
Puritz
,
J. B.
,
Guo
,
X.
,
Hare
,
M.
,
He
,
Y.
,
Hillier
,
L. W.
,
Jin
,
S.
,
Liu
M.
,
Lotterhos
,
K. E.
,
Minx
,
P.
and
Modak
,
T.
, et al. 
(
2023
).
A second unveiling: Haplotig masking of the eastern oyster genome improves population-level inference
.
Mol. Ecol. Resour
.
Robinson
,
M. D.
and
Oshlack
,
A.
(
2010
).
A scaling normalization method for differential expression analysis of RNA-seq data
.
Genome Biol.
11
,
R25
.
Robinson
,
M. D.
,
Mccarthy
,
D. J.
and
Smyth
,
G. K.
(
2010
).
edgeR: a Bioconductor package for differential expression analysis of digital gene expression data
.
Bioinformatics
26
,
139
-
140
.
Rockman
,
M. V.
(
2008
).
Reverse engineering the genotype–phenotype map with natural genetic variation
.
Nature
456
,
738
-
744
.
Royston
,
J. P.
(
1982
).
An Extension of Shapiro and Wilk's W test for normality to large samples
.
Appl. Stat.
31
,
115
.
Rybovich
,
M.
,
Peyre
,
M. K. L.
,
Hall
,
S. G.
and
Peyre
,
J. F. L.
(
2016
).
Increased temperatures combined with lowered salinities differentially impact oyster size class growth and mortality
.
J. Shellfish Res.
35
,
101
-
113
.
Sanford
,
E.
and
Kelly
,
M. W.
(
2011
).
Local adaptation in marine invertebrates
.
Annu. Rev. Mar. Sci.
3
,
509
-
535
.
Schmidtko
,
S.
,
Stramma
,
L.
and
Visbeck
,
M.
(
2017
).
Decline in global oceanic oxygen content during the past five decades
.
Nature
542
,
335
-
339
.
Seebacher
,
F.
,
White
,
C. R.
and
Franklin
,
C. E.
(
2015
).
Physiological plasticity increases resilience of ectothermic animals to climate change
.
Nat. Clim. Change
5
,
61
-
66
.
Shapiro
,
S. S.
and
Wilk
,
M. B.
(
1965
).
An analysis of variance test for normality (complete samples)
.
Biometrika
52
,
591
-
611
.
Sirovy
,
K. A.
,
Johnson
,
K. M.
,
Casas
,
S. M.
,
La Peyre
,
J. F.
and
Kelly
,
M. W.
(
2021
).
Lack of genotype-by-environment interaction suggests limited potential for evolutionary changes in plasticity in the eastern oyster, Crassostrea virginica
.
Mol. Ecol.
30
,
5721
-
5734
.
Sirovy
,
K.
,
Casas
,
S.
,
La Peyre
,
J.
,
Kelly
,
M.
(
2022
).
Data from: population-specific responses in eastern oysters exposed to low salinity in the northern Gulf of Mexico
. Dryad, Dataset.
Swam
,
L. M.
,
Couvillion
,
B.
,
Callam
,
B.
,
La Peyre
,
J. F.
and
La Peyre
,
M. K.
(
2022
).
Defining oyster resource zones across coastal Louisiana for restoration and aquaculture
.
Ocean Coast. Management
225
,
106178
.
Thongda
,
W.
,
Zhao
,
H.
,
Zhang
,
D.
,
Jescovitch
,
L. N.
,
Liu
,
M.
,
Guo
,
X.
,
Schrandt
,
M.
,
Powers
,
S. P.
and
Peatman
,
E.
(
2018
).
Development of SNP panels as a new tool to assess the genetic diversity, population structure, and parentage analysis of the eastern oyster (Crassostrea virginica)
.
Mar. Biotechnol.
20
,
385
-
395
.
Varney
,
R. L.
,
Galindo-Sánchez
,
C. E.
,
Cruz
,
P.
and
Gaffney
,
P. M.
(
2009
).
Population genetics of the eastern oyster Crassostrea virginica (Gmelin, 1791) in the Gulf of Mexico
.
J. Shellfish Res.
28
,
855
-
864
.
Via
,
S.
and
Lande
,
R.
(
1985
).
Genotype-environment interaction and the evolution of phenotypic plasticity
.
Evolution
39
,
505
-
522
.
Volety
,
A. K.
,
Haynes
,
L.
,
Goodman
,
P.
and
Gorman
,
P.
(
2014
).
Ecological condition and value of oyster reefs of the Southwest Florida shelf ecosystem
.
Ecol. Indicators
44
,
108
-
119
.
West-Eberhard
,
M. J.
(
2003
).
Developmental Plasticity and Evolution
.
New York
:
Oxford University Press
.
Widdows
,
J.
(
1985
).
Physiological procedures
. In:
The Effects of Stress and Pollution on Marine Animals
(ed.
B. L.
Bayne
,
D. A.
Brown
,
K.
Burns
,
D. R.
Dixon
,
A.
Ivanovici
and
D. R.
Livingstone
et al. ), pp. 161-178. New York: Praeger Scientific.
Williford
,
D.
,
Anderson
,
J.
,
Karel
,
W.
and
Olsen
,
Z.
(
2021
).
Phylogeography, population structure, and historical demography of black drum in North America
.
N. Am. J. Fish. Manag.
41
,
1020
-
1039
.
Wilson
,
C.
,
Scotto
,
L.
,
Scarpa
,
J.
,
Volety
,
A.
,
Laramore
,
S.
and
Haunert
,
D.
(
2005
).
Survey of water quality, oyster reproduction and oyster health status in the St. Lucie Estuary
.
J. Shellfish Res.
24
,
157
-
165
.
Woltmann
,
S.
,
Stouffer
,
P. C.
,
Bergeon Burns
,
C. M.
,
Woodrey
,
M. S.
,
Cashner
,
M. F.
and
Taylor
,
S. S.
(
2014
).
Population genetics of Seaside Sparrow (Ammodramus maritimus) subspecies along the Gulf of Mexico
.
PLoS ONE
9
,
e112739
.
Wright
,
R. M.
,
Aglyamova
,
G. V.
,
Meyer
,
E.
and
Matz
,
M. V.
(
2015
).
Gene expression associated with white syndromes in a reef building coral, Acropora hyacinthus
.
BMC Genomics
16
,
371
.
Wund
,
M. A.
(
2012
).
Assessing the impacts of phenotypic plasticity on evolution
.
Integr. Comp. Biol.
52
,
5
-
15
.
Xiong
,
X.
,
Cao
,
Y.
,
Li
,
Z.
,
Jiao
,
Y.
,
Du
,
X.
and
Zheng
,
Z.
(
2021
).
Transcriptome analysis reveals the transition and crosslinking of immune response and biomineralization in shell damage repair in pearl oyster
.
Aquacult. Rep.
21
,
100851
.
Zhao
,
X.
,
Yu
,
H.
,
Kong
,
L.
and
Li
,
Q.
(
2012
).
Transcriptomic responses to salinity stress in the Pacific oyster Crassostrea gigas
.
PLoS One
7
,
e46244
.
Zhou
,
S.
,
Campbell
,
T. G.
,
Stone
,
E. A.
,
Mackay
,
T. F. C.
and
Anholt
,
R. R. H.
(
2012
).
Phenotypic plasticity of the Drosophila transcriptome
.
PLoS Genet.
8
,
e1002593
.
Zhou
,
C.
,
Song
,
H.
,
Feng
,
J.
,
Hu
,
Z.
,
Yu
,
Z.-L.
,
Yang
,
M.-J.
,
Shi
,
P.
,
Li
,
Y.-R.
,
Guo
,
Y.-J.
and
Zhang
,
T.
(
2021
).
RNA-Seq analysis and WGCNA reveal dynamic molecular responses to air exposure in the hard clam Mercenaria mercenaria
.
Genomics
113
,
2847
-
2859
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information