
==== Front
Cell Genom
Cell Genom
Cell Genomics
2666-979X
Elsevier

S2666-979X(24)00199-X
10.1016/j.xgen.2024.100605
100605
Article
Crosstalk between epitranscriptomic and epigenomic modifications and its implication in human diseases
Li Chengyu 127
Chen Kexuan 127
Fang Qianchen 12
Shi Shaohui 12
Nan Jiuhong 12
He Jialin 12
Yin Yafei 3
Li Xiaoyu 3
Li Jingyun 4
Hou Lei 5
Hu Xinyang 23
Kellis Manolis 6
Han Xikun han001@mit.edu
6∗
Xiong Xushen xiongxs@zju.edu.cn
128∗∗
1 The Second Affiliated Hospital & Liangzhu Laboratory, Zhejiang University School of Medicine, Hangzhou 311121, China
2 State Key Laboratory of Transvascular Implantation Devices, The Second Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou 311121, China
3 The Second Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, Zhejiang 310058, China
4 Sir Run Run Shaw Hospital, Zhejiang University School of Medicine, Hangzhou 310016, China
5 Department of Medicine, Biomedical Genetics Section, Boston University, Boston, MA 02118, USA
6 Computer Science and Artificial Intelligence Lab, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
∗ Corresponding author han001@mit.edu
∗∗ Corresponding author xiongxs@zju.edu.cn
7 These authors contributed equally

8 Lead contact

08 7 2024
14 8 2024
08 7 2024
4 8 10060510 1 2024
17 4 2024
14 6 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
Summary

Crosstalk between N6-methyladenosine (m6A) and epigenomes is crucial for gene regulation, but its regulatory directionality and disease significance remain unclear. Here, we utilize quantitative trait loci (QTLs) as genetic instruments to delineate directional maps of crosstalk between m6A and two epigenomic traits, DNA methylation (DNAme) and H3K27ac. We identify 47 m6A-to-H3K27ac and 4,733 m6A-to-DNAme and, in the reverse direction, 106 H3K27ac-to-m6A and 61,775 DNAme-to-m6A regulatory loci, with differential genomic location preference observed for different regulatory directions. Integrating these maps with complex diseases, we prioritize 20 genome-wide association study (GWAS) loci for neuroticism, depression, and narcolepsy in brain; 1,767 variants for asthma and expiratory flow traits in lung; and 249 for coronary artery disease, blood pressure, and pulse rate in muscle. This study establishes disease regulatory paths, such as rs3768410-DNAme-m6A-asthma and rs56104944-m6A-DNAme-hypertension, uncovering locus-specific crosstalk between m6A and epigenomic layers and offering insights into regulatory circuits underlying human diseases.

Graphical abstract

Highlights

• Identification of directional crosstalk between m6A and epigenomic modifications

• Crosstalk between m6A and DNA methylation shows a global cross-tissue consistency

• Differential genomic localization between forward and reverse regulatory directions

• 2,036 GWAS loci interpreted by the crosstalk among m6A, DNA methylation, and H3K27ac

Li et al. utilized the quantitative trait loci of m6A, DNA methylation, and H3K27ac as genetic instruments to delineate locus-specific, directional maps of the crosstalk between m6A and the two fundamental epigenomic traits. The genetic crosstalk maps were further integrated with GWASs to uncover and interpret disease variants that are dependent on such m6A-epigenome crosstalk mechanisms.

Keywords

m6A crosstalk
Mendelian randomization
GTEx
GWAS
colocalization
H3K27ac
DNA methylation
Published: July 8, 2024
==== Body
pmcIntroduction

Genome-wide association studies (GWASs) have cataloged over 500,000 associations between genetic variants and over 8,000 human traits.1,2,3 However, more than 93% of the disease-relevant variants are located in non-coding regions, with no direct effect on protein sequences.4,5,6 The molecular mechanisms underlying their contributions to human diseases remain largely unclarified, impeding the precise diagnosis of and therapeutics for personalized medicine.7,8 To address this challenge, studies have evaluated the quantitative genetic traits of gene expression, namely expression quantitative trait loci (eQTLs), in human tissues and cells to pinpoint the putative causal variants and their target genes, thereby providing mechanistic insights for the GWAS variants at the gene expression level.9,10

However, cellular regulatory processes involve highly complex layers beyond gene expression. The mechanisms by which GWAS risk variants exert their effects encompass multiple layers of molecular regulations. In addition to eQTLs, various types of QTLs for other molecular biomarkers, including DNA methylation (DNAme; mQTLs), multiple types of histone modification (haQTLs), chromatin accessibility (caQTLs), and N6-methyladenosine (m6A) RNA modification (m6A-QTLs), have been utilized to delineate the relationships across molecular levels, gene expression, and complex diseases.10,11,12,13,14,15,16,17,18,19,20,21,22 Understanding the cross-layer regulatory pathways across different molecular biomarkers is essential for bridging the mechanistic gap from genetics to diseases.

In fact, an increasing number of studies have revealed crosstalk between different epigenetic molecular levels. For instance, several recent studies have reported cross-layer interactions between m6A RNA modification (epitranscriptomic level) and various epigenomic layers.23,24,25,26,27,28,29,30,31,32,33 m6A has been established as the most prevalent and functionally essential chemical modification in message RNA (mRNA), leading to extensive attention in the field of epitranscriptomics.34,35,36,37,38,39 This epitranscriptomic mark has been well established for its roles in regulating molecular fates of mRNA, spanning mRNA processing, degradation, nuclear export, translation, and beyond.35,37,38,40,41,42 Accordingly, m6A plays crucial roles in various biological processes, such as genetic diseases, cancer, development, and stem cell differentiation.43,44,45 Notably, the installation process of m6A modification occurs concomitantly with transcription, a process referred to as co-transcriptional methylation, which provides a molecular basis for the crosstalk between m6A and epigenomes.46,47

Recent studies discovered that H3K36me3 histone modification and DNAme modulate the m6A methylation process, suggesting an “epigenome-to-m6A” regulation.32,33 Conversely, m6A has also been demonstrated to regulate epigenomic regulators, including chromatin accessibility, DNAme, and various types of histone modifications, representing an “m6A-to-epigenome” regulation.27,28,29,30,31 However, these studies each proposed a global directionality for the crosstalk without considering the actual regulatory mode at specific loci for different tissues/cell types, hindering further functional interrogation. Moreover, these studies have primarily focused on the regulatory roles of such crosstalk in stem cell differentiation and cancer progression, but its functional importance in human complex diseases has been untapped.27,28,29,30,31,32,33

Mendelian randomization (MR) is a statistical approach that has been widely used for inferring the causal relationships between genetically associated molecular traits and human diseases.48,49 MR relies on quantified genetic associations, including molecular QTLs and GWAS summary statistics.48,49 Recent studies have mapped the genetic associations for m6A, H3K27ac, and DNAme, generating m6A-QTLs, haQTLs, and mQTLs across various human tissues and cell lines.17,21,22,50,51,52 The availability of these QTL resources presents a unique opportunity to interrogate the causal directionality of the crosstalk between m6A and epigenomes in both tissue-specific and locus-specific contexts. In addition, integrating the genetically based crosstalk maps with GWASs bridges the mechanistic gap between genetic variants and human complex diseases.

Here, we employ genetic approaches to investigate the crosstalk between m6A and epigenomic layers, including H3K27ac and DNAme. We first design a bidirectional MR framework, which is further combined with colocalization analysis, to delineate the regulatory directionality between m6A and epigenomes. We identify, in total, 153 interaction loci between m6A and H3K27ac and 66,508 loci between m6A and DNAme. The crosstalk loci between m6A and the two epigenomic modifications show global directionality consistency across tissues, suggesting the robustness of these cross-layer interactions. We observe that the DNAme sites being regulated by m6A are enriched in enhancers and transcription start sites (TSSs) and depleted in repressed chromatin regions. Mechanistically, we predict RNA-binding protein (RBP) and transcription factor (TF) pairs that potentially mediate the crosstalk between m6A and epigenomes. Biologically, we identify 20 disease loci across various brain disorders that are dependent on the crosstalk between m6A and H3K27ac in brain and 1,767 and 249 disease loci in lung and muscle traits that are dependent on the crosstalk between m6A and DNAme in lung and muscle, respectively.

Results

Crosstalk directionality between m6A and epigenome elucidated by a bidirectional MR framework

To study the regulatory direction for the crosstalk between m6A and epigenome, we curated m6A-QTLs for m6A, mQTLs for DNAme, and haQTLs for H3K27ac assayed in brain, lung, muscle, and heart21,50,51 (Table S1). We designed a bidirectional MR framework, which utilizes mQTLs and haQTLs as instrumental variables to identify the epigenome-to-m6A regulation, as well as m6A-QTLs as instrumental variables to elucidate the m6A-to-epigenome regulation (Figures 1A and 1B; STAR Methods). The heterogeneity tests based on heterogeneity in dependent instruments (HEIDI) were performed to mitigate potential biases in the causal inference arising from pleiotropic single-nucleotide polymorphisms (SNPs).53 We carried out the colocalization analysis to further increase the robustness of the identified crosstalk54 (Figure 1B; STAR Methods). Considering the tissue specificity of these QTL types,21,22,50,51 we evaluated the crosstalk directionality between m6A and epigenomes for each tissue type individually.Figure 1 Study design

(A) Illustration of crosstalk directionality between m6A and epigenomes in a bidirectional MR framework.

(B) Datasets and statistical methods of the current study.

See also Table S1 and STAR Methods.

Utilizing this framework, we identified 58,379 DNAme-to-m6A regulation loci in lung and 3,396 in muscle, with 19,256 (33.0%) and 1,097 (32.3%) of these loci additionally validated in the colocalization analysis, respectively (Tables S2, S3, and S4; Figures S1A–S1C). We show an example where a DNAme locus (cg10788408) may regulate the m6A level of the TNFSF13 gene, which encodes a tumor necrosis factor ligand cytokine involved in the positive regulation of forced expiratory volumn in the first second (FEV1)55 (Figure S2A). In the reverse direction, we found 4,053 m6A-to-DNAme pairs in lung and 680 pairs in muscle, with 1,173 (28.9%) and 120 (17.6%) of these loci further supported by colocalization, respectively (Tables S2, S3, and S4; Figures S1A–S1C). For instance, cg07535628 methylation is predicted to be regulated by an m6A modification in the CD151 gene, which is highly expressed in lung tissue and exerts pleiotropic effects in multiple lung traits, including lung cancer, asthma, influenza, and idiopathic pulmonary fibrosis56 (Figure S2B). Notably, 2,396 (3.6%) crosstalk pairs were significant in both regulatory directions, suggesting potential bidirectional regulations at these genomic loci (Table S2; Figure S1A). Intriguingly, we noted one such bidirectional crosstalk at the NLRP1 locus, which is a NOD-like receptor that specifically expresses in lung epithelia and functions as a sensor of pathogenic coronavirus, including severe acute respiratory syndrome coronavirus 257 (Figure S2C).

