# Methods

## Samples and experimental design

DNA methylation was profiled in *Crassostrea gigas* (Pacific oyster) under a
two-factor factorial design crossing **ploidy** (diploid, 2N; triploid, 3N)
with **seawater pH** (ambient/high; acidified/low), with six biological
replicates per treatment combination (n = 24 total; 6 × [2N-high, 2N-low,
3N-high, 3N-low]). Samples were labelled 2H-1…2H-6 (2N/high), 2L-1…2L-6
(2N/low), 3H-1…3H-6 (3N/high) and 3L-1…3L-6 (3N/low).

## Library preparation and alignment

Whole-genome bisulfite sequencing libraries were prepared and sequenced, and
reads were aligned upstream of this analysis with Bismark. Per-sample,
strand-collapsed CpG methylation calls were exported as Bismark cytosine-report
coverage files (`*.CpG_report.merged_CpG_evidence.cov`), each a six-column
BED-like table (chromosome, start, end, methylation percentage, count
methylated, count unmethylated). These 24 coverage files were the input to the
present analysis. Positions are reported on the *C. gigas* reference assembly
cgigas_uk_roslin_v1 (NCBI RefSeq accession **GCF_902806645.1**); coverage files
and the annotation GFF share the same RefSeq sequence identifiers (e.g.
NC_047559.1; mitochondrion NC_001276.1).

## Software environment

All analyses were performed in R (v4.5.3) using Bioconductor (v3.22) packages
**methylKit** (v1.36.0), **DSS** (v2.58.0), **bsseq** (v1.46.0),
**GenomicRanges** (v1.62.1), **genomation** (v1.42.0) and **rtracklayer**
(v1.70.1), with data handling in **data.table** and figures in **ggplot2** /
**patchwork** / **pheatmap**. Random-number-dependent steps used a fixed seed
(`set.seed(1)`). Parallelised steps used the number of cores noted below.

## Reading, quality control and filtering

Coverage files were imported with methylKit `methRead()`
(`pipeline = "bismarkCoverage"`, `header = FALSE`, `context = "CpG"`,
`mincov = 1`, `assembly = "cgigas_roslin_v1"`), assigning each sample a
treatment code (2N-high = 0, 2N-low = 1, 3N-high = 2, 3N-low = 3). Per-sample
quality control summarised CpG counts, coverage (mean, median, fraction of CpGs
at ≥5× and ≥10×) and the methylation-level distribution. Each library covered
~9.5–10.9 million CpGs; mean CpG coverage ranged 6.5–12.5×, with 54–74% of CpGs
at ≥5× in every sample. The methylation landscape was bimodal (mosaic), typical
of invertebrate genomes (≈83% of CpGs ≤10% methylation, ≈5% ≥90%; mean ≈10–12%).
No sample was excluded.

CpGs were filtered with `filterByCoverage()` to remove positions below **5×
coverage** (`lo.count = 5`) and above the **99.9th coverage percentile**
(`hi.perc = 99.9`) to guard against PCR/multimapping artefacts, then
coverage-normalised across samples with `normalizeCoverage(method = "median")`.
Between 54% and 74% of CpGs per sample were retained (5.2–8.0 million).

## Sample structure

The filtered objects were merged with `unite(destrand = FALSE)`, requiring a CpG
to be covered in all 24 samples, yielding **1,953,620 common CpGs**. Per-sample
percent-methylation matrices (`percMethylation()`) were used to compute pairwise
Pearson correlations, a Ward.D2 hierarchical clustering of 1 − correlation, and
a principal-component analysis (`prcomp`, centred and scaled) on the 200,000
most-variable CpGs.

## Global methylation

Weighted global methylation per sample was computed over the filtered CpGs as
100 × Σ(methylated counts) / Σ(coverage). The effect of the design on global
methylation was tested with a two-way ANOVA,
`aov(global_methylation ~ ploidy * pH)`, reporting the ploidy main effect, pH
main effect and ploidy × pH interaction.

## Differential methylation (per-CpG DMLs)

