Eastern oysters in the northern Gulf of Mexico are facing rapid environmental changes and can respond to this change via plasticity or evolution. Plasticity can act as an immediate buffer against environmental change, but by buffering the effects of selection it could 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 investigate the roles of plastic and evolved responses to environmental variation and Perkinsus marinus infection in Crassostrea virginica by using a common garden experiment with 80 oysters from 6 families outplanted at two field sites naturally differing in salinity. We use a combination of growth data, P. marinus infection intensities, and 3’ RNA sequencing (TagSeq) to identify the effect of genotype (G), environment (E), and genotype-by-environment interaction (GxE) on the oyster’s response to site. As one of first studies to characterize the joint effects of genotype and environment on transcriptomic and physiological profiles in a natural setting, we demonstrate that C. virginica has a highly plastic response to environment and that this response is parallel among genotypes. This suggests that C. virginica may be able to buffer the immediate impacts of future environmental changes by altering gene expression and physiology, although, this may have limited capacity to evolve additional plasticity due to a lack of genetic variation in plasticity (GxE).

2.1 Oyster conditioning and outplanting

In May 2016, we collected adult C. virginica oysters by dredging at three Louisiana, USA estuaries; Vermilion Bay (29° 36’ 39.99’’N, 92° 3’ 19.70’’W), Sister Lake (29° 12’ 50.70’’N, 90° 56’ 3.12’’W), and Calcasieu Lake (29° 51’ 00’’N, 93° 18’ 40’’W). The adults were transported to the Louisiana Department of Wildlife and Fisheries Michael C. Voisin Oyster Hatchery in Grand Isle, LA (29° 14' 20.3" N , 90° 00' 11.2" W) and placed into off-bottom mesh cages for acclimation. In October 2016, after 6 months of acclimation, the oysters were spawned at the MCV oyster hatchery. Individuals were spawned in a 3 male x 2 female cross resulting in 6 full-sibling families. Unfortunately, hatchery records for the population of origin for the parents used in this cross were lost, however, oysters from all three of these collection sites are genetically similar (Johnson & Kelly, 2020). Oyster spat were reared in an upwelling system, individually tagged, and outplanted in one of three adjustable long-line mesh bags at both the Grand Isle Hatchery farm and near the Louisiana Universities Marine Consortium (LUMCON) (29° 15' 12.6" N, 90° 39' 45.9" W) on February 20^th, 2017. These sites naturally differ in salinity regimes, with LUMCON representing a low salinity site and Grand Isle representing a medium-high salinity site. Oysters were assessed for mortality and cleaned of epibionts approximately every 3 months over a 14-month period. Temperature and salinity conditions were measured by the closest USGS monitoring station (USGS 073802516 Barataria Pass at Grand Isle, LA) or by the consortium itself for LUMCON.

2.2 Sample collection

On April 24th, 2018, at the end of the 14-month outplant, 40 of the tagged individuals were sampled at each site. Shell height was measured from shell umbo to distal edge using a digital caliper (Mituyoto USA, Aurora, Illinois). For gene expression analysis, approximately 0.5 cm² piece of gill tissue was dissected in the field from each individual and preserved with Invitrogen RNAlater. The remaining tissue was blotted and placed in a pre-weighed 50 ml test tube to measure wet meat weight. P. marinus infection intensities were evaluated by adding 0.22 µm filtered seawater (20 ppt) at a concentration of ~0.4 g ml^-1 and homogenizing the oyster tissue in each 50 ml test tube. One ml of the oyster homogenate was used to measure the number of P. marinus parasites g^-1 oyster wet tissue using the whole-oyster procedure (J. F. La Peyre, Casas, & Supan, 2018).

2.3 Morphometrics and infection intensity

We compared final height, wet weight, and infection intensities between sites and families. We used the Shapiro-Wilk normality test (Royston, 1982; Shapiro & Wilk, 1965) to determine if the height, weight, and infection intensity data followed a normal distribution. Oyster height and weight data were normally distributed after square-root transformation (Shapiro-Wilk test p > 0.05). The height and weight data were fit with a linear mixed-effect model with bag as a random factor using the R package lme4 (Bates, Mächler, Bolker, & Walker, 2015). A two-way ANOVA was performed using the car package for both height and weight (Fox & Weisberg, 2011). The infection intensity data was not normally distributed and was therefore analyzed using the nonparametrical Kruskal-Wallis Rank Sum test (Hollander, Wolfe, & Chicken, 2013; Kruskal & Wallis, 1952). P. marinus infection intensities differed greatly by outplant site, so to reduce environmental noise we also analyzed infection intensity at each site separately using the Kruskal-Wallis test. To examine pairwise multiple comparisons we used the post-hoc Dunn’s test with the benjamini-hochberg correction (Benjamini & Hochberg, 1995).

2.4 RNAseq analysis

Total RNA was extracted using a E.Z.N.A.® Total RNA Kit I (Omega BIO-TEK Inc., Norcross, GA, USA) following the manufacturer's instructions. The yield and quantity were initially assessed using a NanoDrop 2000 spectrophotometer. Total RNA extracted from the 80 individuals was sent to the University of Texas at Austin’s Genomic Sequencing and Analysis Facility where RNA quality control was confirmed using a 2100 Agilent Bioanalyzer on a Eukaryote Total RNA Nano chip and libraries were produced using the Tag-Sequencing approach (Meyer, Aglyamova, & Matz, 2011). The resulting 80 libraries were sequenced on two lanes of an Illumina HiSeq 2500 platform, with 100 base pair single-end reads. Sequencing reads were trimmed of adapter sequences using Trimmomatic (version 0.39) (Bolger, Lohse, & Usadel, 2014) and base pairs with quality scores below 30 were removed. The trimmed reads were mapped to the C. virginica reference genome (Gómez-Chiarri, Warren, Guo, & Proestou, 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). Reads were mapped to gene features with the options (--quantMode GeneCounts --outFilterScoreMinOverLread 0.50 --outFilterMatchNminOverLread 0.50) specified to adjust for poly-A tail contamination. A count matrix was generated from the ReadsPerGene.out.tab output from STAR. Genotypes for each individual were called using angsd (version 0.931) to produce an identity-by-state matrix. The filters used for assigning IBS scores included removing sites with allele frequency lower than 0.05, requiring a minimum read mapping quality score of 30, a minimum base mapping quality above 20, and removing SNPs with a p-value > 2e-6. Genotype clusters were identified by plotting the first two axes from a distance-based redundancy analysis with the capscale function in the R program vegan (version 2.5-6).

Changes in gene expression were evaluated with three distinct approaches to ensure robustness. The first approach assessed pairwise changes in gene expression associated with genotype (sire and dam) and environment (outplant site) using the package edgeR (version 3.24.2) (M. D. Robinson, McCarthy, & Smyth, 2010). We used the command (filterByExpr) to filter our genes and the remaining read counts were normalized using a trimmed mean of M-values (TMM) normalization method (Mark D Robinson & Oshlack, 2010). Broad changes in gene expression were first assessed using a principal coordinate analysis (PCoA) conducted using the R program vegan with euclidean distances calculated from log2 +1 transformed normalized counts obtained from the cpm() function in edgeR. These log-transformed counts were also used to test for any significant effects of sire, dam, and outplant on gene expression using a PERMANOVA with 10⁶ permutations. The PERMANOVA was performed using the ‘adonis2’ function within the ‘vegan’ package (version 2.5-6). 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 FDR rates calculated using benjamini-hochberg method (Benjamini & Hochberg, 1995). Functional enrichment of differentially expressed genes was tested using a Fisher’s Exact Test. This method calculates enrichment across three gene ontology (GO) categories; Molecular Function (MF), Biological Processes (BP), and Cellular Component (CC).

The second approach identified genes associated with genotype, environment, or GxE interaction using a PERMANOVA. This was performed using the R function ‘adonis2’ within the ‘vegan’ package (version 2.5-6). We used the log-transformed counts (cpm) from edgeR as our count matrix. The PERMANOVA examined the effect of gene expression~sire + dam + outplant + sire*outplant + dam*outplant for each gene independently with 9,999 permutations. The resulting p-values were corrected for multiple comparisons using the benjamini-hochberg method (Benjamini & Hochberg, 1995). To reduce environmental noise, we also ran a PERMANOVA for each site separately to identify genes associated with shell height and infection. This was done since the two sites greatly differ in infection intensities and salinity regime. Functional enrichment of significant genes was tested using a Fisher’s Exact Test.

The third approach used to assess changes in gene expression was a Weighted Gene Co-Expression Network Analysis (WGCNA) (Langfelder & Horvath, 2008). For all WGCNA analyses we restricted the number of genes to remove lowly expressed features retaining only samples with greater than 5 counts per million in 80% of all samples. WGCNA was run using the 11,214 genes that passed this filter, a soft-threshold of 16, a minimum module size of 30, a signed adjacency matrix, and was correlated to sire (a binary variable for sire 1, 2, and 3), dam (a binary variable for dam 1 and 2) and outplant (a binary variable for Grand Isle or LUMCON). To reduce environmental noise, we also applied a WGCNA for height and infection at each site separately. Functional enrichment of significant modules was tested for using a Fisher’s Exact Test.

2.5 Expression levels, nucleotide diversity, and Ka/Ks Ratios

            To test whether differentially expressed genes showed distinctive patterns of nucleotide diversity or evidence of selection on protein coding sequences we used the average log-transformed counts (cpm) from edgeR as our count matrix. Genes responding to genotype, environment, or infection intensity were compared to the global mean for these statistics using the bootstrap version of the Kolmogorov–Smirnov (KS) test (R function ‘ks.boot’ from package Matching version 4.9-7).

We estimated the nonsynonymous to synonymous substitution (Ka/Ks) ratios for each gene using data from an ongoing project where we have collected low-coverage whole genome data for 20 oysters from each of 9 populations (n=180) in the northern Gulf of Mexico spanning from Port Isabel, Texas to Lake Fortuna, Louisiana (Sirovy et al., In prep) (Table S1). DNA libraries were sequenced on an Illumina HiSeqX platform, with 150 bp paired-end reads. Sequencing reads were trimmed of adapter sequences using Trimmomatic (version 0.39) (Bolger, Lohse, & Usadel, 2014) and base pairs with quality scores below 30 were removed. The trimmed reads were mapped to the C. virginica reference genome (Gómez-Chiarri, Warren, Guo, & Proestou, 2015) using BWA (Burrow’s Wheeler Aligner) (Li & Durbin, 2009) followed by removal of duplicates and indexing using Picard tools software (https://github.com/broadinstitute/picard). We estimated the number of SNP’s using the GATK pipeline (version 4.0.2.1) and annotated those SNP’s using SnpEff (Cingolani et al., 2012). Ka/Ks ratios were then calculated based on the VCF outputs from SnpEff using a script available online at (https://github.com/MerrimanLab/selectionTools/blob/master/extrascripts/kaks.py). Ka/Ks scores were filtered to keep values above 0 and below 3 to limit bias from outliers and significant contrasts were tested down to a ratio of 1.

We generated the nucleotide diversity (Pi) for the coding + intron, upstream, and downstream region using reduced representation bisulfite sequencing data for the same 80 oysters (Johnson et al., In prep). SNP’s were called using program BS-SNPer (Gao et al., 2015) which outputs a VCF file. The VCF’s were normalized, indexed, and converted to BCF’s using BCFtools (version-1.10.2). Nucleotide diversity was then calculated using the program VCFtools (version-0.1.14) with 1000 bp windows and a step size of 1000 bp.

To assess the nucleotide diversity and Ka/Ks levels for the DEGs, we used a series of permutation tests to compare the mean of the DEGs against the genome wide average. For this analysis, the mean nucleotide diversity (Pi) or Ka/Ks ratio was calculated for DEGs from any of the three main approaches. Then, using the same number of significant genes, the nucleotide diversity and Ka/Ks values were shuffled in R using the sample function to generate a new difference in means and this was repeated 10,000 times to generate a null distribution. This allowed us to test if the mean difference in the DE genes was significantly different than the mean difference in the permuted list. Finally, by taking the mean of the absolute value of the permuted list we were able to calculate a p-value based on the number of times the permuted value was greater/less than the observed value.

Nucleotide diversity and Ka/Ks ratios were also compared to average expression levels which had been binned using the R function ‘decile’ within the package StatMeasures (version 1.0). A Kruskal-Wallis Test was performed comparing the nucleotide diversity or Ka/Ks values with the 10 gene expression bins, followed by a Dunn’s test with the benjamini-hochberg correction (Benjamini & Hochberg, 1995).