For the crosstalk between m6A and H3K27ac, we identified, in total, 106 H3K27ac-to-m6A and 47 m6A-to-H3K27ac regulation loci across the four tissues, with 43 (40.6%) and 25 (53.2%) being further supported by colocalization, respectively (Tables S2 and S5; Figures 2A, S3A, and S3B). The relatively fewer significant MR events identified for the crosstalk between m6A and H3K27ac are likely due to the lower statistical power of the haQTLs compared to mQTLs (Table S1). As an example, a H3K27ac locus appears to be negatively regulated by a nearby m6A in the PNMA8B gene, which encodes a paraneoplastic Ma antigen family protein that plays multi-functional roles in paraneoplastic neurological diseases and cancer58,59 (Figure 2B). For the reverse direction in brain tissue, we found a signal indicating H3K27ac-to-m⁶A regulation for the m6A level of the gene STXBP1, which is an essential gene involved in multiple brain developmental disorders, including epilepsy, autism, and intellectual disability60,61,62 (Figure 2C).Figure 2 Crosstalk between m6A and epigenomes

(A) Manhattan plot shows the genomic position (x axis) and the −log10p value (y axis) of all the tested m6A-H3K27ac bidirectional MR pairs in the brain. MR results of m6A-to-H3K27ac are displayed on the top, and results of H3K27ac-to-m6A are shown on the bottom. Significant MR findings are highlighted in scarlet red and cyan with solid or hollow dots, where scarlet red means a positive effect size, cyan means a negative effect size, solid circles represent results further supported by colocalization analysis, and hollow circles represent results with only MR evidence. The m6A-modified genes of these crosstalk are labeled in italics. Sizes of the highlighted circles correspond to the values of MR effect sizes.

(B) An example of m6A-to-H3K27ac crosstalk near the PNMA8B gene in the brain. The box shows the 25th–75th percentiles (i.e., the interquartile range [IQR]), the line shows the median, and the whiskers show 1.5× IQR.

(C) An example of H3K27ac-to-m6A crosstalk near the STXBP1 gene in the brain.

(D) Cross-tissue effect size comparison of m6A-to-DNAme MR pairs between lung and muscle for the significant MR pairs identified in the lung. Pearson correlation coefficients were calculated using standard errors of MR effect sizes in the muscle as weights.

(E) Cross-tissue effect size comparison of m6A-to-DNAme MR pairs between lung and muscle for the significant MR pairs identified in the muscle. Pearson correlation coefficients were calculated using standard errors of MR effect sizes in the lung as weights.

See also Figures S1–S5 and Tables S2 and S5.

We additionally performed a summary data-based MR test, which is capable of distinguishing a pleiotropic model from a linkage model,53 and compared the results with the crosstalk maps obtained from the two-sample MR method used above (Figures S4A and S4B). The effect sizes between the two MR approaches are highly consistent (p < 0.001 across all tissues, R ranges from 0.969 to 1), suggesting the robustness of different MR methods and the reliability of the crosstalk we identified.

Tissue specificity of the crosstalk between m6A and epigenome

We next investigated the tissue specificity of the crosstalk between m6A and epigenomes. For the m6A-to-DNAme regulation, we found seven loci that were significant in both lung and muscle tissues. While the number of overlaps is limited, the regulatory directionality of the crosstalk was consistent across all seven loci between the two tissues (Figures 2D and 2E). Overall, the m6A-to-DNAme regulation shows significant consistency between lung and muscle, for both of the regulatory pairs identified as significant MR loci in lung (p = 2.7 × 10−21, R = 0.62; Figure 2D) and muscle (p = 6.1 × 10−3, R = 0.55; Figure 2E), based on a weighted correlation analysis that accounts for the standard errors of effect sizes. As a comparison, when focusing on the tissue specificity of the DNAme-to-m6A regulation, we observed weaker consistency for the regulatory pairs between lung and muscle (R = 0.095 and R = 0.072 for the loci significant in lung and in muscle under a false discovery rate < 0.1, respectively; Figures S5A and S5B). Of note, we observed increased correlation coefficients and directionality consistency between the two tissues when using more stringent thresholds to define significant crosstalk based on the MR analysis, suggesting that the observed consistency is likely an underestimation caused by statistical power (Figures S5C–S5F). Given the relatively limited number of significant interaction pairs between m6A and H3K27ac, we merged the interaction loci across the four tissues for consistency analysis. We identified, in total, six m6A-to-H3K27ac regulation pairs and two H3K27ac-to-m6A pairs that are shared across tissues (Figures S5G and S5H). These significant tissue-sharing pairs are consistent in their regulatory directionalities between tissues. Overall, the interactions between m6A and H3K27ac show high consistencies between tissues in both regulatory directions (Figures S5G and S5H; R = 0.49 and p = 0.027 for H3K27ac-to-m6A; R = 0.91 and p = 0.005 for m6A-to-H3K27ac pairs). Collectively, the crosstalk between m6A and the two epigenomic modifications demonstrate consistent regulatory directions between tissues, indicating the robustness of such cross-layer interaction.

m6As preferentially regulate DNAmes in enhancer regions

To understand the regional specificity of the crosstalk between m6A and epigenomes, we examined the genomic localization of the DNAmes that are predicted to interact with m6A in both forward and reverse directions based on the annotation of 18 chromatin states63,64 (Figures 3A, 3B, S6A, and S6B).Figure 3 Genomic distribution preference of the crosstalk loci

(A) Genomic enrichments of the DNAme sites being regulated by m6A in the lung against the 18 chromatin states annotations. p values were calculated using a two-sided Fisher’s exact test. Non-significant p values were denoted at the top/bottom of the enrichment bars (the same in B).

(B) Genomic enrichments of the DNAme sites involved in DNAme-to-m6A regulation in the lung against the 18 chromatin states annotations.

See also Figures S6 and S7 and Tables S6 and S7.

We observed that the m6A-regulated DNAme sites were enriched in enhancers and regions near TSSs and depleted for repressive regions (Figures 3A and S6A). Specifically, for lung, we found 1.82- and 1.65-fold of enrichment for genic enhancer2 (EnhG2, p = 3.5 × 10−6) and active enhancer1 (EnhA1, p = 3.4 × 10−15), respectively. This observation therefore supports the role of m6A in regulating gene expression via a mechanism of modulating enhancer activity. For muscle, we again observed significant enrichments of m6A-regulated DNAmes in enhancer regions, including EnhG1 and EnhA1 (Figure S6A). In addition, we noticed moderate enrichments of m6A-regulated DNAmes in the regions flanking TSSs for lung (fold = 1.48 and p = 7.8 × 10−6 for TssFlnkU, fold = 1.33 and p = 4.0 × 10−2 for TssFlnkD; Figure 3A). A recent study has uncovered that the m6A in non-coding RNA can regulate the nearby DNAme and subsequently modulate chromatin state and gene transcription; this enrichment result further indicated that m6A in mRNA may also involve in such regulations.65 Accordingly, these DNAme loci that are modulated by m6A are underrepresented in repressed regions including heterochromatin (fold = 0.70, p = 3.9 × 10−6) and ZNF genes and repeats (fold = 0.66, p = 2.6 × 10−3) compared to all DNAmes (Figure 3A). Notably, although m6A has been revealed to promote heterochromatin formation in mouse embryonic stem cells, the mechanism is mostly through the crosstalk between m6A and H3K9me3 rather than DNAme,28 suggesting a diverse mechanism of m6A in regulating chromatin states. In the reverse regulatory direction, where DNAme influences m6A, we found that these DNAme loci show a widespread enrichment across active regions, including TSSs, regions flanking TSSs, transcription, and enhancers, and, meanwhile, a global depletion across repressive regions (Figures 3B and S6B). This pattern of global enrichment is consistently observed in both lung and muscle tissues, although the significance levels are compromised in muscle, likely due to the relatively smaller number of sites.

Collectively, our enrichment analysis revealed that the DNAmes involved in the “DNAme-to-m6A” regulation show a global enrichment in open regions, whereas the DNAmes in the “m6A-to-DNAme” direction display a trend of specific enrichment in enhancer regions, implying a regulatory mode where m6A influences gene expression potentially through modulating enhancer activity.

We also examined the genomic distribution of m6A sites that are involved in the crosstalk between m6A and DNAme. For the m6A sites involved in m6A-to-DNAme regulation, we found no significant preference of distribution in different mRNA regions (Figures S6C and S6D). For the m6A sites in the DNAme-to-m6A regulation, we observed an enrichment in the transcription terminal site (TTS) region (fold = 1.14, p = 0.012) and a depletion in coding sequence (CDS; fold = 0.98, p = 0.013) for lung (Figure S6E) and an enrichment in the TTS region for muscle (fold = 1.31, p = 0.046; Figure S6F). Therefore, the robust enrichment in the TTS region in both tissues suggested that DNAme might preferentially influence the m6A sites in the 3′ end of the mRNA molecules.

Potential regulator pairs mediating the crosstalk between m6A and epigenome

Chemical modifications to RNA, DNA, and histone are installed and removed by the corresponding writers and erasers, respectively, and are recognized and modulated by their reader proteins. Therefore, we envision that the mechanism underlying the crosstalk between m6A and epigenome is likely dependent on the interactions between relevant regulators. In fact, recent studies have indicated the mechanistic roles of the readers/writers of m6A, DNAme, and histone modification in mediating such crosstalk.27,28,29,30,31,32,33 For instance, studies have identified a CFL1-METTL3-m6A-YTHDC2-MLL1 regulatory axis underlying m6A-H3K4me3 crosstalk and an m6A-YTHDC2-TET1 axis underlying m6A-DNAme crosstalk.65,66

Despite these mechanisms exemplified, the regulators involved in the crosstalk between m6A and H3K27ac/DNAme still lack a systematic investigation. Considering the various novel RBPs recently reported to potentially contribute to m6A manipulation,21,22 we next sought to prioritize RBP-TF pairs that potentially mediate the crosstalk by integrating the binding site enrichment and protein-protein interaction (PPI) maps (STAR Methods).

We first carried out enrichment analyses separately for the m6A and DNAme loci predicted to exhibit crosstalk effects by utilizing the data from the enhanced crosslinking and immunoprecipitation as well as the chromatin immunoprecipitation data available from ENCODE,64,67 respectively. For the m6A-to-DNAme crosstalk pairs in lung, we identified 5 RBPs enriched for m6A, including YTHDF2, YTHDC2, and ATXN2, and 6 TFs for DNAme loci, such as AGO1, AGO2, and POLR2A (Table S6). For the reverse direction of DNAme-to-m6A pairs, we characterized 43 TFs enriched for DNAme and 29 RBPs for m6A loci (Table S6). Notably, we observed an enrichment of YTHDC2 binding for the m6A sites that potentially regulate DNAme, aligning with the m6A-YTHDC2-TET1 regulatory axis reported by a recent study (Figure S7A).65 By further utilizing the PPI map curated by STRING,68 we prioritized the TFs that show interaction evidence with known m6A writers or erasers, resulting in 13 TF-m6A enzyme pairs that potentially mediate the regulation of DNAme on m6A, including CBFA2T3-RBM15, DNMT1-METTL3, DNMT1-ALKBH1, and POLR2A-METTL3/METTL14/WTAP (Figure S7A; Table S7), expanding the potential regulators involved in the crosstalk between m6A and DNAme.

