Mechanisms and disease associations of haplotype-dependent allele specific DNA methylation: Methyl-seq data for the identification of hap-ASM

Haplotype-dependent allele-specific methylation (hap-ASM) can impact disease susceptibility, but maps of this phenomenon using stringent criteria in disease-relevant tissues remain sparse. Here we apply array-based and Methyl-seq approaches to multiple human tissues and cell types, including brain, purified neurons and glia, T lymphocytes, and placenta, and identify 795 hap-ASM differentially methylated regions (DMRs) and 3,082 strong methylation quantitative trait loci (mQTLs), most not previously reported. More than half of these DMRs have cell type-restricted ASM, and among them are 188 hap-ASM DMRs and 933 mQTLs located near GWAS signals for immune and neurological disorders. Targeted bis-seq confirmed hap-ASM in 12/13 loci tested, including CCDC155, CD69, FRMD1, IRF1, KBTBD11, and S100A*-ILF2, associated with immune phenotypes, MYT1L, PTPRN2,CMTM8 and CELF2, associated with neurological disorders, NGFR and HLA-DRB6, associated with both immunological and brain disorders, and ZFP57, a trans-acting regulator of genomic imprinting. Polymorphic CTCF and transcription factor (TF) binding sites are over-represented among hap-ASM DMRs and mQTLs, and analysis of the human data, supplemented by cross-species comparisons to Macaca mulattamacaques, indicates that CTCF and TF binding likelihood predicts the strength and direction of the allelic methylation asymmetry. These results show that hap-ASM is highly tissue-specific; an important trans-acting regulator of genomic imprinting is regulated by this phenomenon; variation in CTCF and TF binding sites is an underlying mechanism in primary tissues, and maps of hap-ASM and mQTLs reveal regulatory sequences underlying supera- and sub-threshold GWAS peaks in immunological and neurological disorders. We used the Agilent SureSelect Methyl-seq DNA hybrid capture kit, followed by deep Nextgen bis-seq. In this protocol, targeted regions (total of 3.7M CpGs) include RefGenes, promoter regions, CpG islands, CpG island shores, shelves, and DNAse I hypersensitive sites. DNA was sheared to an average size of 250 bp and bisulfite converted with the EZ DNA methylation kit (Zymo). Paired-end reads (150 bp) were generated with an Illumina HiSeq2000 sequencer. One of the brain samples (5984) was relatively under-represented in the library, so additional sequences were generated on a MiSeq sequencer to improve the coverage for this sample. After trimming for low-quality bases (Phred score<30) and reads with a length <40 bp with TrimGalore, the reads were aligned to the human genome (GRCh37) using Bismark and duplicate reads were removed using Samtools. Methylation calling was performed using Bismark and and SNP calling was performed with BisSNP with default settings and dbSNP137 as references. ASM calling was performed with Bismark, after separating the valid SNP-containing reads by allele.