Identifying and Analyzing Differentially Methylated Regions.

DMR identification and annotation followed previously published aapproaches (65, 66). The FastQC program (67) was used to assess data quality. Reads were trimmed to remove adapters and low-quality bases using Trimmomatic (68). The reads for each MeDIP sample were mapped to the P. mexciana (48) genome using Bowtie2 (69) with default parameter options. The mapped read files were then converted to sorted BAM (Binary Sequence Alignment/Map) files using SAMtools (70).

Differential coverage between sulfidic and nonsulfidic populations was calculated using the MEDIPS R package (71). P value from edgeR (72) was used to determine the significance of the difference between the two groups for each 100-bp genomic window. Windows with an edgeR P value less than a specified threshold (P < 10⁻⁷) were considered the initial start of the DMR. DMR edges were extended until no genomic window with a P value less than 0.1 remained within 1,000 bp of the DMR. The extended DMR overlap compared the DMRs, with at least one 100-bp window with P < 10⁻⁷ from one comparison, with the genomic windows (100-bp regions) in a second comparison. Windows that had a P value <0.05 in the second comparison were considered overlapping. The Ensembl database (73), accessed with the biomaRt R package (74), was used to annotate DMRs. Genes that were overlapping a DMR, including 10 kb on either side of the DMR, were input into a Kyoto Encyclopedia of Genes and Genomes pathway search (75, 76) to identify relevant associated pathways. The DMR-associated genes were sorted into functional groups using information provided by the DAVID (Database for Annotation, Visualization and Integrated Discovery) (77) and Panther (78) databases incorporated into an internal curated database (https://www.skinner.wsu.edu/). DMR-associated gene correlations present in published literature were further analyzed using Pathway Studio software (version 12.2.1.2: Database of functional relationships and pathways of mammalian proteins; Elsevier). R code computational tools are available at GitHub (https://github.com/skinnerlab/MeDIP-seq) and https://www.skinner.wsu.edu.