Considering the complexity of the m6A and DNAme dynamics, we further searched for potential regulator pairs without restricting to known writers or erasers.69,70 We therefore calculated the frequency of the binding pairs for each potential RBP-TF combination that corresponds to the m6A-DNAme regulatory pairs. For the m6A-to-DNAme regulation in lung, the top RBP-TF pairs include ATXN2-POLR2A, PTBP1-POLR2A, ATXN2-CTCF, and YTHDF2-POLR2A. Among these pairs, the PPIs of PTBP1-POLR2A, ATXN2-AGO2/1, and YTHDF2-AGO2 are curated by the STRING database with experimental evidence (Figure S7B; Table S7). Likewise, we also identified TF-RBP pairs that may contribute to the regulation from DNAme to m6A in lung, including POLR2A-HNRNPC, CTCF-HNRNPC, POLR2A-YTHDF2, and POLR2A-PTBP1 (Figure S7C). For the identified regulator pairs, we further performed in silico validations to assess the enrichment of the overlaps between RBP and TF binding sites in the context of mediating the crosstalk between m6A and DNAme (STAR Methods). We observed significant enrichments of binding site pairing for 25 of the top 30 (83.3%) RBP-TF pairs underlying m6A-to-DNAme regulation and 24 (80%) of the TF-RBP pairs underlying DNAme-to-m6A regulation (Figures S7D and S7E). Additionally, we performed experimental validations for three regulators that we predicted, including ATXN2, YTHDF2, and RAD21 (STAR Methods). Upon the knockdown of ATXN2 and YTHDF2, which were predicted to regulate DNAme via interactions with multiple TFs, we observed significant decreases in DNAme levels in a predicted target, cg06491548 (Figure S7F). Likewise, we found a significant increase of the modification level at a targeting m6A locus upon the depletion of RAD21, which is a TF predicted to regulate m6A via interacting with HNRNPC and CELF2 (Figure S7G).

H3K27ac-m6A crosstalk underlying complex diseases

Genetics studies have separately revealed the vital roles of m6A and histone modification as bridges to connect genetic variants to human diseases.14,21,22,50,71,72 More recently, accumulating evidence has demonstrated that the crosstalk between m6A and histone plays essential roles in regulating chromatin states and subsequently modulating important biological processes, including cell stemness, differentiation, and drug resistance.27,28,29,30,32,33 We therefore asked whether the crosstalk between m6A and H3K27ac can provide an in-depth mechanism for GWAS variants, using a multi-trait colocalization strategy that integrates the genetic associations across the three layers (H3K27ac, m6A, and complex diseases).73 Here, we mainly focused on the crosstalk between m6A and H3K27ac that we identified in the brain given its relatively higher statistical power.

We identified 15 m6A-to-H3K27ac regulatory pairs in brain that showed multi-colocalization effects with genetic loci in 20 GWAS traits and 14 H3K27ac-to-m6A pairs in 30 GWAS traits (Tables S8 and S9; Figure 4A). Among them, 23 (33.3% of 69) GWAS loci are of sub-threshold significance (5 × 10−8 < p < 1 × 10−4), including the SGIP1 locus associated with neuroticism and the MAPK8IP1P2 locus linked to multiple brain disorders (Figure 4A). Albeit not reaching genome-wide significance (p < 5 × 10−8), sub-threshold GWAS loci could potentially contribute to the etiology and heritability of diseases.74,75 Our results provide mechanistic annotations for these GWAS loci at the layers of m6A and H3K27ac and, meanwhile, increase the confidence for these sub-threshold GWAS loci. Indeed, studies have shown that the sub-threshold GWAS gene SGIP1 plays an essential role underlying emotionality and mood phenotypes and that MAPK8IP1P2 is highly expressed in brain tissues and involved in sleep duration.76,77,78 As an example, we show that the SNP rs2271397 associated with bipolar disorder is the lead variant mediating the crosstalk between H3K27ac and a m6A modification within the LINGO1 gene (Figure 4B), which plays an essential role in the central nervous system and is relevant to neurodegeneration.79,80 We observed a posterior probability of 0.86 for this variant from the multi-trait colocalization across m6A-QTLs, haQTLs, and bipolar disorder GWAS loci in brain. These findings together suggest a regulatory circuit from m6A to H3K27ac that may link the risk of bipolar disorder to this locus (Figure 4B). Additionally, we also performed multi-trait colocalization incorporating these three genetic layers in lung, muscle, and heart tissues and identified potential contributions of m6A-H3K27ac crosstalk regulation in tissue-relevant GWAS traits, such as allergic and respiratory diseases in lung and cardiovascular disease in muscle and heart (Figure S8A).Figure 4 m6A-H3K27ac crosstalk explains the genetic risks of human diseases

(A) The genomic position (y axis) and GWAS −log10p value (x axis) of the lead SNPs linking m6A-H3K27ac bidirectional crosstalk to human traits/diseases from multi-trait colocalization analysis in the brain. For SNPs related to multiple traits/diseases, the −log10p value corresponds to the most significant GWAS locus. The lead SNPs, annotated genes, and traits/diseases are described in text, with red, yellow, and gray representing the directions of m6A-to-H3K27ac, H3K27ac-to-m6A, and bidirectional, respectively.

(B) An example of m6A-to-H3K27ac crosstalk in the LINGO1 gene in brain tissue, elucidating its potential role in bipolar disorder. The lead SNP rs2271397 is the shared causal variant between the bipolar disorder GWAS and m6A-to-H3K27ac layer, with a posterior probability of 0.86. The box shows the 25th–75th percentiles, the line shows the median, and the whiskers show 1.5× IQR.

See also Figure S8 and Tables S8 and S9.

Crosstalk between m6A and DNAme contributes to human complex diseases

The crosstalk between m6A and DNAme has also been revealed to play essential roles in multiple key biological processes.25,26,31,65,81 We therefore investigated whether such crosstalk contributes to genetic diseases in lung and muscle tissues, where we have the m6A-DNAme crosstalk pairs available for carrying out integrative genetic analysis.

We first evaluated the global heritability enrichments of m6A-DNAme crosstalk for GWAS loci using stratified linkage disequilibrium score regression in lung and muscle82,83 (STAR Methods; Tables S8 and S10). In lung, the significant DNAme-to-m6A regulatory pairs were enriched in 9 lung-relevant traits, including asthma, chest pain, wheeze, and multiple other respiratory-related traits (nominal p < 0.05; Figure 5A). In the reverse direction of the m6A-to-DNAme pairs, the only enriched GWAS trait is forced vital capacity (FVC) (Figure S8B). In muscle, the m6A-to-DNAme regulatome is enriched for injury of muscle and tendon, and the DNAme-to-m6A regulatome is enriched for tense, sore, and aching muscles (Figures S8C and S8D). The relatively fewer GWAS traits enriched in muscle are likely due to the smaller number of significant MR pairs we identified (4,076 in muscle vs. 62,432 in lung). We next focused on the specific GWAS loci that are explainable by the crosstalk between m6A and DNAme. We mapped 1,767 and 249 loci that show strong multi-colocalization posterior probabilities for the crosstalk pairs in lung and muscle, respectively (Tables S11 and S12).Figure 5 Interpretation of human diseases through m6A-DNAme bidirectional crosstalk

(A) Enrichment of DNAme-to-m6A crosstalk in the lung against lung-related traits/diseases by stratified linkage disequilibrium score regression (S-LDSC). The enrichment is calculated as the ratio between the proportion of heritability and the proportion of SNPs reported by S-LDSC. Enrichment p values were from S-LDSC via Z score calculation. Error bars represent the standard error of the S-LDSC enrichment.

(B) The genomic position (y axis) and GWAS −log10p value (x axis) of the lead SNPs linking m6A-DNAme bidirectional crosstalk to human traits/diseases from multi-trait colocalization analysis in the lung. For SNPs related to multiple traits/diseases, the −log10p value corresponds to the most significant GWAS locus. The details of annotated text adjacent to the lead SNPs are similar to that in Figure 4A, with red, yellow, and gray representing the direction of m6A-to-DNAme, DNAme-to-m6A, and bidirectional, respectively.

(C) An example of DNAme-to-m6A crosstalk in the gene ACBD3-AS1 in lung tissue, which shows genetic colocalization with asthma, with the lead SNP rs3768410 and a posterior probability of 0.94. The box shows the 25th–75th percentiles, the line shows the median, and the whiskers show 1.5× IQR.

(D) An example of DNAme-to-m6A crosstalk in the gene TYRO3 in muscle tissue, which shows genetic colocalization with diastolic blood pressure, with the lead SNP rs7174099 and a posterior probability of 0.71.

See also Figure S8 and Tables S8, S10, S11, and S12.

In lung tissue, we identified 372 disease loci that are potentially dependent on m6A-to-DNAme regulation, delineating the mechanisms for GWAS loci associated with peak expiratory flow (PEF), chest pain, shortness of breath, and more (Figure 5B; Table S11). The SNP rs10889792 associated with PEF appears to mediate the regulation from m6A to DNAme in the SCMH1 locus, a polycomb protein that has been implicated in lung function.84 As for the DNAme-to-m6A pairs, we uncovered 1,729 GWAS loci associated with asthma, lung-related functions such as FEV1, FVC, and PEF, and other allergic diseases (Figure 5B; Table S11). For instance, the SNP rs3768410 is the lead variant where we observed a high posterior probability (0.94) from multi-colocalization analysis across asthma GWASs, DNAme (cg04213163), and m6A modification (ACBD3-AS1 gene) (Figure 5C). We also noted several disease loci that were colocalized with bidirectional regulatory pairs, and the involved disease phenotypes included chest pain, wheeze, and FEV1/FVC, implying a complex regulatory mode for these disease variants (Figure 5B). The FEV1/FVC-associated SNP rs738628 links the bidirectional crosstalk between DNAme and m6A in EP300, which encodes a histone acetyltransferase, suggesting potentially a more complex regulation in this locus.85,86 In addition to lung-related diseases, the m6A-DNAme crosstalk identified in lung shows colocalization effects with other GWAS traits, including schizophrenia, depression, coronary heart disease, and more. This is likely attributed to the tissue-sharing crosstalk between m6A and DNAme that is conserved across different tissues.

