Characterizing the Epigenetic and Transcriptomic Responses to Perkinsus marinus Infection in the Eastern Oyster Crassostrea virginica

Eastern oysters in the northern Gulf of Mexico are routinely infected with the protistan parasite Perkinsus marinus, the cause of the disease commonly known as dermo. Recent experimental challenges among Atlantic coast populations have identified both resistant and susceptible genotypes using comparative transcriptomics. While controlled experimental challenges are essential first assessments, expanding this analysis to field reared individuals provides an opportunity to identify key genomic signatures of infection that appear both in the laboratory and in the field. In this study we combined reduced representation bisulfite sequencing with 3′ RNA sequencing (Tag-seq) to describe two molecular phenotypes associated with infection in oysters outplanted at a common garden field site. These combined approaches allowed us to examine changes in DNA methylation and gene expression for a large number of individuals (n = 40) that developed infections during the course of a common garden outplant experiment. Our epigenetic analysis of DNA methylation identified significant changes in gene body methylation associated with increasing infection intensity, across genes associated with immune responses. There was a smaller transcriptomic response to increasing infection intensities with 32 genes showing differential expression; however, only 40% of these genes were found to also be differentially methylated. While there was no clear pattern between direction of differential methylation and gene expression, there was a significant effect of percent methylation on gene-by-gene expression levels and the coefficient of variation in gene body methylation between treatments. These results show that in C. virginica, heavily methylated genes have high levels of gene expression with low levels of variation. Comparing our differential expression results with previously published experimental P. marinus challenges identified overlapping expression patterns for genes associated with C1q-domain-containing and V-type proton ATPase proteins. Through our comparative transcriptomic approach using field reared individuals and co-expression network analysis we have also been able to identify a network of genes that change in expression in response to infection. These combined analyses provide evidence for a conserved response to P. marinus infections across infection intensities and suggest that DNA methylation may not be a reliable predictor of differential gene expression in long-term infections.

Eastern oysters in the northern Gulf of Mexico are routinely infected with the protistan parasite Perkinsus marinus, the cause of the disease commonly known as dermo. Recent experimental challenges among Atlantic coast populations have identified both resistant and susceptible genotypes using comparative transcriptomics. While controlled experimental challenges are essential first assessments, expanding this analysis to field reared individuals provides an opportunity to identify key genomic signatures of infection that appear both in the laboratory and in the field. In this study we combined reduced representation bisulfite sequencing with 3 RNA sequencing (Tag-seq) to describe two molecular phenotypes associated with infection in oysters outplanted at a common garden field site. These combined approaches allowed us to examine changes in DNA methylation and gene expression for a large number of individuals (n = 40) that developed infections during the course of a common garden outplant experiment. Our epigenetic analysis of DNA methylation identified significant changes in gene body methylation associated with increasing infection intensity, across genes associated with immune responses. There was a smaller transcriptomic response to increasing infection intensities with 32 genes showing differential expression; however, only 40% of these genes were found to also be differentially methylated. While there was no clear pattern between direction of differential methylation and gene expression, there was a significant effect of percent methylation on gene-by-gene expression levels and the coefficient of variation in gene body methylation between treatments. These results show that in C. virginica, heavily methylated genes have high levels of gene expression with low levels of variation. Comparing our differential expression results with previously published experimental P. marinus challenges identified overlapping expression patterns for genes associated with C1q-domain-containing and V-type proton ATPase proteins. Through our comparative transcriptomic approach using field reared individuals and co-expression network analysis we have also been able to identify a network of genes

