DOI: 10.5061/dryad.wpzgmsbmt
title: Data from: Lack of genotype-by-environment interaction suggests limited potential for evolutionary changes in plasticity in the eastern oyster, Crassostrea virginica

description: 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).
description: <b><i>2.1 Oyster conditioning and
 outplanting</i></b> In May 2016, we collected adult
 <i>C. virginica </i>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<sup>th</sup>, 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.
 <b><i>2.2 Sample
 collection</i></b> 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<sup>2</sup> piece of gill
 tissue was dissected in the field from each individual and preserved with
 Invitrogen RNA<i>later</i>. The remaining tissue was blotted
 and placed in a pre-weighed 50 ml test tube to measure wet meat weight.
 <i>P. marinus</i> infection intensities were evaluated by
 adding 0.22 µm filtered seawater (20 ppt) at a concentration of ~0.4 g
 ml<sup>-1</sup> and homogenizing the oyster tissue in each 50
 ml test tube. One ml of the oyster homogenate was used to measure the
 number of <i>P. marinus</i> parasites
 g<sup>-1</sup> oyster wet tissue using the whole-oyster
 procedure (J. F. La Peyre, Casas, & Supan, 2018).
 <b><i>2.3 Morphometrics and infection
 intensity</i></b> 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 &gt;
 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). <i>P. marinus</i> 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). <b><i>2.4 RNAseq
 analysis</i></b> 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 <i>C.
 virginica</i> 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 &gt; 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.<b> The first approach</b> 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<sup>6</sup> 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). <b>The second approach </b>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. <b>The third approach</b> 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. <b><i>2.5 Expression levels, nucleotide diversity, and Ka/Ks Ratios</i></b>             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 <i>C. virginica</i> 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).