In muscle tissue, we uncovered 58 genetic loci that show strong multi-colocalization effects with the m6A-to-DNAme regulatory pairs, explaining GWAS traits such as high blood pressure and pulse rate (Figure S8E; Table S12). One example is the SNP rs56104944, which has been previously reported to contribute to hypertension by regulating m6A in the HSPA4 gene that encodes a heat shock protein.21 Here, we further suggested a more comprehensive mechanism that additionally involves the DNAme level nearby. We also noticed that the variant rs12866090, associated with atrial fibrillation, likely mediates the crosstalk between m6A and DNAme in CUL4A. This gene encodes a core protein of the E3 ubiquitin ligase complex and has been reported to protect oxidative-induced cardiomyocyte apoptosis.87 In the opposite regulatory direction, we identified 209 genetic loci that could be explained by the mechanism of the DNAme-to-m6A regulation path (Figure S8E; Table S12). For example, the hand grip strength-associated SNP rs874885 is indicated to mediate the DNAme-to-m6A regulation in the LDB1 gene, which has been revealed to play a central role in modulating cardiac lineage differentiation.88 rs7174099 is a diastolic blood pressure-associated variant that shows a multi-colocalization effect with DNAme and the m6A level in the gene TYRO3 (posterior probability = 0.71; Figure 5D), which encodes a tyrosine-kinase receptor that plays a key role in endothelial-related hypertension and multiple aspects of cardiovascular pathology.89,90 The integrative analysis of m6A-DNAme crosstalk and GWAS variants collectively provides a more complete mechanistic path that bridges two molecular layers at the RNA and DNAme levels with human complex traits and diseases.

Discussion

In the current study, we identify genome-wide crosstalk between m6A and H3K27ac/DNAme, providing comprehensive cross-layer interaction maps between epitranscriptomic and epigenomic marks. With the crosstalk maps established, we dive into the tissue specificity and the genomic distribution preference of the crosstalk loci, as well as characterizing regulator pairs that potentially mediate such crosstalk. We further integrate the crosstalk maps between m6A and the two epigenomic modifications with GWASs, bridging the multi-layer mechanistic gap for human disease-associated genetic variants.

The existing QTL-based disease interpretation has been largely centered on eQTLs, with an underlying assumption that most disease variants exert their effects through altering gene expression.91 The QTL-based integrative analyses have therefore been primarily focused on eQTLs in contrast to other QTL types, including protein QTLs, caQTLs, haQTLs, and more.9,10 Despite these biologically intuitive integrations with eQTLs, the crosstalk between other QTL types has been largely neglected, causing potential mechanistic gaps in interpreting disease variants. Recent studies have reported the crosstalk between m6A, the most prevalent and functionally essential epitranscriptomic mark, and multiple epigenomic layers, including histone modifications, DNAme, and chromatin accessibility.23,24,25,26,27,28,29,30,31,32,33 Such crosstalk plays vital roles in various fundamental biological processes, including maintaining stem cell self-renewal and differentiation, driving cancer progression, and regulating developmental processes. These findings strongly demonstrate the presence and significance of such regulatory modes that involve multiple epigenetic layers. Therefore, it is of great need and importance to leverage the QTL resources for investigating the crosstalk between m6A and epigenomic modifications.

To this end, we establish systematic regulatory maps with causal directions for the crosstalk between m6A and the epigenomic layers of DNAme and H3K27ac by utilizing the genetic variations of these three molecular traits. Existing studies on the crosstalk were exclusively focused on a single regulatory direction across the whole genome, reporting either m6A-to-epigenome regulation, or vice versa.23,24,25,26,27,28,29,30,31,32,33 Our study identifies the crosstalk in each specific locus using the corresponding QTLs, therefore enabling a locus-specific identification of crosstalk, with the regulatory direction unambiguously delineated. We achieve such locus-specific regulation identification by employing a bidirectional MR framework, which was further supported by heterogeneity test and genetic colocalization. With the QTL resources of m6A, DNAme, and H3K27ac across multiple tissues, our approach allows for the interrogation and comparison of the crosstalk between m6A and epigenomic layers in different contexts. Moreover, a unique and apparent advantage of the MR-based approach is that we can directly link the m6A-epigenome crosstalk to human genetic diseases, therefore building up comprehensive regulatory circuits underlying these diseases. Of note, while in this study, we focused on the m6A-epigenome crosstalk in the GTEx cohort given the data availability, this strategy can be extended to other molecular traits in different populations where full summary statistics are available.

Our results shed light on the regional specificity of the crosstalk between m6A and DNAme by examining the enrichment of the DNAme loci showing crosstalk effects in the chromatin states. Notably, to ensure that we capture the specific distribution preference of the m6A-interacted DNAme loci rather than the general distribution property of DNAme, we used all the DNAme loci that were tested for crosstalk effects as background. We found that the distribution enrichments for chromatin states were different between the DNAme in the forward and reverse regulatory directions with m6A modification. While the crosstalk loci in both directions show a global enrichment in the open chromatin regions, the DNAme loci that are regulated by m6A tend to be strongly enriched in the enhancer regions. Moreover, the regional distributions of the crosstalk loci are highly consistent between different tissues for both regulatory directions, indicating the robustness of such crosstalk-based regulations. Despite these global enrichments observed, the loci with m6A and epigenome crosstalk are widespread across the genome, indicating the sensitivity of our approach in identifying locus-specific crosstalk.

Recent studies investigating the crosstalk between m6A and epigenomes have revealed how these epigenetic modifications from different layers cooperatively regulate multiple processes, particularly stem cell differentiation for both human and mouse.23,24,25,26 However, whether and how m6A-epigenome crosstalk mechanisms contribute to human genetic disease have not been investigated to the best of our knowledge. In this study, we forge a new path to investigate the crosstalk between m6A and epigenomes by using genetic variation across molecular traits and disease phenotypes. The interpretation of disease variants using QTLs has been primarily centered on mRNA expression, with very few studies diving into the regulatory paths that cover multiple epigenetic layers. Our study, by design, starts from the relationship between m6A and epigenomic modifications and subsequently moves into GWAS variant interpretations based on the established molecular crosstalk foundations. This allows us to construct a comprehensive regulatory circuit encompassing multiple molecular features. As a result, we systematically elucidated genetic variants associated with various complex diseases, including neurodegeneration, schizophrenia, asthma, congenital heart disease, and more. Our results collectively established the crucial roles of m6A-epigenome crosstalk in complex human diseases, extending the functional significance of such crosstalk beyond its known roles in modulating stem cell differentiation, drug resistance, and cancer progression.

Overall, the utilization of genetic approaches and resources for interrogating the crosstalk spanning different epigenetic layers provides a mechanistic view toward comprehensively understanding the regulatory flows from genetic variants to diseases. The resulting mechanistic insights can facilitate the search for potential diagnostic and therapeutic targets for human complex diseases at multiple regulatory layers.

Limitations of the study

Despite these advancements, our study has the following limitations. First, while we assembled the largest QTLs for m6A and epigenomic biomarkers, the sensitivity of crosstalk identification is constrained by the sample sizes available for each dataset, especially for the m6A-QTL given the experimental difficulty for m6A profiling.21,22 We believe that there will be substantially more loci uncovered as the sample size and the statistical power increase for the QTLs of m6A, DNAme, and histone modification. The increased statistical power will also facilitate the interpretation of more GWAS variants. Second, in our MR framework, we investigated the genetically predicted associations between m6A and epigenome, where the crosstalk independent of genetic regulation may not be captured. Therefore, it will be helpful to integrate both genetic-based approaches and functional experiment validations, such as the typical writer/eraser depletion-based assays, for a comprehensive mapping of crosstalk and identification of underlying mechanisms. Third, the tissues/cell types and the epigenomic modification types are available for five primary human tissues at the current stage. As the QTL resources of different molecular traits continue to accumulate in diverse tissues/cell types, we will gain more comprehensive knowledge of the specificity of the crosstalk between m6A and various types of epigenomic regulation. Lastly, the current m6A-QTLs were confined to mRNA molecules due to the limitation of the methylated RNA immunoprecipitation sequencing methods utilized; therefore, future research is warranted to interrogate the m6A sites in the non-mRNA types, such as enhancer RNA and long non-coding RNA. While recent studies have reported the crosstalk between epigenome and m6A in various non-coding RNA species, we envision that an optimized m6A profiling method that captures a full spectrum of m6A-QTLs will help construct a more comprehensive crosstalk map between m6A and epigenomes.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Antibodies	
	
rabbit polyclonal anti-m6A	ABclonal	Cat# A17924; RRID: AB_2770239	
	
Chemicals, peptides, and recombinant proteins	
	
Dulbecco’s modified eagles’s medium	VivaCell	Cat# C3113-0500	
fetal bovine serum	Gibco	Cat# 10270-106	
RPMI-1640	VivaCell	Cat# C3010-0500	
Puromycin	Gibco	Cat# A11138-03	
1M Tris-HCl 7.5	Invitrogen	Cat# 15567-027	
IGEPAL CA-630	MP Biomedicals	Cat# 198596	
Murine RNase Inhibitor	Vazyme	Cat# R301	
5M NaCl	Invitrogen	Cat# AM9760G	
TRIzon reagent	CWBIO	Cat# CW0580S	
	
Critical commercial assays	
	
Lipofectamine 3000 reagent	Thermo Fisher Scientific	Cat# L3000015	
FastPure Cell/Tissue Total RNA Isolation Kit	Vazyme	Cat# RC101-01	
One Step TB Green PrimeScrip RT-PCR Kit	Takara	Cat# RR066	
Universal Genomic DNA Kit	CWBIO	Cat# CW2298M	
EZ DNA Methylation Kit	Zymo Research	Cat# D5002	
ChamQ Universal SYBR qPCR Master Mix	Vazyme	Cat# Q711	
Dynabeads mRNA DIRECT purification kit	Thermo Fisher Scientific	Cat# 61011	
NEBNext Magnesium RNA Fragmentation Module	NEB	Cat# E6150S	
RNA Clean & Concentrator-5 kit	Zymo Research	Cat# R1016	
Direct-zol RNA Microprep Kit	Zymo Research	Cat# R2060	
	
Deposited data	
	
m6A-QTL summary statistics for brain, lung, muscle, and heart tissues from the eGTEx project	Xiong et al.21	https://compbio.mit.edu/m6AQTLs/m6A.QTL.Xiong.etal/SumStat/	
haQTL summary statistics for brain, lung, muscle, and heart tissues from the eGTEx project	Hou et al.50	https://zenodo.org/records/7992724	
mQTL summary statistics for lung and muscle tissues from the eGTEx project	Oliva et al.51	https://gtexportal.org/home/datasets	
Code used in this study	This manuscript	https://doi.org/10.5281/zenodo.11481934	
	
Experimental models: Cell lines	
	
Human: A549	ATCC(CRL-185)	N/A	
Human: HEK293T	ATCC(CRL-11268)	N/A	
	
Oligonucleotides	
	
Oligonucleotides used in this study for experimental validation	Table S13	N/A	
	
Recombinant DNA	
	
Plasmid: pLKO.1	Moffat et al.92	Addgene Plasmid # 10878	
Plasmid: psPAX2	Trono Lab Packaging and Envelope Plasmids (unpublished)	Addgene Plasmid # 12260	
Plasmid: pMD2.G	Trono Lab Packaging and Envelope Plasmids (unpublished)	Addgene Plasmid # 12259	
	
Software and algorithms	
	