INTRODUCTION
The relationship between tissue-specific DNA methylation and gene expression is an area of interest to a broad research community as these two molecular phenotypes have been proposed as mechanisms of environmental responsiveness and used as biomarkers to describe interactions between organisms and their environment (Roberts and Gavery, 2012;Dixon et al., 2014;Putnam et al., 2016;Hawes et al., 2018;Schmid et al., 2018;Eirin-Lopez and Putnam, 2019;Rubi et al., 2019). The utility of these metrics relies heavily on a well-developed understanding of how a stimulus will modulate either gene expression or DNA methylation levels, and how methylation could in turn modulate gene expression. As a result, the last decade has seen a dramatic increase in the number of studies using whole transcriptome sequencing to describe how non-model species modify their transcriptomes in response to developmental cues, diseases, and shifts in the abiotic environment (Alvarez et al., 2015;De Wit et al., 2018;Jones et al., 2019;Proestou and Sullivan, 2020). Interest in the functional role of DNA methylation has also grown in recent years due to its potential contributions to transgenerational plasticity (Kappeler and Meaney, 2010;Herman and Sultan, 2016;Gavery and Roberts, 2017;Ryu et al., 2018). These studies together have further fueled interest in the mechanisms that control transcriptome variation, with the majority of studies focusing on either DNA methylation or chromatin modifications (Clark et al., 2018;Weinhold, 2018;Eirin-Lopez and Putnam, 2019).
Oysters in the genus Crassostrea are both an ecologically and environmentally important keystone species distributed across broad geographic ranges. These widely dispersed species inhabit highly variable environments leading to species with a high degree of phenotypic plasticity . As such, oysters are an ideal organism for exploring how DNA methylation influences phenotypic plasticity. Recent research on this genus has focused on describing the degree to which gene expression contributes to phenotype or how DNA methylation directs gene expression levels. These studies have found strong evidence that environmental variation and/or disease state can have a predictable influence on gene expression (Yan et al., 2017;Jones et al., 2019;Proestou and Sullivan, 2020). In a similar manner, a growing body of literature has provided ample evidence that DNA methylation can shape the transcriptomic phenotype (Rivière et al., 2013;Olson and Roberts, 2014;Song et al., 2017). In the Pacific oyster (C. gigas) high levels of methylation in gene bodies are associated with elevated expression, while lowly methylated genes showed higher levels of transcriptomic plasticity (Gavery and Roberts, 2013;Olson and Roberts, 2014).
In the eastern oyster (C. virginica) DNA methylation is more divergent between populations than single nucleotide polymorphism-based estimates of divergence, suggesting a nongenetic (i.e., environmental) influence on DNA methylation (Johnson and Kelly, 2020). Together, these studies provide some evidence of an association between DNA methylation and transcriptomic expression levels; however, the functional consequences of changes in DNA methylation are still largely unknown in bivalves.
In a similar manner, shifts in gene expression are often essential for responding to infectious diseases in bivalves (Renault et al., 2011;Moreira et al., 2012;Rosani et al., 2015). Previous research investigating the transcriptomic responses of C. virginica challenged with Perkinsus marinus found significant changes in the regulation of genes involved in immune defense (Tanguy et al., 2004;Wang et al., 2010). Although previous studies have shown that C. virginica regulates gene expression in response to P. marinus, a major gap remains regarding whether there is a population-specific transcriptomic response to dermo in the northern Gulf of Mexico. Previous research has shown that oyster populations vary in mortality rates in response to dermo, though the molecular basis for this is still unclear Leonhardt et al., 2017;La Peyre et al., 2019). A recent study looking at the effect of P. marinus infection on the global gene expression pattern of resistant vs susceptible C. virginica families found strong acute responses to infection with over 3,000 differentially expressed transcripts with resistant oysters upregulating genes involved with peptidase inhibitor activity and regulation of proteolysis (Proestou and Sullivan, 2020). This response diminished over 28 days after which only 21 differentially expressed genes were observed within the dermo resistant family while the dermo susceptible family differentially expressed over 2,000 genes. These large differences in response to infection in controlled settings provide a base to begin exploring how responses to infections differ in the field.
To this end we present a comparison of 40 individuals with differing levels of P. marinus infections following 14 months of outplant at a common-garden field site. We investigated both transcriptomic and epigenetic signatures of these infections to better understand how light vs heavy intensity infections influence molecular phenotypes. Using a combination of reduced representation bisulfite sequencing and 3 RNA sequencing (TAGseq) we are able to describe two molecular responses to infection intensity in a common garden study. These results also allow us to compare the molecular phenotype of farmraised oysters with P. marinus infections to recent results from a controlled infection (Proestou and Sullivan, 2020). Through a common garden outplant of hatchery reared juveniles we were able to compare changes in DNA methylation and gene expression for a large number of individuals (n = 40) that were identified as infected with P. marinus, but at varying intensities. This allowed us to use a weighted gene co-expression network analysis to identify a collection of genes positively correlated with infection. Furthermore, we were able to explore how changes in DNA methylation related to changes in gene expression on a gene-by-gene basis.