Pairwise differential methylation was assessed with methylKit
`calculateDiffMeth()` (logistic regression with overdispersion correction;
`mc.cores = 4`), computing methylation differences as treatment (group 1) minus
control (group 0). For each contrast the relevant samples were subset with
`reorganize()` and re-united (`destrand = FALSE`) before testing. Differentially
methylated loci (DMLs) were defined as **q < 0.01 and |methylation difference|
≥ 15%** (q-values from SLIM correction as implemented in methylKit). Three
contrasts were run: (i) pH within diploids (2N low vs high), (ii) pH within
triploids (3N low vs high) and (iii) the ploidy main effect (3N vs 2N, on the
united 24-sample object).

## Ploidy × pH interaction test

Because methylKit does not provide an interaction term, the central
ploidy × pH interaction was tested per CpG on the 1,953,620 common CpGs with a
binomial generalised linear model,
`glm(cbind(methylated, unmethylated) ~ ploidy * pH, family = binomial())`. For
each CpG the Wald-test p-value for the interaction coefficient
(`ploidy3N:pHlow`) was extracted; p-values were Benjamini–Hochberg adjusted
across all CpGs. Computation was chunked and parallelised with
`parallel::mclapply(mc.cores = 6)`. The interaction effect size was defined as
the difference of level-differences from group weighted means:
[(3N low − 3N high) − (2N low − 2N high)]. Interaction DMLs were called at
**q < 0.05 and |interaction effect| ≥ 15%**. Each significant CpG was classified
by the sign of the pH response in each ploidy as a *directional reversal*
(response of opposite sign between ploidies) or a *magnitude change*
(same sign, different amplitude).

## Differentially methylated regions (DMRs)

As a region-level, memory-efficient complement, the genome was tiled into
non-overlapping 1-kb windows with
`tileMethylCounts(win.size = 1000, step.size = 1000, cov.bases = 3)` (requiring
≥3 covered CpGs per window). The three contrasts above were repeated on the
tiled object with `calculateDiffMeth()`, and DMRs were called at the same
thresholds used for DMLs (q < 0.01, |methylation difference| ≥ 15%).

## Genomic-feature annotation

DMLs were annotated against the NCBI RefSeq GFF for GCF_902806645.1
(`*_genomic.gff.gz`), imported with `rtracklayer::import()`. Gene, exon, CDS,
mRNA and lncRNA features were extracted as GRanges; introns were derived as
`setdiff(reduce(genes), reduce(exons))` and promoters as the 1-kb region
upstream of each gene's transcription start. Each DML (and, as background, all
tested CpGs) was assigned to a single feature using a fixed priority hierarchy:
**CDS > exon > intron > promoter > intergenic**. Feature proportions were
compared between each DML set and the genome-wide CpG background.

## GO enrichment

DMLs were mapped to genes by genomic overlap (`GenomicRanges::findOverlaps`).
Gene-to-GO annotations were obtained by joining each gene's numeric NCBI GeneID
(parsed from the GFF `Dbxref`) to transcript identifiers and then to Gene
Ontology terms via the project's `geneid2go` table; GO term names were taken
from the accompanying GO-slim annotation. Term-level enrichment was tested with
the hypergeometric distribution,
`phyper(k − 1, N, N_universe − N, n_selected, lower.tail = FALSE)`, where the
universe was restricted to genes overlapping at least one tested CpG and terms
were required to contain ≥5 annotated genes; p-values were Benjamini–Hochberg
adjusted, and terms with q < 0.05 were considered enriched. Fold enrichment was
computed as observed/expected gene counts.

## Data and code availability

Sample metadata was reconstructed from the analysis code accompanying the
coverage files (the original sample sheet was not available with the public
data); it encodes ploidy and pH only, and additional covariates (e.g. tank,
family, sequencing batch), if recorded, were not modelled. All intermediate
tables (per-sample QC, global methylation, DML/DMR result tables, feature
distributions, GO enrichment) and figures are provided as supplementary
artifacts.