DOI: 10.5061/dryad.wh70rxwq9
title: Data From: Population-specific responses in eastern oysters exposed to low salinity in the northern Gulf of Mexico

description: Eastern oysters, Crassostrea virginica, are facing rapid environmental
 changes in the northern Gulf of Mexico and can respond to these changes
 via plasticity or evolution. Plastic responses can immediately buffer
 against environmental changes, although this buffering may impact the
 organism’s ability to evolve in subsequent generations. While plasticity
 and evolution are not mutually exclusive, the relative contribution and
 interaction between them remain unclear. In this study, we investigate the
 roles of plastic and evolved responses to low salinity in C. virginica
 using a common garden experiment with four populations exposed to two
 salinities. We use three transcriptomic analyses (edgeR, PERMANOVA, and
 WGCNA) combined with physiology data to identify the effect of genotype
 (population), environment (salinity), and genotype-by-environment
 interaction on both whole organism and molecular phenotypes. We
 demonstrate that variation in gene expression is mainly driven by
 population, with relatively small changes in response to salinity. In
 contrast, the morphology and physiology data reveal that salinity has a
 larger influence on oyster performance than the population of origin. All
 analyses lacked signatures of genotype-by-environment interaction, and in
 contrast to previous studies, we find no evidence for population-specific
 responses to low salinity. However, individuals from the highest salinity
 estuary displayed highly divergent gene expression from other populations,
 which could potentially drive population-specific responses to other
 stressors. Our findings suggest that C. virginica largely rely on
 plasticity in physiology to buffer the effects of low salinity, but that
 these changes in physiology do not rely on large persistent changes in
 gene expression.
