## Interaction between methylation and expression across genes To understand how methylation may influence expression of distal genes, we conducted a sparse partial least squares discriminant analysis (sPLS-DA) using the mixOmics package in R [(Rohart et al. 2017)](https://paperpile.com/c/3sGCWx/yJJV). Prior to analysis, genes with methylation and expression variance less than five were removed, as features with near-zero variance could skew the analysis. Average methylation and expression data for each gene was filtered. Genes with methylation and expression variance less than five were removed.The sPLS-DA was run keeping the 25 top genes from both methylation (X) and expression (Y) datasets that explained any treatment differences for either the first or second component of the sPLS-DA (keepX = c(25,25), keepY = c(25,25). The contribution of the top 25 X and Y variables in driving differences by treatment was assessed (selectVar) for each component. Pearson’s correlations between the X and Y variables on each component were assessed using a cim plot (cim). Relationships between X and Y variables were further analyzed using network plots for Pearson’s correlations > 0.8 (network). A sPLS-DA was conducted for female and male oysters separately to determine how methylation may influence expression of distal genes. The top 25 genes that explained pH treatment differences for each component of the sPLS-DA ([Supplementary Table 5](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Fem.txt)). For female oysters, methylation drivers with the highest loadings for the first component included MAP kinase-interacting serine/threonine-protein kinase 1-like (-0.34) and zinc finger protein 236-like (-0.30) and THAP domain-containing protein 10-like (0.60) for the second component ([Supplementary Table 5](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Fem.txt)). RNA-binding protein Nova-1-like (0.44) and sorting nexin-33-like (0.43) had the strongest loadings for expression drivers on the first component, and adenylate kinase 7-like (0.39) and mitochondria-eating protein-like (0.35) had the strongest loadings on the second component ([Supplementary Table 5](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Fem.txt)). A clustered image map of Pearson’s correlations between driver genes for methylation and expression datasets revealed one group of genes that was predominantly positively correlated, and another that was mostly negatively correlated (Figure 5A). A similar trend was found in a network plot of all driver genes (Figure 5B). When considering Pearson’s correlations > 0.8, two independent networks were identified, with one consisting of driver genes positively correlated with eachother, and another with negatively correlated driver genes. Genes with signal transduction, cell organization and biogenesis, cytoskeletal activity, and enzyme regulator activity functions were found in the negatively correlated cluster associated with the first component of the sPLS-DA ([Supplementary Table 5](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Fem.txt)). The positively correlated cluster was associated with the second sPLS-DA component, had genes involved in RNA metabolism, developmental processes, nucleic acid binding, and transcription regulatory activity ([Supplementary Table 5](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Fem.txt)). For male oysters, various uncharacterized proteins had the highest loadings for methylation drivers in the first component, with hemicentin-2-like (0.43) having the highest loading for methylation drivers in the second component ([Supplementary Table 6](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Male.txt)). Abhydrolase domain-containing protein 2-like (0.40) and AN1-type zinc finger protein 4-like (0.32), and E3 ubiquitin-protein ligase TRIM71-like (0.35) had the highest loading for expression drivers in the first and second components, respectively ([Supplementary Table 6](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Male.txt)). Two clusters of driver genes were identified in clustered image maps and network plots, and each cluster was associated with a sPLS-DA component. The first component was associated with a network consisting of both positive and negative correlations (Figure 5C). This network had genes involved in cell organization and biogenesis, developmental processes, and cytoskeletal activity ([Supplementary Table 6](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Male.txt)). Another network with predominantly positive correlations was associated with the second sPLS-DA component, and had cell organization and biogenesis, developmental processes, and signal transduction activity genes (Figure 5D; [Supplementary Table 6](https://github.com/sr320/ceabigr/blob/main/output/gene-methylation-expression-mixomics/selectedVariables-Annot-Male.txt)).