Polyploidy modulates the DNA-methylation response to ocean acidification in Crassostrea gigas

Whole-genome bisulfite (BS-seq) analysis of diploid vs triploid Pacific oysters under high and low pH

2×2 factorial design · ploidy (2N/3N) × pH (high/low) · n = 6 per cell · 24 samples · genome: cgigas_uk_roslin_v1 (GCF_902806645.1) · analysis: methylKit + per-CpG GLM interaction test

1. Study design & samples

Bisulfite alignment (Bismark) was completed upstream; this analysis began from per-sample CpG-report coverage files (.cov). Sample metadata was reconstructed from the provided code snippet (the original sample_metadata.csv is not on the public server). The design is balanced: 6 replicates in each of four groups.

2. Coverage & methylation QC

Each library covers ~9.5–10.9 M CpG positions. Mean CpG coverage ranged 6.5–12.5×; 54–74% of CpGs reached ≥5× in every sample. No sample was an outlier for exclusion. The methylation landscape is the classic invertebrate mosaic / bimodal pattern — most CpGs are near-unmethylated (~83% ≤10%) with a smaller highly-methylated mode (~5% ≥90%), and mean methylation is low (~10–12%) and remarkably uniform across all 24 samples.

3. Sample structure

Genome-wide, samples do not cluster by treatment. Pairwise correlations are uniformly high (0.85–0.95) and PC1 explains only 13.2% of variance. The mild separation in the correlation heatmap tracks sequencing coverage, not ploidy or pH. This is expected and is the key motivation for the per-CpG differential testing that follows: treatment effects are real but localized to a small fraction of loci, so they do not dominate global structure.

4. Global methylation: the interaction emerges

Weighted global methylation per sample was tested with a 2×2 ANOVA. The ploidy×pH interaction is significant (p = 0.018) while neither main effect is (ploidy p = 0.13, pH p = 0.63). The pattern is a crossover: acidification raises methylation in diploids (10.43 → 10.80%) but lowers it in triploids (11.16 → 10.63%).

5. Differential methylation — main contrasts

Per-CpG differential methylation used methylKit's logistic regression (calculateDiffMeth), calling DMLs at q < 0.01 and |methylation difference| ≥ 15%. The two within-ploidy pH contrasts reveal the flip directly: diploids respond to acidification with predominantly hypermethylation, triploids with predominantly hypomethylation — and triploids respond at more loci.

6. Ploidy × pH interaction — the central test

methylKit has no interaction term, so a per-CpG binomial GLM cbind(methylated, unmethylated) ~ ploidy * pH was fitted across all 1,953,620 CpGs covered in every sample, extracting the Wald test on the ploidy3N:pHlow coefficient (BH-adjusted). At q < 0.05 and |interaction effect| ≥ 15%, 17,053 CpGs are interaction-significant. Classifying each by how the pH response differs between ploidies:

• 12,271 loci — acidification methylates in 2N but demethylates in 3N (reversal)
• 4,097 loci — the opposite reversal (demethylates in 2N, methylates in 3N)
• 361 + 316 loci — same direction but amplified in one ploidy

16,368 of 17,053 (96%) are directional reversals, not mere magnitude differences. Polyploidy fundamentally re-wires which loci gain and lose methylation under acidification.

7. Region-level DMRs (1 kb tiles)

As a memory-safe, region-level complement, the genome was tiled into 1 kb windows (≥3 CpG) and the same contrasts re-run. DMRs reproduce the CpG-level direction flip cleanly.

8. Genomic-feature & GO annotation

DMLs were annotated against the NCBI RefSeq GFF (GCF_902806645.1) with a hierarchical assignment (CDS > exon > intron > promoter > intergenic). Across all contrasts, DMLs are strongly depleted in intergenic space (6–7% vs 25% of background CpGs) and enriched in gene bodies (introns 50–56% vs 36%; CDS 28–35% vs 26%) — the canonical signature of functional invertebrate gene-body methylation.

DML-associated genes were mapped to GO terms via the study's geneid2go annotation and tested for enrichment (hypergeometric, background = all genes overlapping a tested CpG). Genes carrying the ploidy-dependent (interaction) signal are enriched for cytoskeleton and motor function (dynein complex, microtubule plus-end binding, actomyosin organization, centrosome localization), signalling (GTPase-activator activity, TGF-β receptor signalling, Ras GTPase binding), and transcriptional regulation. Top interaction terms:

9. Interpretation

The data give a consistent, multi-scale answer to the study question: polyploidy does not just dampen or amplify the epigenetic response to ocean acidification — it reverses its direction at thousands of loci. The signal is coherent from the genome-wide average (ANOVA interaction p = 0.018), through per-CpG DMLs (17k interaction sites, 96% reversals), to 1 kb regions, and it concentrates in gene bodies of cytoskeletal, signalling and transcriptional-regulation genes. Diploids and triploids appear to deploy fundamentally different methylation programs to cope with low pH, which may underlie ploidy differences in stress physiology and survival.

Methods summary

• Input: 24 Bismark CpG-report .cov files; ~10 M CpGs/sample.
• Load/QC: methylKit methRead (bismarkCoverage); filter lo.count=5, hi.perc=99.9; normalizeCoverage(median).
• Common set: unite(destrand=FALSE) → 1,953,620 CpGs in all 24 samples.
• DML: calculateDiffMeth (logistic regression); q < 0.01, |diff| ≥ 15%.
• Interaction: per-CpG binomial GLM cbind(Cs,Ts) ~ ploidy*pH, Wald test on interaction term, BH-adjusted; q < 0.05, |effect| ≥ 15%.
• DMR: tileMethylCounts (1 kb, step 1 kb, ≥3 CpG) + calculateDiffMeth.
• Annotation: NCBI RefSeq GFF GCF_902806645.1; GO via study geneid2go; hypergeometric enrichment, BH.