description: <strong>Oyster collection and
 conditioning</strong> In December 2017 and
 January 2018, subtidal adult <em>C. virginica</em> oysters
 were collected from two estuaries in Louisiana, USA; Calcasieu Lake (CL;
 29° 50′ 58′′ N, 93° 17′ 1′′ W) and Vermilion Bay (VB; 29° 34′ 47′′ N, 92°
 2′ 4′′ W). Additionally, in August 2018 oysters were collected from two
 estuaries in Texas, USA; Packery Channel (PC; 27° 37′ 38′′ N, 97° 13′ 59′′
 W) and Aransas Bay (AB; 28° 7′ 38′′ N, 96° 59′ 8′′ W). In August 2018, all
 four broodstock populations were naturally induced to spawn at the Auburn
 University Shellfish Laboratory (AUSL) in Dauphin Island, Alabama, as
 described in Marshall et al. (2021). Oyster spat were maintained in
 upwelling nursery systems until ~6mm in shell height, and subsequently
 deployed in bags at an AUSL-permitted grow-out site at Bayou Sullivan,
 Alabama (30° 21′ 52′′ N, 88° 12′ 57′′ W) before being moved in March 2019
 to the Grand Bay Oyster Park (GBOP), Alabama (30˚ 22′ 15′′ N, 88˚ 19’ 0”
 W) for further growth (see Marshall et al., 2021 for environmental data).
 In February 2020, we transferred the oysters (~18 months old) from each
 stock from GBOP to the static systems at Louisiana State University
 Agricultural Center Animal and Food Sciences laboratory building (AFL) in
 Baton Rouge, Louisiana. Oysters were scrubbed, tagged,
 and placed into one of two 400-L tanks with recirculating artificial
 seawater (Crystal Sea Marinemix, Marine Enterprises International,
 Baltimore, MD, USA) adjusted to a salinity of 20 ppt and temperature of
 20°C. Each population had ~10 oysters in each of the two tanks. The
 temperature was gradually adjusted at a rate of 1.0°C every 2 days until
 the experimental temperature of 25°C was reached. Salinity in each tank
 was gradually adjusted at a rate of 3 units every 2–3 days until the
 target salinities were reached (6ppt or 12ppt). We used artificial
 seawater or carbon-filtered and aerated freshwater to adjust the salinity.
 Oysters were fed daily at ~5% of their meat dry weight with Shellfish Diet
 1800® (Reed Mariculture Inc, Campbell, CA, USA). After
 3 weeks of acclimation at each salinity, we began measuring clearance
 rates and oxygen consumption rates. We had limited access to the
 laboratory due to COVID restrictions, which increased the time oysters
 spent in these tanks and increased overall mortality. This increased
 mortality led to an uneven sampling between the populations and salinities
 (Table S1). In July 2020 we measured the final shell height, shell width,
 shell weight, and dry meat weight of the surviving oysters for each
 population at each salinity (Table S1). Shell height was measured from
 shell umbo to distal edge and shell width was measured as the maximum
 distance between the two valves when closed (Galtsoff, 1964). Dry meat
 weight was measured after drying at 70 °C for 48 h. For gene expression
 analysis, we sampled 5 oysters per population per salinity (n=40).
 Approximately 0.5 cm<sup>2</sup> piece of gill tissue was
 dissected and immediately frozen in liquid nitrogen and placed in a -80˚C
 freezer until later extraction. <strong>Measuring
 clearance rates</strong> The clearance rate,
 defined as the volume of water completely cleared of suspended particles
 per unit of time, was quantified using a static system (Widdows, 1985).
 Oysters were individually placed in 1 L beakers filled with 1 L of 0.5-μm
 filtered seawater and left undisturbed for ~ 60 min. Shellfish Diet 1800®
 was added to each beaker to give an initial concentration of 3 ×
 10<sup>4</sup> cells mL<sup>−1</sup>. Particle
 concentrations (&gt;5 μm) were collected every 30 seconds until
 particle count declined to 50% of original values and was measured using a
 PAMAS Model S4031GO particle counter (PAMAS Partikelmess-und
 Analysesysteme GMBH, Rutesheim, Germany). We used air stones to reduce
 algal sedimentation. Beakers with algal cells and oyster shells inside
 were used as controls. Individual oyster clearance rate
 (CRi) was calculated according to the equation CRi =
 [(<em>b</em> – <em>b'</em>) × Vol(L) × 60 min
 h<sup>−1</sup>], where <em>b </em>is the slope of
 the linear regression between the natural logarithm of cell concentration
 (cells mL<sup>−1</sup>) and time (min) in a beaker,
 <em>b′</em> is the slope for the control beaker, and Vol is
 the volume of seawater in the beaker (1 L). This method follows Cranford
 et al. (2016), in which a generalization of Coughlan’s (1969) method is
 developed to integrate the values of intermediate particle concentration
 samplings, which aims to reduce the impact of outliers. Accordingly, only
 linear regressions with an <em>r<sup>2</sup></em>
 &gt; 0.90 were included in these calculations. CRi measurements are in
 units of L h<sup>−1</sup>. 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
 =(?<sub>std</sub>/?<sub>exp</sub>)<sup>?</sup>×CRi, where CRh is the clearance rate standardized by shell height, Hstd is the shell height of the standard oyster (80 mm), Hexp is the shell height of the experimental oyster, CRi is the individual oyster clearance rate, and <em>b </em>is the allometric exponent for shell height set to 1.78 (Cranford et al. 2011). CRh measurements are in units of L h<sup>−1</sup> 80 mm<sup>−1</sup>. <strong>Measuring oxygen consumption rates</strong> To measure the oxygen consumption rate of fed oysters, known as the routine oxygen consumption rate (Bayne, 2017; Gosling, 2015), individual oysters were placed in 915-mL acrylic chambers filled with 0.5-μm filtered seawater from their respective tank. We used a ProOBOD probe (YSI Incorporated, Yellow Springs, OH, USA) to measure dissolved oxygen and temperature recorded every 30 s by a YSI Multilab 4010-3 m (YSI Inc.) connected to the probe. The readings started once the oyster’s valves opened and ran continuously for up to 90 min or until dissolved oxygen fell below 70% saturation. Individual oyster oxygen consumption rate (OCRi) was calculated as OCRi = [(<em>b</em> – <em>b'</em>) × Vol (L) × 60 min h<sup>−1</sup>], where b is the slope of the linear regression of oxygen concentration (mg L−1) versus time, <em>b′</em> is the slope for the control, and Vol is the volume of water (in L) in the chamber (chamber minus oyster volume). OCRi measurements are in units of mg O<sub>2</sub> h<sup>−1</sup>. Only linear regressions with an <em>r<sup>2</sup></em> &gt; 0.95 were included in the calculations. Oxygen consumption rates were also standardized by dry meat weight (OCRw), specifically to a standard oyster of 1 g dry meat weight. We used the equation OCRw =(?<sub>std</sub>/?<sub>exp</sub>)? ×OCRi, where W<sub>std</sub> is the dry meat weight of the standard oyster (1 g), W<sub>exp</sub> is the dry meat weight of the experimental oyster, and <em>b</em> is the allometric exponent, 0.58 (Casas et al., 2018). OCRw measurements are in units of mg O<sub>2</sub> h<sup>−1</sup> g<sup>−1</sup>. <strong>Statistical analyses</strong> We examined the effect of population, salinity, and population x salinity interaction on clearance rates (CRh), oxygen consumption rates (OCRw), and dry meat weight. We used a Shapiro-Wilk normality test (Shapiro & Wilk, 1965; Royston, 1982) to determine if the data followed a normal distribution and used either a log or square root transformation if needed. A two-way ANOVA was performed using the R function “aov” in the <em>stats </em>package (version 3.6.2) (p-value &lt; 0.05). To examine pairwise multiple comparisons we used the post-hoc Dunn’s test with the Benjamini-Hochberg correction (p-value &lt; 0.05) (Benjamini & Hochberg, 1995). <strong>Gene expression analysis; sampling and initial processing</strong> 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 <em>C. virginica</em> reference genome (Gómez-Chiarri et al., 2015) with known haplotigs removed (https://github.com/jpuritz/OysterGenomeProject/tree/master/Haplotig_Masked_Genome) using the single pass option for STAR RNA-seq aligner (version 2.6.0a) (Dobin et al., 2013). Three samples displayed poor sequencing quality (i.e., the lowest number of reads, poorest mapping rates, and isolated in PCA) and were subsequently removed for all downstream analyses. A count matrix was generated from the ReadsPerGene.out.tab output from STAR. To ensure statistical robustness, we used three different approaches to examine gene expression: edgeR, a PERMANOVA test, and a Weighted Gene Co-Expression Network Analysis (WGCNA). <strong>Gene expression analysis; edgeR</strong> We first assessed pairwise changes in gene expression associated with population and salinity using the package edgeR (version 3.5.2) (Robinson et al., 2010). We filtered all genes using the command (filterByExpr) and the remaining read counts were normalized using a trimmed mean of M-values (TMM) normalization method (Robinson & Oshlack, 2010). Broad changes in gene expression were first assessed using a principal component analysis (PCA) conducted using the R program vegan with euclidean distances calculated from log2 +1 transformed normalized counts obtained from the cpm() function in edgeR. Pairwise differential expression was measured using the genewise negative binomial generalized linear model implemented in the edgeR function glmQLFit, and significantly differentially expressed genes (DEGs) were identified based on False Discovery Rates (FDR) calculated using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995). Functional enrichment of differentially expressed genes was tested using a Fisher’s Exact Test (p-value &lt; 0.05). All Fisher’s Exact Tests were conducted using the scripts originally developed by Wright et al., 2015 and available at: (https://github.com/z0on/GO_MWU/blob/master/GO_MWU.R). This method uses a binary input (DEGs = 1, non-DEGs = 0) to calculate if there is enrichment in the DEGs across three gene ontology (GO) categories: Molecular Function (MF), Biological Processes (BP), and Cellular Component (CC). <strong>Gene expression analysis; PERMANOVA</strong> We used a permutational multivariate analysis of variance test (PERMANOVA) to identify genes associated with population, salinity, or population x salinity interaction. We performed two PERMANOVA analyses to examine the effect on total gene expression and on a per gene basis. For both PERMANOVA analyses, log-transformed counts were evaluated using the R function ‘adonis2’ within the ‘vegan’ package (version 2.5-6). The PERMANOVA examined the effect of gene expression~ population + salinity + population*salinity with 10<sup>6</sup> permutations (total expression) or 99,999 permutations (per-gene expression) (p-value &lt; 0.05). The difference in the number of permutations between the total expression and per-gene analysis is the result of the higher computational requirements for the per-gene analysis. For the per-gene PERMANOVA, the p-values were corrected for multiple comparisons using the Benjamini-Hochberg method (Benjamini & Hochberg, 1995). Functional enrichment of significant genes was tested using a Fisher’s Exact Test (p-value &lt; 0.05). <strong>Gene expression analysis; WGCNA</strong> We used a Weighted Gene Co-Expression Network Analysis (WGCNA) (Langfelder & Horvath, 2008) to find clusters of genes with highly correlated expression patterns associated with population, salinity, oxygen consumption rate, final shell height, shell width, and dry meat weight (n = 37). We also ran a separate WGCNA analysis for individuals with clearance rate measured (n=23). For all WGCNA analyses, we restricted the number of genes to remove lowly expressed features by only using samples with greater than 5 counts per million in 80% of all samples. WGCNA was run using the 7,575 genes (7,533 genes for clearance rate) that passed this filter, a soft-threshold of 16, and a minimum module size of 30. Functional enrichment of significant modules was tested using a Fisher’s Exact Test (p-value &lt; 0.05).