Pleiotropic influence of DNA methylation QTLs on physiological and ageing traits

Figures & data

Figure 1. Graphical summary. Genetic variation in the trans-meQTL hotspot, meQTL.5a (black and pale-brown represent the genotype variants of meQTL.5a), influences CpG methylation at multiple developmental genes. DD in meQTL.5a is associated with higher methylation in the downstream targets. As these CpGs also gain methylation with aging, the DD genotype has a more “aged methylome.” The variation in methylation could contribute to genotype-dependent variation in lifespan. (Emojis downloaded from https://dryicons.com).

Figure 2. Overview of methylation QTL (meQTLs) in the liver. (a) plot of the genome-wide peak LOD score for the 27,966 CpGs, and distance between the CpG and the maximum LOD marker. Yellow: D allele (DBA/2J genotype) has positive additive effect; blue: B allele (C57BL/6J genotype) has positive additive effect. (b) meQTL location (x-axis; from chromosomes 1–19) and counts of meQTLs with LOD ≥ 3.5. (c) the genome graph plots location of the genome-wide peak QTL marker (x-axis), and location of the linked CpG (y-axis). Shows only the 3921 CpGs that map at LOD ≥ 3.5. Relative enrichment or depletion in (d) genomic location, and (e) predicted chromatin states for CpGs that map as cis-meQTLs (black), and as trans-meQTLs (grey). Asterisks denote hypergeometric enrichment p < 0.001. (f) enrichment or depletion in differentially methylated CpGs among the cis- and trans-modulated CpGs. (Asterisks denote hypergeometric enrichment p < 0.001; hash denotes p = 0.003) (g) Portion of the peak interval in the chromosome 5 meQTL hotspot: meQTL.5a (from UCSC Genome browser GRCm38/mm10).

Table 1. Liver meQTL hotspots.

Name	Broad location (peak)	Liver meQTLs	Hotspot type	Module QTL	Liver eQTLs	me/eQTL overlap
meQTL.2a	Chr2:102–112 (107 Mb)	72 (41 cis)	cis methyl		80 (13 cis)	Them7 (cg10774906);Mpped2 (cg07667286, cg00811894)
meQTL.2b	Chr2:141–151 (144–149 Mb)	164 (15 cis)	trans methyl	Green (2092)	27 (13 cis)	Rem1 (cg23754359, cg25361894) Rin2 (cg12687767)
meQTL.4a	Chr4:114–124 (119 Mb)	39 cis	cis methyl		43 (25 cis)	Faah (cg08815464, cg15610892, cg19641802, cg06184921)
meQTL.5a	Chr5:110–120 (115 Mb)	535 (33 cis)	trans methyl	Blue (5067); Black (1087)	376 (36 cis)	Tenm3 (trans; cg24399106), Clcn3 (trans; cg16842643)
meQTL.7a	Chr7:132–142 (138 Mb)	59 (9 cis)	trans methyl		197 (7 cis)	Mgmt (cg00046614, cg11711038, cg23272565)
meQTL.14a	Chr14:17–27 (23 Mb)	51 (32 cis)	cis methyl		5 cis	-
meQTL.14b	Chr14: 41–51 (46–49 Mb)	184 (19 cis)	trans methyl	Greenyellow (998)	97 (31 cis)	Ddhd1 (cg17695612; cg10924987);Prmt5 (cg24106188)
distal Eaa19^1	Chr19:42–52 (48 Mb)	103 (100 cis)	cis methyl	Royalblue (62); Lightgreen (1761)	274 (43 cis)	Crtac1; Ldb1; Psd; Col17a1

Note: ¹Epigenetic age acceleration.

Figure 3. Genetics of co-methylation CpG networks. QTL maps for four module eigengenes (MEs) are shown. Mapping was done using a linear mix model. The horizontal dashed red line marks a relatively lenient threshold of –log₁₀p = 3.5. The Blue, black and Lightgreen modules share suggestive overlapping QTLs in meQtl.5a. Eaa19 has a strong cis-regulatory effect on the Royalblue module. The chromosome 2 peak for the black module is proximal to meQTL.2b in Table 1.

Figure 4. Enriched functional pathways genomic regions enrichment analysis of (a) CpGs that are trans-modulated by meQTL.5a, and (b) CpGs that are members of the Blue module. Gene ontology and pathway enrichment among (c) transcripts, and (d) proteins that map to meQTL.5a in liver.

Table 2. Candidate genes in meQtl.5a.

Symbol	Brief description of annotated function	Chr	Mb	Missense variants (SIFT score)^1
Tchp	Apoptotic process; negative regulator of cell growth; cytoskeleton	5	114.7	rs29566070 (1)
Git2	G-protein coupled receptor protein signalling pathway; related to brain development	5	114.7	rs32133813 (0.86–1)
1500011B03Rik		5	114.8
Oasl2	Purine nucleotide biosynthetic process; immune signalling	5	114.9	rs29822904 (0.49–0.97); rs32142001 (0)
Gm13822		5	114.9
Hnf1a	Liver-enriched transcription factor; development & growth; MODY3	5	114.9	rs33234601 (0.14)
Acads	Mitochondrial flavoprotein; acyl-CoA dehydrogenase family; fatty acid beta-oxidation pathway	5	115.1
Coq5	Mitochondrial co-enzyme; methylation and ubiquinone biosynthetic process; CoQ10 biosynthesis pathway	5	115.3
Srsf9	Pre-mRNA splicing factor; mRNA export	5	115.3	rs33739429 (0)
Gatc	Glutaminyl-tRNA synthase; mitochondrial	5	115.3	rs33338640 (0.21)
Triap1	p53 binding; DNA damage response; negative regulation of apoptosis; phospholipid transport	5	115.3	rs33338640 (0.21)
Sirt4	Sirtuin member; mitochondrial; metabolic processes; mono-ADP-ribosyltransferase	5	115.5	rs46787798 (0.58); rs6400038 (0.32)
Pxn	Cytoskeletal; angiogenesis; transforming growth factor beta receptor signalling	5	115.5	rs50194001 (0.39); rs52040466 (0.87); rs46615100 (0.34); rs33590215 (1); rs50879465 (1); rs47873388 (0.03); rs33728337 (0.28); rs33892383 (1)
Rplp0	Ribosomal protein	5	115.6	rs52016292 (1)
Rab35	GTpase activity; neuron projection; mitochondrial	5	115.6	rs48405889 (0)
Ccdc64	Small GTPase binding; dynactin binding; golgi to secretory granule transport; neuron projection	5	115.6	rs47577059 (1)
Cit	Serine/threonine-protein kinase; cell division; central nervous system development	5	115.8	rs48791426 (1); rs47954950 (0.01); rs48893178 (0.05); rs48063202 (0)
Prkab1	AMP-activated protein kinase; energy sensing; metabolism	5	116.0	rs46168068 (0)
Hspb8	Heat shock protein; unfolded protein response; Charcot-Marie-Tooth disease	5	116.4

Note ¹The SIFT scores provided within parenthesis predicts the effect of missense mutation on protein function; lower scores are more likely to be deleterious. SIFT scores for the missense variants were obtained from the Ensembl Browser.

Figure 5. Interaction networks that connect the trans-meQTLs to the candidate genes in meQTL.5a. (a) the trans-modulated CpGs and candidate genes in meQTL.5a form a highly connected and functionally enriched protein–protein interaction network. HNF1A (yellow triangle) is the most connected candidate in the central sub-network, and member of the enriched MODY (maturity onset diabetes of the young) pathway (red nodes). (b) CREBBP is the hub gene in the meQTL-based network, and a CpG in its exon is trans-modulated by meQtl.5a (top). Its mRNA (bottom of mirrored Manhattan plot) has a weak peak in meQTL.5a. Crebbp is located on chromosome 16 (red triangle). (c) the CpG in the 5’UTR of Hnf1a is cis-modulated, and its expression has a weak cis-eQTL. Strong positive correlation between the CpGs (d), and between the mRNAs (e) of Hnf1a and Crebbp. Weak inverse correlation between (f) the Hnf1a mRNA and Crebbp methylation, and (g) between the mRNA and methylation of Hnf1a.

Figure 6. Primary protein–protein interaction partners of HNF1A and their trans-modulation (a) the network shows the top 10 interaction partners of HNF1A based on protein–protein interactions. CpGs in CREBBP, FOXA2 (HNF3B), GATA4, HNF4A, ONECUT1, and PCBD2 are trans-modulated by the meQTL.5a locus, making Hnf1a a prime candidate. (b) at the transcriptomic level, expression of these genes in the liver are also highly intercorrelated. Only correlations |r| > 0.45 are shown (line thickness conveys strength of correlation; blue: negative; red: positive). Mirrored meQTL (top), and eQTL (bottom) for two example gene: (c) Foxa2, and (d) Gata4. Red triangles mark the location of genes.

Figure 7. Joint modulation of CpGs by genetic and non-genetic variables (a) CpGs with trans-meQTLs in meQTL.5a are enriched in differentially methylated CpGs (DMCs) associated with aging, diet, and body weight. Asterisks denote hypergeometric enrichment p ≤ 0.001; hash denotes p = 0.006. (b) trans-CpGs that are increased in methylation by the D allele in meQTL.5a gain methylation with age (positive regression estimate on the y-axis) while those with higher methylation for the B allele lose methylation with age. Dashed red horizontal line indicate Bonferroni p ≤ 0.05 for DMC. (c) CpG in the Jarid2 5’UTR gains methylation with age and has higher methylation in BXDs with the DD genotype in meQTL.5a. (d) QTL plots for the Jarid2 CpG (top Manhattan plot), and mRNA (bottom). Red triangle marks the location of Jarid2. (e) allele dependent increase in methylation at the meQTL.5a trans-CpGs due to HFD. (f) allele dependent association with body weight for the trans-CpGs. (g) Overall mean methylation of the trans-modulated CpG is higher for BXDs with DD genotype in meQTL.5a for both control diet (CD; 0.35 ± 0.02 for BB; 0.38 ± 0.03 for DD; pair-wise p < 0.0001) and high fat diet (HFD; 0.36 ± 0.03 for BB; 0.38 ± 0.03 for DD pair-wise p < 0.005). (h) allele dependent increase in mean methylation for the CpGs with trans-meQTL in meQTL.5a. (i) body weight is slightly higher for BXDs with BB in meQTL.5a for both CD (27 ± 6 g for BB; 26 ± 5 g for DD; pair-wise p < 0.1) and HFD (47 ± 13 g for BB; 42 ± 13 g for DD; pair-wise p < 0.02).

Table 3. Multivariable regressions for mean methylation and weight.

Outcome	Predictors	Estimate	Std Error	t Ratio	p
Mean methylation of trans-CpGs^1	rs29733222[B]	−0.011	0.002	−7.27	<.0001
	Age (days)	6.4E–05	8.4E–06	7.64	<.0001
	Diet[CD]	−0.008	0.0021	−3.84	0.0002
	Body weight	−0.0002	0.0002	−1.45	0.15
Body weight^2	rs29733222[B]	1.517	0.5316	2.85	0.005
	Age (days)	0.006	0.0030	2.17	0.03
	Diet[CD]	−8.94	0.5415	−16.51	<.0001

Note ¹(Mean methylation for meQTL.5a trans-CpGs) ~ genotype + age + diet + weight, where diet is control diet (CD) or high fat diet (HFD), and genotype is BB (72 on CD, 50 on HFD) or DD (105 on CD, 63 on HFD) for marker rs29733222 in meQTL.5a. ²Weight ~ genotype + age + diet.

Figure 8. Pleiotropic influence on meQTL.5a on body weight at young age and lifespan (a) body weight at 6 months (mos) from a separate cohort of BXD mice show higher mean weight for strains with BB genotype in meQTL.5a for control diet (CD; 26 ± 6 g for BB; 25 ± 5 g for DD; pair-wise p < 0.004) and high fat (HFD; 34 ± 9 g for BB; 31 ± 8 g for DD; pair-wise p < 0.0001). Samples numbers: 383 BB and 503 DD for CD; 392 BB and 457 DD for HFD. (b) Kaplan-Meir survival plots by genotype at meQTL.5a for CD mice (399 BB and 510 DD). Median lifespan in days (MedianLS) for the genotypes shown (log-rank p = 0.008). (c) similar survival plot for HFD shows no significant difference between genotypes (402 BB and 462 DD). (d) Kaplan-Meir survival after adjusting lifespan for 6 mos weight. Adjusted median lifespan in days shown for each genotype-by-diet below the graph. Within each diet, the BB genotype has longer lifespan compared to DD (based on pairwise comparison, log-rank p < 0.0001 for CD, and p = 0.03 for HFD). (e) model depicting horizontal pleiotropic influence on CpG methylation and weight, and vertical pleiotropic influence on lifespan mediated by CpG, which are also under the influence of ageing and diet.

Table 4. Multivariable regressions for body weight at 6 months of age and strain lifespan.

Outcome	Term	Estimate	Std Error	t Ratio	p
Weight_6M^1	rs29733222[B]	0.98	0.17	5.78	<.0001
	Diet[CD]	−3.38	0.17	−20.11	<.0001
Lifespan^2	rs29733222[B]	14	4.41	3.08	0.002
	Diet[CD]	26	4.82	5.43	<.0001
	Weight_6M	−4.0	0.62	−6.34	<.0001

Note ¹(Weight at 6 months) ~ genotype + diet. ²Lifespan ~ genotype + diet + weight at 6 months, where diet is control diet (CD) or high fat diet (HFD), and genotype is BB (399 on CD, 402 on HFD) or DD (510 on CD, 462 on HFD) for marker rs29733222 is meQTL.5a.

Table 5. Phenotypes that map to meQtl.5a.

GN number^1	Phenotype	Category	PMID	-log10(p)	QTL peak location
17439	Rickettsia tsutsugamushi susceptibility of both sexes at 6–12 weeks-of-age	Immune	6774020	4.42	Peak at 104 @rs32034514
17266	Brain activity and coherence of electrical field oscillations at 160 Hz in L2/3 of the primary whisker motor cortex	Nervous system	24686563	4.03	Peak at ~ 107 Mb @ rs29681689
15092	Total fat content measured by Fourier Transform Infrared Spectroscopy in liver at 140 days, males, fed high fat diet fed from 4 weeks on	Metabolism	23758785	3.94	Peak from 110–114 Mb
15091	Saturated fat content measured by Fourier Transform Infrared Spectroscopy in liver at 140 days, males, fed high fat diet feeding from 4 weeks on	Metabolism	23758785	3.62	Peak from 110–114 Mb
15094	Ratio of lipid to protein content in liver at 140 days, males, fed high fat diet from 4 weeks on	Metabolism	23758785	3.44	Peak from 110–114 Mb
14786	Proliferation of BrdU-labelled cells in subgranular zone, 1 h BrdU injection, unadjusted data	Nervous system	24640950	3.43	Peak at 117 Mb @rs29728022
17265	Brain activity and coherence of electrical field oscillations at 159 Hz in L2/3 of the primary whisker motor cortex	Nervous system	24686563	3.31	Peak at ~ 107 Mb @ rs29681689
16783	Ratio of total branched-chain amino acid/total amino acid	Metabolites	22939713; 30709776	3.23	Peak at 110 Mb @rs49420585
17245	Brain activity and coherence of electrical field oscillations, coherence at 139 Hz between local field potentials (LFP) at two sites (0.3 mm apart) in L2/3 of the primary whisker motor cortex in awake 6	Nervous system	24686563	3.19	Peak at ~ 107 Mb @ rs29681689
16774	Ratio of total branched-chain amino acid/Alanine_CD	Metabolites	22939713; 30709776	3.1	Peak at 110 Mb @rs49420585

Note ¹Search carried out using Hnf1a as the search key in https://systems-genetics.org/; the traits ID can be used to retrieve the BXD strain level data from www.genenetwork.org.

Data availability statement

The data used in this study are available from NCBI Gene Expression Omnibus (GEO accession ID GSE199979). Genotype data are provided as Supplemental Data.