Population-specific responses in eastern oysters exposed to low salinity in the northern Gulf of Mexico

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.

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.

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 cm² piece of gill tissue was dissected and immediately frozen in liquid nitrogen and placed in a −80°C freezer until later extraction.

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×10⁴ cells ml⁻¹. 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 (CR_i, in l h⁻¹) was calculated according to the equation CR_i=[(b−b′)×V×60], where b is the slope of the linear regression between the natural logarithm of cell concentration (cells ml⁻¹) 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 r²>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 CR_h=(H_std/H_exp)^b×CR_i, where CR_h (l h⁻¹ 80 mm⁻¹) is the clearance rate standardized by shell height, H_std is the shell height of the standard oyster (80 mm), H_exp is the shell height of the experimental oyster, CR_i is the individual oyster clearance rate, and b is the allometric exponent for shell height set to 1.78 (Cranford et al., 2011).

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 O₂ h⁻¹) was calculated as V̇_O2,i=[(b−b′)×V×60], where b is the slope of the linear regression of oxygen concentration (mg l⁻¹) 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 r²>0.95 were included in the calculations. Oxygen consumption rates were also standardized by dry meat mass (V̇_O2,m, in mg O₂ h⁻¹ g⁻¹), specifically to a standard oyster of 1 g dry meat mass. We used the equation V̇_O2,m=(M_std/M_exp)^b×V̇_O2,i, where M_std is the dry meat mass of the standard oyster (1 g), M_exp is the dry meat mass of the experimental oyster, and b is the allometric exponent, 0.58 (Casas et al., 2018).

We examined the effect of population, salinity and the population×salinity interaction on clearance rates (CR_h), 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).

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.

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 log₂+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).

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 10⁶ 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).

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).

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).

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.

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).

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.

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).

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.

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.

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.

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.

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.