Oysters
In May 2016, adult oysters (C. virginica) were collected by dredging from each of two estuaries; Vermilion Bay, LA (29 • 35 7.26 N, 91 • 0\2 33.92 W) and Calcasieu Lake, LA (29 • 50 58.02 N, 93 • 17 1.32 W). These oysters 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 common garden acclimation. In October 2016, after 5 months of acclimation, oysters were spawned at the MCV oyster hatchery. Oyster spat were reared in an upwelling system, individually tagged, and outplanted in one of three adjustable long-line mesh bags at the Grand Isle Hatchery farm on February 20, 2017. Oysters within each bag were monitored for mortality and cleaned of epibionts approximately every 3 months over a 14-month period.

Sample Collection
On April 24, 2018, after 14 months at the Grand Isle outplant site, 50 individuals were haphazardly sampled. Shell height of each individual was measured from shell umbo to distal edge using a digital caliper (ABS Coolant Proof Calipers, Mituyoto Corporation, Japan). Approximately 1 cm 2 piece of gill tissue was sampled in the field from each individual and preserved with either Invitrogen RNAlater (gene expression) or 95% ethanol (DNA methylation). The remaining whole animal was placed in a pre-weighed 50 ml test tube in order to measure wet meat weight. Infection intensities were enumerated by adding 0.22 µm filtered seawater (20 ppt) at a concentration of ∼0.4 g ml −1 and homogenizing the oyster meat in each 50 ml test tube. One ml of the oyster homogenate was used to measure the number of P. marinus hypnospores g −1 oyster wet tissue using the wholeoyster procedure (La Peyre et al., 2018). Oyster infections were classified as either very-light (≤1,000 parasites g −1 wet tissue), light (1,001-10,000 parasites g −1 wet tissue), moderate (10,001-100,000 parasites g −1 wet tissue) or moderate-heavy (>100,000 parasites g −1 wet tissue). Of the 50 individuals sampled we sequenced 40 individuals that spanned the four different infection categories ( Table 1).

DNA Methylation Analysis
DNA was extracted using the OMEGA E.Z.N.A. Tissue DNA Kit (D3396-01; Omega bio-tek) with a 2 min RNase A digestion to remove co-purified RNA. DNA purity was assessed based on 260/280 and 260/230 ratios using a nanodrop spectrophotometer (ND1000; Thermofisher Scientific). Presence of high molecular weight DNA was confirmed using a 1.5% agarose gel, and DNA concentration was verified using a Qubit 3.0 Fluorometric dsDNA BR assay kit (Q32850; Life Technologies). The epiGBS library preparation followed previously published methods (Van Gurp et al., 2016;Johnson and Kelly, 2020). Briefly, a total of 500 ng of purified genomic DNA was double digested using the two frequent cutter enzymes AseI and NsiI (NEB-R0127L and NEB-R0526L; Van Gurp 2016). Digested DNA was ligated to custom y-yoked methylated sequencing adapters using a T4 DNA ligase (B9000S; New England Biolabs) with additional rATP to ensure ligation of custom adapters (Glenn et al., 2019). The adapter ligated DNA was bisulfite converted in a 96 well plate using the Zymo Research EZ DNA Methylation-Lightning kit (D5031; Zymo Research) with a 15 min L-desulphonation step. This bisulfite converted DNA was tagged and amplified with Illumina adapters using 16 cycles of PCR. Amplified libraries were size selected to 300-600 base-pairs (bp) using the Zymo Research Select-A-Size DNA clean & concentrator (D4080; Zymo Research). Size selection was confirmed using the Agilent Bioanalyzer DNA high sensitivity chip (5067-4626; Agilent Technologies). Libraries were pooled and sequenced by NovoGene Inc. (R) with a 10% PhiX spike-in on a full flow cell of the Illumina HiseqX with 100 bp paired-end reads.
The epiGBS sequencing reads were adapter trimmed and base pairs with a phred score less than 30 were removed using Trimmomatic (version 0.39) (Bolger et al., 2014). Trimmed reads were mapped to the reference genome (NCBI GCF_002022765.2) and CpG methylation was called using the software package bismark (v0.19.0) (Krueger and Andrews, 2011). The bismark commands used in the mapping allowed for 1 mismatch in a seed alignment of 10 with a minimum alignment score setting of −0.6 (-score_min L, 0, −0.6). These settings having previously been used in this species and were selected to account for genomic variations between C. virginica collected from the northern Gulf of Mexico (nGOM, this study) and the disease-resistant inbred line from the United States East Coast used for the construction of the reference genome (Gómez-Chiarri et al., 2015;Johnson and Kelly, 2020). CpG methylation was extracted from the non-deduplicated mapped reads using the bismark command bismark_methylation_extractor with the following commands; -ignore_r2 2, -bedGraph, -zero_based, -no_overlap, -cytosine_report, and -report. Differential methylation was conducted on CpG features using the bismark coverage files with methylation across both strands merged and analyzed using the R program MethylKit (v.1.2.4) (Akalin et al., 2012). Methylated regions were identified using a tiled window approach with a tile size of 1,000 bp (1 kb) and a step size of 1 kb. The 1 kb regions were filtered using the filterByCoverage command to require coverage greater than 10× in at least 8 of the 20 individuals.
Pair-wise differential methylation was measured between individuals with moderate or moderate-heavy (parasites g −1 >10,000, n = 15) and very-light (parasites g −1 ≤1,000, n = 16) infection intensities. For these analysis, one sample in the high infection group was not included as a result of poor sequence quality. Significantly differentially methylated regions (DMRs) Mean shell height and mean wet meat weight are also reported with standard deviations.
were identified based on a minimum percent methylation difference between groups of 20% and an adjusted P-value (qvalue) less than or equal to 0.05 (Mathers et al., 2019). Functional enrichment of methylation was assessed using two methods. The first approach explored enrichment in methylation amongst all samples tested using a Mann-Whitney U-test to identify enriched ontologies between highly and lowly methylated genes. This test used mean methylation among all individuals for each gene feature with at least a single 1 kb region overlapping a gene's promoter, gene body or downstream region; and, calculated enrichment across each gene ontology (GO) category (MF, Molecular Function; BP, Biological Process, and CC, Cellular Component). The second approach tested for functional enrichment between genes that showed differentially methylated regions using a Fisher's Exact test for each GO category. For both analyses we used a background list consisted of all genes with GO annotations for which DNA methylation was measured (n = 8,624 of genes passing filter with GO annotations).