TwoSampleMR v0.5.6	Hemani et al.93	https://mrcieu.github.io/TwoSampleMR/	
SMR	Zhu et al.53	https://yanglab.westlake.edu.cn/software/smr/	
bedtools v2.31.0	Quinlan et al.94	https://bedtools.readthedocs.io/en/latest/	
coloc v5.2.2	Giambartolomei et al.54	https://chr1swallace.github.io/coloc/	
S-LDSC v1.0.1	Bulik-Sullivan et al.82	https://github.com/bulik/ldsc	
moloc v0.1.0	Giambartolomei et al.73	https://github.com/clagiamba/moloc	
	
Other	
	
Protein A/G Magnetic Beads	Selleck	Cat# B23201	

Resource availability

Lead contact

Further information and requests for resources should be directed to the lead contact, Xushen Xiong (xiongxs@zju.edu.cn).

Materials availability

This study did not generate new unique reagents.

Data and code availability

Code for the data processing and statistical genetics analysis is accessible at the github repository (https://github.com/xiongxslab/LiC_et_al_MR_crosstalk) and Zenodo (https://doi.org/10.5281/zenodo.11481934). All the summary statistics data of QTLs are available publicly, with detailed information available in the key resources table and Table S1.

Experimental model and subject details

Epigenome and epitranscriptome datasets

To identify the crosstalk between epitranscriptome (m6A RNA modification) and epigenome (DNA methylation and H3K27ac histone modification), we collected the summary statistics of m6A-QTL for m6A, haQTL for H3K27ac and mQTL for DNA methylation.21,50,51 The availability of these summary statistics is detailed in Table S1. Information including ethnicity, sample size, and number of QTL loci are provided. Briefly, the summary statistics of m6A-QTLs and haQTLs are available for brain, lung, muscle, and heart tissues, and mQTLs for lung and muscle tissues. The datasets were all from European ancestry. The sample sizes of the m6A-QTL dataset range from 32 to 53, resulting in the availability of 302–539 genetically-associated m6A methylations. The sample sizes of the haQTL dataset range from 66 to 113, resulting in the availability of 537 to 2,162 genetically-associated H3K27ac loci. The sample size of the mQTL dataset is 190 for the lung and 42 for the muscle, which generated 158,503 and 14,131 genetically-associated DNA methylation loci, respectively.

Method details

Bidirectional two-sample MR for crosstalk directionality inference

We implemented a bidirectional two-sample Mendelian randomization analysis to infer the putative causal directionality between epigenome and epitranscriptome. For the MR tests from m6A to epigenome, we selected the top independent m6A-QTL(s) at each m6A locus as genetic instrumental variables (IVs). For the reverse MR tests from DNA methylation and H3K27ac to m6A, we chose the top independent mQTL(s) and haQTL(s) as IVs, respectively. To select IVs from QTL summary statistics, we excluded the SNPs with MAF lower than 0.05 to ensure the robustness of MR tests. To ensure the independence of the IVs for each locus, we applied linkage disequilibrium (LD) clumping to iteratively remove the SNPs in high LD with the lead variants (r2 < 0.01 within 100 kb window size, using GTEx version 8 reference panel95). Of note, the selection of IVs was purely based on their independence and significance as QTLs in exposure traits, without considering their significance in outcome traits.

We performed the Wald ratio MR tests for the loci with a single IV, while the inverse variance weighted (IVW) method was used for the MR tests with at least two IVs.96 To account for multiple testing, we applied the Benjamini-Hochberg correction for p values and used the threshold of FDR <0.1 to prioritize significant causal relationships. We conducted the MR tests using the R package TwoSampleMR (v0.5.6)93.

Using SMR and HEIDI tests for verification

To ensure the reliability of the crosstalk between m6A and epigenome we identified based on the two-sample MR approach, we further carried out the SMR (summary data-based Mendelian randomization) test.53 Although none of the crosstalk loci identified by the two-sample MR method reached the significance threshold (BH-adjusted p value <0.1) based on the SMR test because of the limited statistical power, we found that the effect sizes between the two-sample MR and SMR are highly consistent. This indicated the reliability of the crosstalk loci that we identified and the robustness of the effect sizes of the crosstalk.

We further performed the HEIDI (heterogeneity in dependent instruments) test, which can detect whether the associations between the epigenomic and epitranscriptomic layers are due to a shared genetic variant or a linkage model.53 We filtered out those two-sample MR findings with PHEIDI < 0.01.

Cross-tissue effect sizes consistency of bidirectional MR findings

We compared the effect sizes and the directionalities of the crosstalk loci between m6A and epigenome. For the m6A-DNAme crosstalk pairs, we carried out the cross-tissue consistency analysis using the significant m6A-DNAme MR pairs identified in lung and in muscle, respectively. To understand the potential effect caused by the power issue, we tried different p value thresholds (FDR = 0.1, 0.05, 0.01, and 0.001) in the DNAme-to-m6A MR analysis when defining the significant crosstalk loci and used them for this cross-tissue consistency analysis. A weighted regression analysis was performed to calculate the correlation coefficients by taking the standard errors of MR effect sizes in the other tissue (compared to the primary tissue) into consideration. Given the limited number of the significant m6A-H3K27ac crosstalk loci identified, we conducted the same consistency analysis for the merged interaction loci across the four tissues.

Genomic distribution enrichment of DNA methylation interacted with m6A

To examine the genomic distribution preference for the crosstalk between m6A and DNA methylation, we downloaded the ChromHMM epigenomic annotation from Roadmap.63,97 The E096 (lung) and E108 (muscle) annotations were utilized for the enrichment analysis of the crosstalk identified in lung and muscle, respectively. The intersection between DNA methylation loci and the chromatin state annotations was performed using bedtools (v2.31.0) intersect sub-command.94 The enrichment fold was calculated as the proportion of the m6A-interacted DNA methylation loci in a specific functional region over the proportion of all the MR tested DNA methylation loci annotated in this region. The p value of enrichment was calculated using the two-sided Fisher’s exact test in fisher.test function in R.

Identification of potential regulators underlying m6A and DNA methylation crosstalk

To uncover potential regulators involved in the crosstalk between m6A and DNA methylation, we identified the RBPs whose binding sites enrich for the involved m6A and TFs for the involved DNA methylation. Binding sites of 220 RBPs were obtained from the eCLIP datasets curated by the POSTAR3 database, and binding sites of 340 TFs were downloaded from the ENCODE version 2 and 3 ChIP data.64,67 For each molecular layer (m6A or DNA methylation) under each crosstalk directionality (“m6A-to-DNAme” or “DNAme-to-m6A”), the enrichment fold for each specific protein (RBP or TF) was calculated as the proportion of significant MR loci overlapped over the proportion of all the genetically-associated loci overlapped. The right-sided Fisher’s exact test was performed to examine the enrichment effect. We used BH-adjusted p value <0.05 as a threshold to define enriched RBPs for m6A loci and enriched TFs for DNA methylation.

Next, we utilized two strategies to identify potential regulator pairs in mediating the crosstalk between m6A and DNA methylation. In the first strategy, we characterized the enriched RBPs that showed protein-protein interactions (PPIs) with known writers/erasers of DNA methylation, thereby establishing the potential “m6A-to-DNAme” regulatory axis. Likewise, we recognized the enriched TFs that showed PPIs with known m6A writers/erasers to explain “DNAme-to-m6A” regulatory pairs. The PPIs were based on the curation of the STRING database68 and literature evidence. For the second strategy, we looked for RBP-TF pairs without requiring one of them to be a known enzyme, aiming to uncover novel regulator pairs potentially involved in the crosstalk between m6A and DNA methylation. For the RBPs and TFs that showed binding site enrichments in the crosstalk between m6A and DNA methylation loci that we identified above, we again used PPI information to prioritize RBP-TF pairs, which resulted in 30 RBP-TF pairs for the m6A-to-DNAme regulation, and 1,239 TF-RBP pairs for the DNAme-to-m6A regulation, respectively.

To validate the regulator pairs, we calculated the enrichment of the overlap between RBP and TF binding sites underlying the m6A and DNA methylation crosstalk using the eCLIP-seq and ChIP-seq data we described above. For each RBP-TF pair (or vice versa), we asked whether the overlap between the RBP binding sites (eCLIP-seq peaks) and TF binding sites (ChIP-seq peaks) underlying the significant m6A-to-DNAme MR pairs are enriched compared to all the overlapped site pairs used for MR tests. For each m6A and DNA methylation pair, the m6A overlapped by an RBP binding site and the DNA methylation overlapped by a TF binding site were calculated as a molecular pair mediated by the corresponding RBP and TF pair. The fold enrichment for each regulator pair is calculated as follows:Foldenrichment=#SignificantexposurelociboundbyR1+#SignificantoutcomelociboundbyR2#TestedexposurelociboundbyR1+#TestedoutcomelociboundbyR2−(#SignificantexposurelociboundbyR1+#SignificantoutcomelociboundbyR2)#Significantexposureloci+#Significantoutcomeloci#Testedexposureloci+#Testedoutcomeloci−(#Significantexposureloci+#Significantoutcomeloci)

where R1 refers to the RBP and R2 refers to the TF in a prioritized RBP-TF pair mediating the m6A-to-DNAme regulation, while R1 refers to the TF and R2 refers to the RBP in a prioritized TF-RBP pair mediating the opposite regulation. The p value of enrichment was calculated using the right-sided Fisher’s exact test.

Experimental validation of the regulators mediating the m6A and DNA methylation

HEK293T cells were cultured in Dulbecco’s modified eagles’s medium supplemented with 10% fetal bovine serum. A549 cells were cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum. shRNA targeting sequences were chosen from The RNAi Consortium collection (https://www.sigmaaldrich.cn/CN/en/semi-configurators/sirna) and cloned into pLKO.1 vector (10878, Addgene). The shRNA targeting sequences used were listed in Table S13. The shRNA lentivirus was produced using Lipofectamine 3000 reagent (L3000015, Thermo Fisher Scientific) with shRNA plasmid, psPAX2 (12260, Addgene) and pMD2.G (12259, Addgene) co-transfecting into HEK293T cells according to the manufacturer’s instructions. After 48 h, the supernatant was collected and applied to A549 cells. After 96 h, stable gene-knockdown cell lines were selected with 0.15 μg/mL puromycin for 5 days. Total RNA was extracted from cell pellet using the FastPure Cell/Tissue Total RNA Isolation Kit (RC101-01, Vazyme), and the relative target gene expression level was then evaluated with qRT-PCR using One Step TB Green PrimeScrip RT-PCR Kit (RR066, Takara). Primers used for qRT-PCR were listed in Table S13.

The following steps were used for quantifying DNA methylation level. Genomic DNA was extracted from cell pellet using the Universal Genomic DNA Kit (CW2298M, CWBIO) according to the manufacturer’s instructions. Then 2 μg genomic DNA was bisulfite-converted using EZ DNA Methylation Kit (D5002, Zymo Research) according to the manufacturer’s instructions and eluted with 40 μL H2O. The DNA methylation level was subsequently calculated based on the difference of Ct values of two qPCR reactions using ChamQ Universal SYBR qPCR Master Mix (Q711, Vazyme), with primers of one specific for methylated and another one for unmethylated DNA. Primers used for qPCR were listed in Table S13.

The following steps were used for quantifying m6A level. Total RNA was extracted from cell pellet using the FastPure Cell/Tissue Total RNA Isolation Kit (RC101-01, Vazyme) according to the manufacturer’s instructions. Further, 50 μg total RNA was purified with the Dynabeads mRNA DIRECT purification kit (61011, Thermo Fisher Scientific) and 1 μg mRNA from that was fragmented to 200 nucleotide fragments by NEBNext Magnesium RNA Fragmentation Module (E6150S, NEB) at 94°C for 5 min and purified with RNA Clean & Concentrator-5 kit (R1016, Zymo Research). 1/10 of fragmented RNA was kept as input. For m6A immunoprecipitation, 0.25 mg Protein A/G Magnetic Beads (B23201, Selleck) was linked with 2 μL rabbit anti-m6A antibody (ABclonal, A17924, 1 mg/mL) for 4 h at 4°C. Then, approximately 9/10 of fragmented RNA was bound to antibody-coupled beads in 200 μL of Binding buffer (150 mM NaCl, 10 mM Tris-HCl 7.5 and 0.1% IGEPAL CA-630) supplied with 1U/μL Murine RNase Inhibitor (R301, Vazyme) for 4 h at 4°C. The beads-antibody-RNA was washed twice successively by Binding buffer, low-salt Wash buffer (50 mM NaCl, 10 mM Tris-HCl 7.5 and 0.1% IGEPAL CA-630) and high-salt Wash buffer (500 mM NaCl, 10 mM Tris-HCl 7.5 and 0.1% IGEPAL CA-630) supplied with 0.1U/μL Murine RNase Inhibitor and then extracted using Direct-zol RNA Microprep Kit (R2060, Zymo Research) with 300 μL TRIzon reagent. Locus-specific m6A-IP-qPCR was performed in IP and input samples to determine the relative m6A level using One Step TB Green PrimeScrip RT-PCR Kit (RR066, Takara). Primers used for qPCR were listed in Table S13.

Colocalization between molecular QTLs

We carried out colocalization analyses between m6A-QTL and haQTL and between m6A-QTL and mQTL using a Bayesian-based method in R package coloc (v5.2.2)54. For each two-sample MR tested molecular pair, beta (effect size) and varbeta (the square of standard error) were used to calculate the posterior probabilities under five different assumptions (namely H0 to H4). The H0 represents no causal variant identified, H1 and H2 represent a causal effect for only one of the two traits, H3 represents two distinct causal variants, and H4 represents colocalized genetic effects. While both two traits we used were pre-selected significant genetic loci, we relied on the posterior probabilities for assumptions H3 (PP3) and H4 (PP4) to define significant colocalization events. We selected results with PP4 > 0.1 and PP4/(PP3 + PP4) > 0.75 as the significant colocalization molecular pairs.

GWAS summary statistics collection

We collected GWAS summary statistics of brain-, lung-, muscle-, and heart-relevant diseases or complex traits from GWAS ATLAS and IEU OpenGWAS database for the multi-trait colocalization and the heritability enrichment analysis.1,98 The detailed information about the accession number, the ancestry, description of the complex trait, sample size, reference of the GWAS studies were summarized in Table S8. All variants from GWAS datasets were converted to (if initially not) GRCh38 (hg38) coordinates in order to be consistent with our QTL datasets. Missing MAF were imputed from the 1000 Genomes Project phase 3 (1KGP v3) reference panel.99 We unified the effect allele and other allele in both GWAS datasets collected and our QTL datasets, and reversed the signs of effect sizes for the inconsistent loci if needed.

GWAS enrichment analysis

We carried out GWAS heritability enrichment analysis using S-LDSC (v1.0.1)82,83 for the genomic loci that mediate crosstalk between m6A and DNA methylation. We did not perform enrichment analysis for the crosstalk between m6A and H3K27ac given the limited number of loci identified. Baseline model version 1.1 was downloaded from LDSC and utilized for the calculation of heritability enrichment. For the DNAme-to-m6A regulatory loci, we used the DNA methylation loci for the GWAS heritability enrichment analysis. For the m6A-to-DNAme loci, we used the m6A loci for the GWAS heritability enrichment analysis. We performed GWAS enrichment analysis against lung-related diseases and respiratory traits for the crosstalk identified in lung, and muscle-related traits for the crosstalk identified in muscle (see Table S10).

Multi-trait colocalization across molecular QTLs and GWAS diseases

The epitranscriptome-epigenome crosstalk mechanism discovered above (detected via both bidirectional MR and colocalization method) can be further integrated with GWAS datasets to prioritize and interpret genetic risk loci. To this end, we performed multiple trait colocalization analysis across the two molecular-layer QTLs and GWAS datasets, using the GWAS summary statistics we collected. Multi-trait colocalization analyses were carried out using the R package moloc (v0.1.0), an extension of the Bayesian-based colocalization method adopted for multiple traits.73 Considering the relatively smaller sample size for our QTL datasets, effect size (beta) and standard error (se) were used as inputs for molecular QTLs, and p value of summary statistics was used as the input for GWAS. Accounting for the different amounts of loci identified between the crosstalk of m6A-H3K27ac and m6A-DNAme, among all the posterior probabilities of 15 configurations, we chose a threshold of pp_abc >0.1 and pp_abc >0.25 to define the potential sharing variants of the two molecular layers and GWAS data for these two crosstalk modalities above, respectively.

Quantification and statistical analysis

All of the quantitative and statistical methods, strategies, and analyses are described in the relevant sections of the method details or in the relevant figure legends.

Supplemental information

Document S1. Figures S1–S8

Table S1. Details about availability, ethnicity, sample size, and number of QTLs of m6A-QTLs, haQTLs, and mQTLs used in this study, related to Figure 1 and STAR Methods

Table S2. Summary of significant bidirectional MR and colocalization loci between m6A and epigenome, related to Figures 2, S1, and S3

Table S3. Crosstalk pairs between m6A and DNAme in lung, identified by bidirectional MR and colocalization analysis, related to Figures 2 and S1

Table S4. Crosstalk pairs between m6A and DNAme in muscle, identified by bidirectional MR and colocalization analysis, related to Figures 2 and S1

Table S5. Crosstalk pairs between m6A and H3K27ac in brain, lung, muscle, and heart, identified by bidirectional MR and colocalization analysis, related to Figures 2 and S3

Table S6. Enrichment analysis of RBP or TF binding sites in crosstalk-involved m6A or DNAme loci, respectively, in the lung, related to Figures 3 and S7

Table S7. Potential RBP-TF pairs underlying the crosstalk between m6A and DNAme in the lung, related to Figures 3 and S7

Table S8. Summary of the GWASs collected, which are relevant to the tissue types in this study, related to Figures 1, 4, 5, and S8

Table S9. GWAS loci of tissue-relevant human traits/diseases interpreted by the crosstalk between m6A and H3K27ac in brain, lung, muscle, and heart, identified by multi-trait colocalization analysis, related to Figures 4 and S8

Table S10. S-LDSC results of m6A-DNAme crosstalk in lung- and muscle-relevant GWAS traits/diseases, related to Figures 5 and S8

Table S11. GWAS loci of tissue-relevant human traits/diseases interpreted by the crosstalk between m6A and DNAme in lung, identified by multi-trait colocalization analysis, related to Figure 5

Table S12. GWAS loci of tissue-relevant human traits/diseases interpreted by the crosstalk between m6A and DNAme in muscle, identified by multi-trait colocalization analysis, related to Figures 5 and S8

Table S13. The shRNA targeting sequences and primers used for experimental validation, related to STAR Methods and Figure S7

Document S2. Article plus supplemental information

Acknowledgments

We thank the members of the Xiong lab for discussion and suggestions throughout the project. We are grateful for the support from the core facilities and computing platform of Liangzhu Laboratory at Zhejiang University. This work was supported by the 10.13039/501100001809 National Natural Science Foundation of China (32370609 and 92353301 to X.X.), the 10.13039/501100012166 National Key R&D Program of China (2023YFA1800700 to X.X.), and funding from the Liangzhu Laboratory at Zhejiang University and the State Key Laboratory of Transvascular Implantation Devices at Zhejiang University.

Author contributions

This study was designed by C.L., X. Han, and X.X. and directed and coordinated by X. Han and X.X. C.L. carried out the analyses with help from J.N., J.H., L.H., X. Hu, and M.K. and under the supervision of X. Han and X.X. K.C., Q.F., and S.S. performed experimental validations with help from Y.Y., X.L., and J.L. All authors participated in the discussion of the project. C.L., X. Han, and X.X. wrote the manuscript.

Declaration of interests

The authors declare no competing interests.

Supplemental information can be found online at https://doi.org/10.1016/j.xgen.2024.100605.
==== Refs
References

1 Watanabe K. Stringer S. Frei O. Umićević Mirkov M. de Leeuw C. Polderman T.J.C. van der Sluis S. Andreassen O.A. Neale B.M. Posthuma D. A global overview of pleiotropy and genetic architecture in complex traits Nat. Genet. 51 2019 1339 1348 31427789
2 Wang J. Huang D. Zhou Y. Yao H. Liu H. Zhai S. Wu C. Zheng Z. Zhao K. Wang Z. CAUSALdb: a database for disease/trait causal variants identified using summary statistics of genome-wide association studies Nucleic Acids Res. 48 2020 D807 D816 31691819
3 Sollis E. Mosaku A. Abid A. Buniello A. Cerezo M. Gil L. Groza T. Güneş O. Hall P. Hayhurst J. The NHGRI-EBI GWAS Catalog: knowledgebase and deposition resource Nucleic Acids Res. 51 2023 D977 D985 36350656
4 Ward L.D. Kellis M. Interpreting noncoding genetic variation in complex traits and human disease Nat. Biotechnol. 30 2012 1095 1106 23138309
5 Tak Y.G. Farnham P.J. Making sense of GWAS: using epigenomics and genome engineering to understand the functional relevance of SNPs in non-coding regions of the human genome Epigenet. Chromatin 8 2015 57
6 Uffelmann E. Huang Q.Q. Munung N.S. de Vries J. Okada Y. Martin A.R. Martin H.C. Lappalainen T. Posthuma D. Genome-wide association studies Nat. Rev. Methods Primers 1 2021 59 10.1038/s43586-021-00056-9
7 Hardison R.C. Blobel G.A. GWAS to therapy by genome edits? Science 342 2013 206 207 24115432
8 Abdellaoui A. Yengo L. Verweij K.J.H. Visscher P.M. 15 years of GWAS discovery: Realizing the promise Am. J. Hum. Genet. 110 2023 179 194 36634672
9 Umans B.D. Battle A. Gilad Y. Where are the disease-associated eQTLs? Trends Genet. 37 2021 109 124 32912663
10 Albert F.W. Kruglyak L. The role of regulatory variation in complex traits and disease Nat. Rev. Genet. 16 2015 197 212 25707927
11 eGTEx Project Enhancing GTEx by bridging the gaps between genotype, gene expression, and disease Nat. Genet. 49 2017 1664 1670 29019975
12 Grundberg E. Small K.S. Hedman Å.K. Nica A.C. Buil A. Keildson S. Bell J.T. Yang T.-P. Meduri E. Barrett A. Mapping cis- and trans-regulatory effects across multiple tissues in twins Nat. Genet. 44 2012 1084 1089 22941192
13 GTEx Consortium The Genotype-Tissue Expression (GTEx) project Nat. Genet. 45 2013 580 585 23715323
14 Sun W. Poschmann J. Cruz-Herrera Del Rosario R. Parikshak N.N. Hajan H.S. Kumar V. Ramasamy R. Belgard T.G. Elanggovan B. Wong C.C.Y. Histone acetylome-wide association study of autism spectrum disorder Cell 167 2016 1385 1397.e11 27863250
15 Wainberg M. Sinnott-Armstrong N. Mancuso N. Barbeira A.N. Knowles D.A. Golan D. Ermel R. Ruusalepp A. Quertermous T. Hao K. Opportunities and challenges for transcriptome-wide association studies Nat. Genet. 51 2019 592 599 30926968
16 Huan T. Joehanes R. Song C. Peng F. Guo Y. Mendelson M. Yao C. Liu C. Ma J. Richard M. Genome-wide identification of DNA methylation QTLs in whole blood highlights pathways for cardiovascular disease Nat. Commun. 10 2019 4267 31537805
17 GTEx ConsortiumLaboratory, Data Analysis &Coordinating Center (LDACC)—Analysis Working GroupStatistical Methods groups—Analysis Working GroupEnhancing GTEx (eGTEx) groupsNIH Common FundNIH/NCINIH/NHGRINIH/NIMHNIH/NIDABiospecimen Collection Source Site—NDRI Genetic effects on gene expression across human tissues Nature 550 2017 204 213 29022597
18 Wu L. Candille S.I. Choi Y. Xie D. Jiang L. Li-Pook-Than J. Tang H. Snyder M. Variation and genetic control of protein abundance in humans Nature 499 2013 79 82 23676674
19 Li Y.I. van de Geijn B. Raj A. Knowles D.A. Petti A.A. Golan D. Gilad Y. Pritchard J.K. RNA splicing is a primary link between genetic variation and disease Science 352 2016 600 604 27126046
20 Park E. Guo J. Shen S. Demirdjian L. Wu Y.N. Lin L. Xing Y. Population and allelic variation of A-to-I RNA editing in human transcriptomes Genome Biol. 18 2017 143 28754146
21 Xiong X. Hou L. Park Y.P. Molinie B. GTEx ConsortiumGregory R.I. Kellis M. Genetic drivers of m6A methylation in human brain, lung, heart and muscle Nat. Genet. 53 2021 1156 1165 34211177
22 Zhang Z. Luo K. Zou Z. Qiu M. Tian J. Sieh L. Shi H. Zou Y. Wang G. Morrison J. Genetic analyses support the contribution of mRNA N6-methyladenosine (m6A) modification to human disease heritability Nat. Genet. 52 2020 939 949 32601472
23 Wei J. He C. Chromatin and transcriptional regulation by reversible RNA methylation Curr. Opin. Cell Biol. 70 2021 109 115 33706173
24 Xu Z. Xie T. Sui X. Xu Y. Ji L. Zhang Y. Zhang A. Chen J. Crosstalk between histone and m6A modifications and emerging roles of m6A RNA methylation Front. Genet. 13 2022 908289
25 Kan R.L. Chen J. Sallam T. Crosstalk between epitranscriptomic and epigenetic mechanisms in gene regulation Trends Genet. 38 2022 182 193 34294427
26 Xu W. Shen H. When RNA methylation meets DNA methylation Nat. Genet. 54 2022 1261 1262 36071174
27 Liu J. Dou X. Chen C. Chen C. Liu C. Xu M.M. Zhao S. Shen B. Gao Y. Han D. He C. N6-methyladenosine of chromosome-associated regulatory RNA regulates chromatin state and transcription Science 367 2020 580 586 31949099
28 Xu W. Li J. He C. Wen J. Ma H. Rong B. Diao J. Wang L. Wang J. Wu F. METTL3 regulates heterochromatin in mouse embryonic stem cells Nature 591 2021 317 321 33505026
29 Liu J. Gao M. He J. Wu K. Lin S. Jin L. Chen Y. Liu H. Shi J. Wang X. The RNA m6A reader YTHDC1 silences retrotransposons and guards ES cell identity Nature 591 2021 322 326 33658714
30 Wei J. Yu X. Yang L. Liu X. Gao B. Huang B. Dou X. Liu J. Zou Z. Cui X.-L. FTO mediates LINE1 m6A demethylation and chromatin regulation in mESCs and mouse development Science 376 2022 968 973 35511947
31 Deng S. Zhang J. Su J. Zuo Z. Zeng L. Liu K. Zheng Y. Huang X. Bai R. Zhuang L. RNA m6A regulates transcription via DNA demethylation and chromatin accessibility Nat. Genet. 54 2022 1427 1437 36071173
32 Huang H. Weng H. Zhou K. Wu T. Zhao B.S. Sun M. Chen Z. Deng X. Xiao G. Auer F. Histone H3 trimethylation at lysine 36 guides m6A RNA modification co-transcriptionally Nature 567 2019 414 419 30867593
33 Wang J. Li Y. Wang P. Han G. Zhang T. Chang J. Yin R. Shan Y. Wen J. Xie X. Leukemogenic chromatin alterations promote AML leukemia stem cells via a KDM4C-ALKBH5-AXL signaling axis Cell Stem Cell 27 2020 81 97.e8 32402251
34 Roundtree I.A. Evans M.E. Pan T. He C. Dynamic RNA modifications in gene expression regulation Cell 169 2017 1187 1200 28622506
35 Zaccara S. Ries R.J. Jaffrey S.R. Reading, writing and erasing mRNA methylation Nat. Rev. Mol. Cell Biol. 20 2019 608 624 31520073
36 Zhao B.S. Roundtree I.A. He C. Post-transcriptional gene regulation by mRNA modifications Nat. Rev. Mol. Cell Biol. 18 2017 31 42 27808276
37 Barbieri I. Kouzarides T. Role of RNA modifications in cancer Nat. Rev. Cancer 20 2020 303 322 32300195
38 Yao B. Christian K.M. He C. Jin P. Ming G.-L. Song H. Epigenetic mechanisms in neurogenesis Nat. Rev. Neurosci. 17 2016 537 549 27334043
39 Li X. Xiong X. Yi C. Epitranscriptome sequencing technologies: decoding RNA modifications Nat. Methods 14 2016 23 31 28032622
40 Frye M. Harada B.T. Behm M. He C. RNA modifications modulate gene expression during development Science 361 2018 1346 1349 30262497
41 Huang H. Weng H. Chen J. m6A modification in coding and non-coding RNAs: roles and therapeutic implications in cancer Cancer Cell 37 2020 270 288 32183948
42 Livneh I. Moshitch-Moshkovitz S. Amariglio N. Rechavi G. Dominissini D. The m6A epitranscriptome: transcriptome plasticity in brain development and function Nat. Rev. Neurosci. 21 2020 36 51 31804615
43 Jiang X. Liu B. Nie Z. Duan L. Xiong Q. Jin Z. Yang C. Chen Y. The role of m6A modification in the biological functions and diseases Signal Transduct. Targeted Ther. 6 2021 74
44 Yang C. Hu Y. Zhou B. Bao Y. Li Z. Gong C. Yang H. Wang S. Xiao Y. The role of m6A modification in physiology and disease Cell Death Dis. 11 2020 960 33162550
45 He P.C. He C. m6A RNA methylation: from mechanisms to therapeutic potential EMBO J. 40 2021 e105977
46 Sendinc E. Shi Y. RNA m6A methylation across the transcriptome Mol. Cell. 83 2023 428 441 36736310
47 Fu Y. Dominissini D. Rechavi G. He C. Gene expression regulation mediated through reversible m6A RNA methylation Nat. Rev. Genet. 15 2014 293 306 24662220
48 Neumeyer S. Hemani G. Zeggini E. Strengthening causal inference for complex disease using molecular quantitative trait loci Trends Mol. Med. 26 2020 232 241 31718940
49 Porcu E. Rüeger S. Lepik K. eQTLGen ConsortiumBIOS ConsortiumSantoni F.A. Reymond A. Kutalik Z. Mendelian randomization integrating GWAS and eQTL data reveals genetic determinants of complex and clinical traits Nat. Commun. 10 2019 3300 31341166
50 Hou L. Xiong X. Park Y. Boix C. James B. Sun N. He L. Patel A. Zhang Z. Molinie B. Multitissue H3K27ac profiling of GTEx samples links epigenomic variation to disease Nat. Genet. 55 2023 1665 1676 37770633
51 Oliva M. Demanelis K. Lu Y. Chernoff M. Jasmine F. Ahsan H. Kibriya M.G. Chen L.S. Pierce B.L. DNA methylation QTL mapping across diverse human tissues provides molecular links between genetic variation and complex traits Nat. Genet. 55 2023 112 122 36510025
52 Zheng Z. Huang D. Wang J. Zhao K. Zhou Y. Guo Z. Zhai S. Xu H. Cui H. Yao H. QTLbase: an integrative resource for quantitative trait loci across multiple human molecular phenotypes Nucleic Acids Res. 48 2020 D983 D991 31598699
53 Zhu Z. Zhang F. Hu H. Bakshi A. Robinson M.R. Powell J.E. Montgomery G.W. Goddard M.E. Wray N.R. Visscher P.M. Yang J. Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets Nat. Genet. 48 2016 481 487 27019110
54 Giambartolomei C. Vukcevic D. Schadt E.E. Franke L. Hingorani A.D. Wallace C. Plagnol V. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics PLoS Genet. 10 2014 e1004383
55 Shrine N. Guyatt A.L. Erzurumluoglu A.M. Jackson V.E. Hobbs B.D. Melbourne C.A. Batini C. Fawcett K.A. Song K. Sakornsakolpat P. New genetic signals for lung function highlight pathways and chronic obstructive pulmonary disease associations across multiple ancestries Nat. Genet. 51 2019 481 493 30804560
56 Wong A.H. Tran T. CD151 in Respiratory Diseases Front. Cell Dev. Biol. 8 2020 64 32117989
57 Planès R. Pinilla M. Santoni K. Hessel A. Passemar C. Lay K. Paillette P. Valadão A.-L.C. Robinson K.S. Bastard P. Human NLRP1 is a sensor of pathogenic coronavirus 3CL proteases in lung epithelial cells Mol. Cell. 82 2022 2385 2400.e9 35594856
58 Pang S.W. Lahiri C. Poh C.L. Tan K.O. PNMA family: Protein interaction network and cell signalling pathways implicated in cancer and apoptosis Cell. Signal. 45 2018 54 62 29378289
59 Schüller M. Jenne D. Voltz R. The human PNMA family: novel neuronal proteins implicated in paraneoplastic neurological disease J. Neuroimmunol. 169 2005 172 176 16214224
60 Epi4K ConsortiumEpilepsy Phenome/Genome ProjectAllen A.S. Berkovic S.F. Cossette P. Delanty N. Dlugos D. Eichler E.E. Epstein M.P. Glauser T. De novo mutations in epileptic encephalopathies Nature 501 2013 217 221 23934111
61 Vatta M. Tennison M.B. Aylsworth A.S. Turcott C.M. Guerra M.P. Eng C.M. Yang Y. A novel STXBP1 mutation causes focal seizures with neonatal onset J. Child Neurol. 27 2012 811 814 22596016
62 Di Meglio C. Lesca G. Villeneuve N. Lacoste C. Abidi A. Cacciagli P. Altuzarra C. Roubertie A. Afenjar A. Renaldo-Robin F. Epileptic patients with de novo STXBP1 mutations: Key clinical features based on 24 cases Epilepsia 56 2015 1931 1940 26514728
63 Ernst J. Kellis M. ChromHMM: automating chromatin-state discovery and characterization Nat. Methods 9 2012 215 216 22373907
64 ENCODE Project Consortium An integrated encyclopedia of DNA elements in the human genome Nature 489 2012 57 74 22955616
65 Sun T. Xu Y. Xiang Y. Ou J. Soderblom E.J. Diao Y. Crosstalk between RNA m6A and DNA methylation regulates transposable element chromatin activation and cell fate in human pluripotent stem cells Nat. Genet. 55 2023 1324 1335 37474847
66 Li R. Zhao H. Huang X. Zhang J. Bai R. Zhuang L. Wen S. Wu S. Zhou Q. Li M. Super-enhancer RNA m6A promotes local chromatin accessibility and oncogene transcription in pancreatic ductal adenocarcinoma Nat. Genet. 55 2023 2224 2234 37957340
67 Zhao W. Zhang S. Zhu Y. Xi X. Bao P. Ma Z. Kapral T.H. Chen S. Zagrovic B. Yang Y.T. Lu Z.J. POSTAR3: an updated platform for exploring post-transcriptional regulation coordinated by RNA-binding proteins Nucleic Acids Res. 50 2022 D287 D294 34403477
68 Szklarczyk D. Gable A.L. Lyon D. Junge A. Wyder S. Huerta-Cepas J. Simonovic M. Doncheva N.T. Morris J.H. Bork P. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets Nucleic Acids Res. 47 2019 D607 D613 30476243
69 Shi H. Wei J. He C. Where, when, and how: context-dependent functions of RNA methylation writers, readers, and erasers Mol. Cell. 74 2019 640 650 31100245
70 Mattei A.L. Bailly N. Meissner A. DNA methylation: a historical perspective Trends Genet. 38 2022 676 707 35504755
71 Tan W.L.W. Anene-Nzelu C.G. Wong E. Lee C.J.M. Tan H.S. Tang S.J. Perrin A. Wu K.X. Zheng W. Ashburn R.J. Epigenomes of human hearts reveal new genetic variants relevant for cardiac disease and phenotype Circ. Res. 127 2020 761 777 32529949
72 Ng B. White C.C. Klein H.-U. Sieberts S.K. McCabe C. Patrick E. Xu J. Yu L. Gaiteri C. Bennett D.A. An xQTL map integrates the genetic architecture of the human brain’s transcriptome and epigenome Nat. Neurosci. 20 2017 1418 1426 28869584
73 Giambartolomei C. Zhenli Liu J. Zhang W. Hauberg M. Shi H. Boocock J. Pickrell J. Jaffe A.E. CommonMind ConsortiumPasaniuc B. Roussos P. A Bayesian framework for multiple trait colocalization from summary association statistics Bioinformatics 34 2018 2538 2545 29579179
74 Wang X. Tucker N.R. Rizki G. Mills R. Krijger P.H. de Wit E. Subramanian V. Bartell E. Nguyen X.-X. Ye J. Discovery and validation of sub-threshold genome-wide association study loci using epigenomic signatures eLife 5 2016 e10557 10.7554/eLife.10557
75 Manolio T.A. Collins F.S. Cox N.J. Goldstein D.B. Hindorff L.A. Hunter D.J. McCarthy M.I. Ramos E.M. Cardon L.R. Chakravarti A. Finding the missing heritability of complex diseases Nature 461 2009 747 753 19812666
76 Dvorakova M. Kubik-Zahorodna A. Straiker A. Sedlacek R. Hajkova A. Mackie K. Blahos J. SGIP1 is involved in regulation of emotionality, mood, and nociception and modulates in vivo signalling of cannabinoid CB1 receptors Br. J. Pharmacol. 178 2021 1588 1604 33491188
77 Hájková A. Techlovská Š. Dvořáková M. Chambers J.N. Kumpošt J. Hubálková P. Prezeau L. Blahos J. SGIP1 alters internalization and modulates signaling of activated cannabinoid receptor 1 in a biased manner Neuropharmacology 107 2016 201 214 26970018
78 Doherty A. Smith-Byrne K. Ferreira T. Holmes M.V. Holmes C. Pulit S.L. Lindgren C.M. GWAS identifies 14 loci for device-measured physical activity and sleep duration Nat. Commun. 9 2018 5257 30531941
79 Zhou Z.-D. Sathiyamoorthy S. Tan E.-K. LINGO-1 and neurodegeneration: pathophysiologic clues for essential tremor Tremor Other Hyperkinet. Mov. 2 2012 tre-02-51-249-1 10.7916/D8PZ57JV
80 de Wit J. Hong W. Luo L. Ghosh A. Role of leucine-rich repeat proteins in the development and function of neural circuits Annu. Rev. Cell Dev. Biol. 27 2011 697 729 21740233
81 Zhang M. Zhai Y. An X. Li Q. Zhang D. Zhou Y. Zhang S. Dai X. Li Z. DNA methylation regulates RNA m6A modification through transcription factor SP1 during the development of porcine somatic cell nuclear transfer embryos Cell Prolif. 57 2024 e13581
82 Bulik-Sullivan B.K. Loh P.-R. Finucane H.K. Ripke S. Yang J. Schizophrenia Working Group of the Psychiatric Genomics ConsortiumPatterson N. Daly M.J. Price A.L. Neale B.M. LD Score regression distinguishes confounding from polygenicity in genome-wide association studies Nat. Genet. 47 2015 291 295 25642630
83 Finucane H.K. Bulik-Sullivan B. Gusev A. Trynka G. Reshef Y. Loh P.-R. Anttila V. Xu H. Zang C. Farh K. Partitioning heritability by functional annotation using genome-wide association summary statistics Nat. Genet. 47 2015 1228 1235 26414678
84 Shrine N. Izquierdo A.G. Chen J. Packer R. Hall R.J. Guyatt A.L. Batini C. Thompson R.J. Pavuluri C. Malik V. Multi-ancestry genome-wide association analyses improve resolution of genes and pathways influencing lung function and chronic obstructive pulmonary disease risk Nat. Genet. 55 2023 410 422 36914875
85 Eckner R. Ewen M.E. Newsome D. Gerdes M. DeCaprio J.A. Lawrence J.B. Livingston D.M. Molecular cloning and functional analysis of the adenovirus E1A-associated 300-kD protein (p300) reveals a protein with properties of a transcriptional adaptor Genes Dev. 8 1994 869 884 7523245
86 Ogryzko V.V. Schiltz R.L. Russanova V. Howard B.H. Nakatani Y. The transcriptional coactivators p300 and CBP are histone acetyltransferases Cell 87 1996 953 959 8945521
87 Ye N. Zhang N. Zhang Y. Qian H. Wu B. Sun Y. Cul4a as a New Interaction Protein of PARP1 Inhibits Oxidative Stress-Induced H9c2 Cell Apoptosis Oxid. Med. Cell. Longev. 2019 2019 4273261
88 Caputo L. Witzel H.R. Kolovos P. Cheedipudi S. Looso M. Mylona A. van IJcken W.F.J. Laugwitz K.-L. Evans S.M. Braun T. The Isl1/Ldb1 Complex Orchestrates Genome-wide Chromatin Organization to Instruct Differentiation of Multipotent Cardiac Progenitors Cell Stem Cell 17 2015 287 299 26321200
89 Peng S. Sun M. Sun X. Wang X. Jin T. Wang H. Han C. Meng T. Li C. Plasma levels of TAM receptors and ligands in severe preeclampsia Pregnancy Hypertens. 13 2018 116 120 30177037
90 McShane L. Tabas I. Lemke G. Kurowska-Stolarska M. Maffia P. TAM receptors in cardiovascular disease Cardiovasc. Res. 115 2019 1286 1295 30980657
91 Connally N.J. Nazeen S. Lee D. Shi H. Stamatoyannopoulos J. Chun S. Cotsapas C. Cassa C.A. Sunyaev S.R. The missing link between genetic association and regulatory function eLife 11 2022 e74970 10.7554/eLife.74970
92 Moffat J. Grueneberg D.A. Yang X. Kim S.Y. Kloepfer A.M. Hinkle G. Piqani B. Eisenhaure T.M. Luo B. Grenier J.K. A lentiviral RNAi library for human and mouse genes applied to an arrayed viral high-content screen Cell 124 2006 1283 1298 16564017
93 Hemani G. Zheng J. Elsworth B. Wade K.H. Haberland V. Baird D. Laurin C. Burgess S. Bowden J. Langdon R. The MR-Base platform supports systematic causal inference across the human phenome eLife 7 2018 e34408 10.7554/eLife.34408
94 Quinlan A.R. Hall I.M. BEDTools: a flexible suite of utilities for comparing genomic features Bioinformatics 26 2010 841 842 20110278
95 GTEx Consortium The GTEx Consortium atlas of genetic regulatory effects across human tissues Science 369 2020 1318 1330 32913098
96 Burgess S. Butterworth A. Thompson S.G. Mendelian randomization analysis with multiple genetic variants using summarized data Genet. Epidemiol. 37 2013 658 665 24114802
97 Roadmap Epigenomics ConsortiumKundaje A. Meuleman W. Ernst J. Bilenky M. Yen A. Heravi-Moussavi A. Kheradpour P. Zhang Z. Wang J. Integrative analysis of 111 reference human epigenomes Nature 518 2015 317 330 25693563
98 Elsworth B. Lyon M. Alexander T. Liu Y. Matthews P. Hallett J. Bates P. Palmer T. Haberland V. Smith G.D. The MRC IEU OpenGWAS Data Infrastructure Preprint at bioRxiv 2020 10.1101/2020.08.10.244293
99 1000 Genomes Project ConsortiumAuton A. Brooks L.D. Durbin R.M. Garrison E.P. Kang H.M. Korbel J.O. Marchini J.L. McCarthy S. McVean G.A. Abecasis G.R. A global reference for human genetic variation Nature 526 2015 68 74 26432245