TAGseq Analysis
Total RNA was extracted using a E.Z.N.A. R Total RNA Kit I (Omega BIO-TEK Inc., Norcross, GA, United States) 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 control was confirmed using a 2100 Agilent Bioanalyzer on a Eukaryote Total RNA Nano chip and libraries were produced using the 3 poly-A-directed mRNAsequencing (TAGseq) method (Meyer et al., 2011). The resulting 40 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 et al., 2014) and base pairs with quality scores below 30 were removed. The trimmed reads were then mapped to the published C. virginica reference genome (Gómez-Chiarri et al., 2015) using the single pass option for STAR RNA-seq aligner (version 2.6.0a) (Dobin et al., 2012). We used the default of 10 allowed mismatches for filtering and allowed for multi-mapping. After mapping, we used HTSeq (version 0.11.2) (Anders et al., 2015) to obtain the number of reads mapped to each gene and allowed for multiple alignments to be counted (-non-unique = all). The counts were sorted by alignment position and based on gene features obtained from the C. virginica genome assembly on NCBI (version GCF_002022765.2).
Changes in gene expression associated with P. marinus infection intensity were tested with two approaches. The first approach assessed pairwise changes in gene expression between individuals with moderate-heavy and very-light infection using the package edgeR (version 3.24.2) . For this analysis genes with fewer than three counts per million mapped reads across 50% (n = 20) of all samples were removed. The remaining read counts were distributed across 17,439 gene features and were normalized using the trimmed mean of M-values (TMM) normalization method (Robinson and 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 interactions between gene expression and infection intensity interactions using a PERMANOVA with 1e 6 permutations. Pairwise differential expression was measured using the genewise negative binomial generalized linear model implemented in the edgeR function glmFit and significantly differentially expressed genes (DEGs) were identified based on FDR rates calculated using benjamini-hochberg method (Benjamini and Hochberg, 1995). Comparisons were made between the lowest infection group (very-light, n = 16) and a combination of the two highest infection groups (moderate and moderate-high, n = 16). Functional enrichment of differential gene expression for each comparison was tested using a rankbased gene ontology analysis with adaptive clustering that uses a Mann-Whitney U-Test to identify enriched ontologies (Wright et al., 2015). For these tests, we tested for enrichment using the log-transformed and signed p-value for each gene in each comparison contrasted using all of the genes that passed the expression filtering for the analysis and had annotated GO terms (n = 11,057).
The second approach used to assess changes in gene expression associated with P. marinus infection intensity was a Weighted Gene Network Analysis (WGCNA). For this component of the analysis we restricted the number of genes to remove lowly expressed features retaining only samples with greater than five counts per million in 75% of all samples (n = 30). This additional filtering was included to remove genes with low counts prior to network analysis. WGCNA was run using the 11,998 genes that passed this filter, a soft-threshold of 12, a minimum module size of 30, a signed adjacency matrix, and was correlated to shell height, meat weight, and infection intensity (a continuous variable of counts/g). Functional enrichment was assessed across 7,398 of the 11,998 genes that had annotated GO terms for the five modules that showed significant correlation with any of the five traits using the GO_MWU methods for WGCNA using eigengene-based module connectivity (kME) as a continuous variable. This method calculates significant enrichment across each gene ontology (GO) category (MF, BP, and CC) for each module, and identifies functional categories in the modules using the Fisher's Exact test. These functional categories were further tested for an association with higher kME values using a within-module MWU test. The product of the resulting p-values from the Fisher's Exact and MWU test was then used to calculate false discovery rates with ten permutations of randomly shuffled significance measures among genes.

Methylation and Gene Expression
An important component of this study was the opportunity to examine the interaction between DNA methylation and transcript expression across many individuals. We tested the relationship between percent methylation for moderate-heavy infected individuals vs very-lightly infected individuals and changes in gene expression, as measured by the pairwise differential gene expression analysis. To do this, we summarized each 1 kb region into three categories: 1 kb regions that were located in (i) gene promoter regions (within 2 kb upstream of the first exon), (ii) gene body regions (first exon to last exon), or (iii) downstream regions (within 2 kb downstream of the last exon). For all comparisons if multiple 1 kb regions were present within the annotated feature (promoter, gene body, or downstream region) the mean p-value, mean q-value, and mean percent difference in methylation was calculated. For each of these three gene regions we explored the association of mean methylation with the log-transformed cpm expression data. For these analyses cpm expression data was distributed into 10 deciles using the decile function available in the R package StatMeasures (v.1.0). Significance of the association between mean methylation and expression level deciles was assessed using a Kruskal-Wallis test along with a Dunn post hoc test and corrected for multiple comparisons using the Benjamini-Hochberg method.

DNA Methylation
RRBS sequencing produced a total of 878 million reads with an average of 21.4 million reads per sample. Trimming of these reads led to an average of 18.1 million reads per sample, and mapping resulted in an average of 84.6% of reads mapping to the reference genome (range: 77.3-89.1%). Genome tiling in methylkit distributed these reads across 43,085 1-kb regions with a minimum count of 10 in a minimum of 9 individuals in either of the 2 groups (very-light or moderate and moderate-heavy). Of these, 10,741 regions were located along intergenic regions (i.e., greater than 2 kb from any gene feature). The remaining 32,344 1kb regions were distributed across 15,543 annotated gene features with 13.1% (n = 2,364 genes) found within promoter regions (i.e., 2 kb upstream of first exon), 71.3% (n = 12,781 genes) found within gene bodies, and 15.6% (n = 2,791 genes) found within 2 kb downstream of a gene's final exon. Functional enrichment testing using the rank-based gene ontology analysis with adaptive clustering among all 40 individuals based on mean methylation identified enriched gene functions among both heavily-and lightly methylated genes (relative to the genome-wide average) in each of the three gene ontology categories (14 MF, 6 BP, and 11 CC). We observed higher methylation of key molecular processes such as RNA metabolism, chromosome organization, and general protein binding. While lightly methylated genes were associated with environmentally responsive categories such as enzyme regulation, immune processes and oxidoreductase activity ( Table 2).

Gene Expression
Transcriptome sequencing using TAGseq produced a total of 189 million reads, with 4.7 million reads per sample. Trimming of those reads led to a final read count of 4.5 million per sample, which is sufficient for this method (Meyer et al., 2011). Star mapping resulted in 80.20% of reads mapping to the reference genome (range: 75.7-83.1%). Filtering mapped reads based on expression resulted in a total of 17,439 gene features with more than 3 counts per million reads in 20 of the 40 samples. This stringent level of filtering was chosen to avoid including too many genes with low count totals. Principle coordinate analysis revealed significant overlap with the three low-infection groups with separation of these from the single high infection group (Figure 1). The PERMANOVA testing expression ∼ infection (parasites g −1 wet tissue wet) found a non-significant interaction between infection and gene expression (p-value = 0.09).
Pairwise differential gene expression analysis identified modest levels of differential gene expression between infection intensities. Gene expression changes associated with moderateheavy vs very-light infections, identified 31 genes up-regulated and 8 genes down-regulated in the moderate-heavy infected individuals (Figure 2). A Mann-Whitney U-Test of the 11,057 genes in the analysis (regardless of DE status) was run Percent of reference shows the number of genes identified as enriched when compared to the total number of genes for that GO term in the background reference list.
FIGURE 1 | Biplot of the first two principle components from principle coordinate analysis using the gene expression data.

WGCNA Results
The WGCNA analysis identified a total of 12 gene modules; of these, only 1 module was significantly correlated with dermo infection intensity ( Figure 3A). Below we discuss the significance of module pink as this module was the only one significantly associated with infection.

Module Pink
Module Pink was comprised of 385 transcripts and was significantly associated with increased intensities of dermo infection. These genes showed a positive correlation between dermo infection and expression levels with the highest infection group showing the strongest association with the module (Figure 3B). Of the 24 genes with the highest module membership scores (KmE > 0.75), 2 have been identified in other studies as being involved in immunity in oysters (Wang et al., 2010;Li et al., 2017). In addition, module pink was found to also contain isoform 1 of the ninjurin-1-like gene that was also identified as differentially expressed in the pairwise comparison. Functional enrichment of this module identified enrichment of supramolecular fiber organization, intracellular, and cellular carbohydrate metabolic process. Comparing the list of genes in module Pink and genes identified as differentially expressed in Proestou and Sullivan (2020) identified 16 genes in both datasets. These genes included two isoforms of a multimerin-2-like gene, a calmodulin-2/4-like gene, a HSP70 12A-like gene, and a prolinerich transmembrane protein. This module was also found to have significant GO enrichment for 8 MF terms, 13 BP terms, and 8 CC terms (Figure 3C).

Changes in DNA Methylation and Gene Expression
We examined the relationship between gene expression and methylation for gene promoters (2 kb up-stream of first exon), gene bodies, downstream regions (2 kb down stream of last exon), and all intragenic regions combined (e.g., region covered by promoter, gene body, and downstream region). After breaking the expression into deciles, we used a Kruskal-Wallis test to compare mean methylation of each decile. This analysis identified a near-significant association between increasing expression and increased promoter methylation (p-value = 0.065), a significant positive interaction association between expression and gene body methylation (p-value < 2.2e −16 ; Figure 4A), and no interaction association between expression and downstream methylation (p-value = 0.11). We also explored the relationship between the mean methylation and the coefficient of gene expression variation for each gene calculated across all samples. All but one ontology was found to be enriched among downregulated transcripts in the moderate-heavily infected individuals. Percent of reference shows the number of genes identified as enriched when compared to the total number of genes for that GO term in the background reference list.
This analysis identified a significant relationship between decreasing methylation and increasing variation in gene expression (p-value < 2.2e −16 ; Figure 4B). However, there was only significant overlap between 2 DMRs and 2 DEGs and a Fisher's Exact test revealed no significant overlap between differential expression and differential methylation (pvalue > 0.05).

DISCUSSION
Our study described how P. marinus infection intensities influence DNA methylation and gene expression in the eastern oyster C. virginica. Through a combination of reduced representation DNA methylation sequencing and 3 -RNA sequencing (TAGseq) we identified significant changes in DNA methylation and subtle changes in gene expression that support previous findings of transcriptomic responses to infection and confirms the relationship between percent DNA methylation and the magnitude and variation in gene expression. The ability to assess these changes across multiple individuals (n = 40) strengthened these observations and encourages the use of gene body percent methylation as a proxy for expression and reinforces a potential role for plasticity. However, this study did not find compelling evidence for a relationship between differential methylation and differential expression. This lack of overlap suggests that the majority of changes in DNA methylation in response to infection might occur at any level of infection (i.e., all individuals are exhibiting abroad infection methylome); or, that these changes in DNA methylation are only necessary for acute responses that gradually return to either a seasonally or environmentally responsive methylation pattern. The lack of evidence for this direct connection may derive from the limited knowledge of when infection was initiated, and mechanisms governing changes to gene body methylation. These data suggest much of the observed variation in methylation does not reliably predict the direction of differential gene expression. This observation however, may also be a reflection of the TAGseq approach that is unable to identify changes in splice variant expression. It is therefore quite possible that the changes in DNA methylation are not associated with changes in expression at the gene level, but these changes could be influencing expression at the isoform level.

Changes in DNA Methylation in Response to P. marinus Infection
When we then look at the ontologies enriched in the genes that were differentially methylated between the moderate-heavy vs very-light infection categories, we see the regions that were hypomethylated in more heavily infected individuals were enriched for immune response genes, oxidative stress genes, and enzyme regulator activity. The hypomethylation of these regions may increase the plasticity of these genes given that lower methylation was also associated with greater variability in expression in our genome-wide analysis. That this occurs as a function of infection provides some evidence that the plasticity of these genes is being promoted. In contrast, regions hypermethylated in infected individuals were enriched for broad organizational categories associated with cellular structure and replication. If hypermethylation leads to canalization of expression, then it is possible that the increase in percent methylation of the gene bodies reinforces the stability of expression among these genes that may otherwise be altered by the activation of immune responses. The hypomethylation in heavily infected individuals of genes involved in apoptosis (caspase-8), pathogen clearance (C1q tumor necrosis factorrelated protein 4), and molecular chaperones (HSP 70 12Alike), suggests that moderate-heavy infections may be promoting higher plasticity of immune response genes (de Lorgeril et al., 2011;Wang et al., 2018). These potential roles are still largely based on the hypothesis that large changes in DNA methylation will drive phenotypic responses. As this still remains untested in bivalves, the community would greatly benefit from controlled studies investigating changes in DNA methylation across time in the same individuals.

Changes in Gene Expression in Response to P. marinus Infection Intensity
Differential gene expression analysis identified only 39 genes that showed differential expression between the moderateheavy (>100,000 parasites g −1 ) and very-light infected (≤1,000 parasites g −1 ) individuals (Figure 2). This low level of differential gene expression is comparable to the dampened response (21 genes DEG) observed in a dermo resistant family when measured 28-days post-exposure in a controlled experiment (Proestou and Sullivan, 2020). The genes that were differentially expressed in both studies (day 28; Proestou and Sullivan, 2020) included 2 forms of the complement C1q tumor necrosis factor (protein 6-like and protein 7-like) and a nuclear polyadenylated RNA-binding protein 3-like. The fact that these two datasets provide similar findings of a small number of genes being differentially expressed suggests two major takeaways; (1) that the individuals observed in this study were likely suffering from longterm infections of P. marinus; and, (2) that the shared differential expression of genes associated with tumor necrosis factors (complement system C1q) means these genes in particular are good candidates for further investigation as potential markers of disease resistance. Assuming that the individuals were suffering from long-term infections is consistent with other observed infection rates at this out-plant site (La Peyre et al., 2018). For the genes identified by both our study and the previous controlled exposure (complement C1q tumor necrosis factor-related protein 7-like and multimerin-2-like) we see an approximate sixfold increase in expression between the moderate-heavy and verylight infection categories. The role of these genes in immune responses is further supported by recent assessments of the functional role of these gene families in bivalve innate immunity (Gerdol et al., 2019). We also observed two versions of a ninjurin 1-like gene that were both up-regulated in the moderate-heavy infected individuals. Ninjurin-like genes have been found to be upregulated in Pacific Blue Shrimp (Litopenaeus stylirostris) that survive infection with Vibrio penaeicidia (de Lorgeril et al., 2005). We also found strong upregulation of two forms of a phytanoyl-CoA dioxygenase domain-containing protein 1-like that is putatively involved in peroxisomal lipid metabolism, a potential signal of a prolonged immune response in Manila clams (Venerupis philippinarum) responding to Perkinsus olseni (Romero et al., 2015). Together, these differential expression data highlight potential gene targets for future selective breeding programs aimed at increasing disease resistance.

WGCNA Module Pink Shows Responses to P. marinus Infection Intensity
WGCNA identified 1 module that was significantly correlated with dermo infection (module pink). Comparing the genes in module pink with the results from Proestou and Sullivan (2020) found 18 genes that were identified as DE after 28 days of infection in the dermo susceptible group and 3 genes identified as DE after 28 days in the dermo resistant group in their study. Of these genes the multimerin-2-like gene appears to be associated with angiogenesis and tumor formation (Khan et al., 2017). The positive association of infection of this gene within module Pink may suggest a role of angiogenesis in healing tissues degraded by P. marinus secreted proteases (La Peyre et al., 1995). In addition, multimerin-2-like genes have been proposed to have similar domains to complement component C1q genes (Gerdol et al., 2019) that were significantly differentially expressed in the above pairwise analysis of gene expression. In addition, two genes within this module, Thymosin beta-4 and kelch-like protein, have been associated with immune defense against pathogens, Thymosin beta-4 is also involved in antibacterial activity in C. gigas (Nam et al., 2015) and was found to facilitate the clearance of Vibrio alginolyticus in C. hongkongensi . Kelch-like proteins have also been shown to be differentially regulated in C. virginica when challenged with dermo (Wang et al., 2010). Our results show that both of these genes are up-regulated as P. marinus infection intensity increases. These findings suggest a potential role for both genes in the immune response of C. virginica infected with dermo. However, further studies are needed to clarify the role these genes play in defending against dermo infection in oysters. This overlap between groups further supports the module pink as a gene network involved in immune responses and provides additional gene targets for better understanding the variability in susceptibility to P. marinus across oyster families. These results also further support the use of weighted gene co-expression network analysis for exploring genomic patterns of expression for disease resistance in marine invertebrates.

Differential Methylation and Differential Expression
Global levels of DNA methylation across gene features revealed some unexpected associations, with limited support for a role of methylation among promoter and downstream regions on gene expression levels. Previous research in the Pacific oyster has found evidence for a putative role of promoter methylation in regulating gene expression during development (Rivière et al., 2013;Rivière, 2014). In our study, we tested this directly and found a nearly significant association (p-value = 0.065) between gene expression and mean promoter methylation. This lack of significant association is surprising, and it is possible that these differences in observations maybe associated with sampling during development (Rivière et al., 2013) or that only some genes are influenced by promoter methylation leading to a nearly significant result. Other more recent studies have reported that methylation increases at the start of the gene body and not upstream in promoter regions . Rather, gene body methylation showed a strong positive correlation with gene expression (Figure 4), a finding that has been previously reported in the Pacific oyster (Gavery and Roberts, 2013;Olson and Roberts, 2014). In our study we found that genes with higher mean methylation have higher overall expression and a lower coefficient of variation of expression. When comparing mean methylation with gene expression, we found that there were four categories ( Figure 4A) of mean methylation associated with increasing gene expression, but, eight categories ( Figure 4B) of mean methylation associated with increasing variation in expression. This suggests that mean methylation is a better predictor of gene variation than it is of the magnitude of gene expression. These results align well with previous research in Pacific oyster where it has been proposed that higher levels of gene body methylation are associated with constitutively expressed genes while lower levels of gene body methylation are associated with "environmentally responsive" gene categories (Roberts and Gavery, 2012;Song et al., 2017). As such, it is plausible that genes showing a decrease in gene body methylation (hypomethylated) are modified in a manner that may lead to an increase in plasticity of that gene. This potential increase in plasticity will require additional studies, ideally consisting of time-course observations of a small number of individuals in order to determine if an environmentally induced change in methylation will lead to future changes in the plasticity of that gene. In addition, future studies should further investigate the role of differential methylation on alternative splicing events. Regardless, these data suggest a relatively minor role of differential DNA methylation on observed differences in gene expression among individuals experiencing chronic exposure to infections.

CONCLUSION
This study sought to understand how chronic exposure to sub-lethal infections with the protistan parasite P. marinus influenced DNA methylation and gene expression in eastern oysters. Our results have confirmed a relationship between DNA methylation and gene expression such that lower levels of DNA methylation are found among genes with higher variations in expression, while highly methylated genes are found to have higher but less variable constitutive levels of expression. Genes with low gene body methylation included those involved in immune system process and enzyme regulator activity. Gene expression analyses found a limited response to infection, similar to the differential expression levels seen among dermo resistant families after 28 days of exposure (Proestou and Sullivan, 2020). The concordant differential expression results between ours and this previous study for C1q tumor necrosis factor genes suggests additional studies should explore the function of these genes in infection resistance. Finally, the significant but limited overlap in genes showing differential methylation and differential expression is a potential indicator of a limited influence of changes in methylation on changes in expression.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the NCBI SRA database under BioProject PRJNA604121.

AUTHOR CONTRIBUTIONS
MK and JL conceived of and designed the outplant study. KJ, MK, JL, SC, and KS sampled the animals for this study. SC and JL measured infection intensities. KJ prepared RRBS libraries and performed differential methylation analysis. KS and KJ analyzed the transcriptomic data and wrote the manuscript with input from SC, JL, and MK. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank the Louisiana Department of Fisheries and Wildlife, Dr. John Supan, and Dr. Brian Callam for collecting and spawning the oysters used in this study. We would also like to thank Bennett Thomas and Michelle Gautreaux for their assistance in RNA extractions, Joanna Griffiths for feedback on data analysis. Finally, we thank Alicia Reigel and Megan Guidry for assistance in sampling the oysters used in this study.