Jacob Bergstedt
email: jacob.bergstedt@ki.se

Sadoune Ait

Kaci Azzou

Kristin Tsuo

Anthony Jaquaniello

Alejandra Urrutia

Maxime Rotival

David T S Lin

Julia L Macisaac

Michael S Kobor

Matthew L Albert

Darragh Duffy

Etienne Patin
email: epatin@pasteur.fr

Lluís Quintana-Murci

Factors Driving DNA Methylation Variation in Human Blood

Keywords: epigenetics, DNA methylation, humans, blood, immune cells, aging, cytomegalovirus, smoking, sex, meQTLs, EWAS, gene-by-environment interactions, cellular deconvolution

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Epigenetic research has facilitated our understanding of the existing links between environmental risk factors, aging, genetic variation and human disease [START_REF] Cavalli | Advances in epigenetics link genetics to the environment and disease[END_REF][START_REF] Michalak | The roles of DNA, RNA and histone methylation in ageing and cancer[END_REF]. Epigenome-wide association studies (EWAS) have shown that DNA methylation (i.e., 5-methylcytosine, 5mC), the most studied epigenetic mark in human populations, is associated with a wide range of environmental exposures along the life course, such as chemicals, air pollution and nutrition [START_REF] Martin | Environmental Influences on the Epigenome: Exposure-Associated DNA Methylation in Human Populations[END_REF], as well as past socioeconomic status [START_REF] Bush | The biological embedding of early-life socioeconomic status and family adversity in children's genome-wide DNA methylation[END_REF][START_REF] Karlsson Linner | An epigenome-wide association study meta-analysis of educational attainment[END_REF][START_REF] Lam | Factors underlying variable DNA methylation in a human community cohort[END_REF][START_REF] Stringhini | Life-course socioeconomic status and DNA methylation of genes regulating inflammation[END_REF]. Changes in DNA methylation have also been associated with non-communicable diseases, such as Parkinson's and Alzheimer's diseases, multiple sclerosis, systemic lupus erythematosus, type 2 diabetes and cardiovascular disease [START_REF] Hwang | The emerging field of epigenetics in neurodegeneration and neuroprotection[END_REF][START_REF] Ling | Epigenetics in Human Obesity and Type 2 Diabetes[END_REF][START_REF] Mazzone | The emerging role of epigenetics in human autoimmune disorders[END_REF][START_REF] Van Der Harst | Translational Perspective on Epigenetics in Cardiovascular Disease[END_REF].

These studies collectively suggest that DNA methylation marks could have tremendous value as a gauge of the exposome and as clinical biomarkers [START_REF] Berdasco | Clinical epigenetics: seizing opportunities for translation[END_REF][START_REF] Wild | Complementing the genome with an "exposome": the outstanding challenge of environmental exposure measurement in molecular epidemiology[END_REF].

However, interpretation of EWAS remains limited. First, because the epigenome of a cell reflects its identity [START_REF] Farlik | DNA Methylation Dynamics of Human Hematopoietic Stem Cell Differentiation[END_REF][START_REF] Roadmap Epigenomics | Integrative analysis of 111 reference human epigenomes[END_REF], a risk factor or a disease that alters cellular composition also alters 5mC levels measured in the tissue [START_REF] Liu | Epigenome-wide association data implicate DNA methylation as an intermediary of genetic risk in rheumatoid arthritis[END_REF]. It is thus necessary to determine if an exposure affects cellular composition or DNA methylation states of cell types, in order to better understand the link between such an exposure, DNA methylation and disease [START_REF] Lappalainen | Associating cellular epigenetic models with human phenotypes[END_REF]. Previous studies have accounted for cellular heterogeneity in blood by using cell sorting experiments, or cellular proportions estimated from 5mC profiles through deconvolution techniques [START_REF] Houseman | DNA methylation arrays as surrogate measures of cell mixture distribution[END_REF][START_REF] Teschendorff | A comparison of referencebased algorithms for correcting cell-type heterogeneity in Epigenome-Wide Association Studies[END_REF], but these approaches focus on a subset of frequent cell types that capture only a part of blood cellular composition. Second, genetic variation and DNA methylation are inextricably linked, as attested by the numerous DNA methylation quantitative trait loci (meQTLs) detected so far [START_REF] Villicana | Genetic impacts on DNA methylation: research findings and future perspectives[END_REF], and genetic variants could confound associations between risk factors, 5mC levels and diseases or traits. Finally, environmental risk factors with a yet-unknown effect on DNA methylation, such as common infections, could also confound associations between other risk factors, DNA methylation and human phenotypes. Thus, a detailed study of the factors that impact DNA methylation at the population level, and the extent to which their effects are mediated by changes in cellular composition, is required to understand the role of epigenetic variation in health and disease.

To address this gap, we generated whole blood-derived DNA methylation profiles at >850,000

CpG sites for 958 healthy adults of the Milieu Intérieur cohort. We leveraged the deep characterization of the cohort, including high-resolution immunophenotyping by flow cytometry 4 [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF][START_REF] Thomas | The Milieu Interieur study -an integrative approach for study of human immunological variance[END_REF], to determine whether and how cellular composition, genetic variation, intrinsic factors (i.e., age and sex) and 139 diverse health-related variables and environmental exposures affect the blood DNA methylome. We performed EWAS adjusted or not for the measured proportions of 16 immune cell subsets, to robustly delineate effects on DNA methylation that are direct, i.e., acting through changes within cells, from those that are mediated, i.e., acting through subtle changes in cellular composition [START_REF] Houseman | Cell-composition effects in the analysis of DNA methylation array data: a mathematical perspective[END_REF]. We find that the largest effects on DNA methylation are due to DNA sequence variation, whereas the most widespread differences among individuals are the result of blood cellular heterogeneity. We also identify latent cytomegalovirus (CMV) infection as a major driver of epigenetic variation and observe an increased dispersion of DNA methylation with aging, suggesting a decrease in the fidelity of the epigenetic maintenance machinery. Finally, we show that a large part of the effects on DNA methylation of aging, smoking, CMV serostatus and chronic low-grade inflammation is due to fine-grained changes in blood cell composition, and characterize the DNA methylation signature of cell-types affected by these factors. This work generates new hypotheses about mechanisms underlying DNA methylation variation in the human population and highlights critical factors to be considered in medical epigenomics studies.

5

RESULTS

Population Variation in DNA Methylation Differs Markedly Across the Genome

To investigate the contributions of genetic and non-genetic factors to population variation in DNA methylation, we quantified 5mC levels at >850,000 CpG sites, with the Illumina Infinium MethylationEPIC array, in the 1,000 healthy donors of the Milieu Intérieur cohort [START_REF] Thomas | The Milieu Interieur study -an integrative approach for study of human immunological variance[END_REF]. The cohort includes individuals of Western European origin, equally stratified by sex (i.e., 500 women and 500 men) and age (i.e., 200 individuals from each decade between 20 and 70 years of age), who were surveyed for detailed demographic and health-related information, including dietary habits, upbringing, socioeconomic status, mental health, past and latent infections, and vaccination and medical histories (Table S1). To adjust for and quantify the blood cell-compositionmediated effect on DNA methylation of genetic factors, intrinsic factors (i.e., age and sex) and environmental exposures, we measured, in all donors, the proportions of 16 major and minor immune cell subsets by standardized flow cytometry, including neutrophils, basophils, eosinophils, monocytes, natural killer (NK) cells, dendritic cells, B cells, CD4 + and CD8 + T cells at four differentiation stages (naive, central memory, effector memory and terminally differentiated effector memory cells) and CD4 -CD8 -T cells [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. Furthermore, all donors were genotyped at 945,213 single-nucleotide polymorphisms (SNPs), yielding 5,699,237 accurate SNPs after imputation [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. After quality control filtering, high-quality measurements of 5mC levels were obtained at 644,517 CpG sites for 958 unrelated individuals (Figure S1; STAR Methods).

We first investigated population variation in DNA methylation across different wellcharacterized chromatin states, using naïve CD4 + T cells as a reference [START_REF] Roadmap Epigenomics | Integrative analysis of 111 reference human epigenomes[END_REF]. We found that CpG sites in transcription start sites (TSS) are typically unmethylated and exhibit the lowest population variance in 5mC levels (Figure 1A,B), suggesting that epigenetic constraints are the strongest in promoters, whereas actively transcribed gene bodies and heterochromatin are highly methylated and also show low population variance. 5mC levels in enhancers and Polycomb-repressed regions are the most variable (Figure 1A,B), suggesting that DNA methylation in these regions are preferentially affected by genetic, intrinsic or environmental factors, or cellular heterogeneity. These results indicate that 5mC measurements from our cohort reproduce the known properties of DNA methylation and show high levels of variation across the epigenome and among individuals.

Widespread Local Genetic Control of DNA Methylation in Whole Blood

Studies of DNA methylation quantitative trait loci (meQTLs) have revealed that DNA sequence variation affects 5mC levels at numerous nearby CpG sites [START_REF] Bonder | Disease variants alter transcription factor levels and methylation of their binding sites[END_REF][START_REF] Hannon | Leveraging DNA-Methylation Quantitative-Trait Loci to Characterize the Relationship between Methylomic Variation, Gene Expression, and Complex Traits[END_REF], a feature that can confound associations between environmental exposures and DNA methylation [START_REF] Lappalainen | Associating cellular epigenetic models with human phenotypes[END_REF]. To account for this, we estimated, for each CpG site, the effects of nearby DNA sequence variation on their 5mC levels (100-Kb window; STAR Methods). We adjusted models on genetic ancestry and blood cell-type proportions, as well as age, sex, smoking status and CMV serostatus, which we have previously shown to affect blood cell composition [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. We found a significant local meQTL for 29.2% of the 644,517 tested CpG sites (n = 188,129; two-stage adjusted P-value, P adj < 0.05; Figure S2; STAR methods).

We detected 1,978 CpG sites with a difference in the proportion of DNA methylation between homozygotes larger than 0.3, indicating that the effect of local meQTLs can be substantial.

We found that CpG sites with a local meQTL are enriched in enhancers (odds ratio [OR] 95% CI: [1.87, 1.95]; Figure 1C; STAR Methods), where the population variance of DNA methylation is the largest (Figure 1B). Conversely, CpG sites with a local meQTL are depleted in TSS and actively transcribed genes (OR CIs: [0.32,0.34] and [0.49, 0.52]), where DNA methylation variance is the lowest and sequence conservation the highest (Figure 1D). While confirming that local meQTLs are enriched in disease and trait associations by genome-wide association studies (GWAS) (Bonder et al., 2017), we found that the enrichment is strongest in enhancers and genic enhancers (enrichment > 1.8; P resampling < 1.0×10 -4 ; Figure 1E) and absent from regions of low sequence conservation, such as heterochromatin. These findings indicate that DNA sequence variants have widespread, strong effects on nearby DNA methylation levels, particularly in regulatory elements, and that local meQTLs are enriched in genetic variants that affect phenotypic variation and disease risk.

Structural Factors and Zinc Finger Proteins are Regulators of DNA Methylation

We investigated the long-range genetic control of DNA methylation, by estimating the effect of genome-wide variants on 5mC levels of a selection of CpG sites, to reduce the burden of multiple testing. We selected the 50,000 CpG sites with the highest residual variance after fitting a model including as predictors: (i) the most associated local meQTL variant, (ii) genetic ancestry, (iii) blood cell proportions and (iv) non-genetic factors affecting blood cell composition (STAR Methods). We found 2,394 independent long-range meQTLs, for 1,816 CpG sites (3.6%) and 1,761 independent variants (Table S2). The effects of long-range meQTLs are generally weaker than those of local meQTLs (Figure S2), yet we found 152 CpG sites with a difference in the proportion of methylation between homozygotes larger than 0.15. As for local meQTLs, remote-effect variants are also enriched in GWAS hits (enrichment = 1.78; P resampling < 6.3×10 -5 ). Interestingly, CpG sites under remote genetic control are enriched in TSS regions and regions associated with ZNF genes, in contrast with CpG sites under local genetic control (Figure 1C,F). Furthermore, we found that remote meQTL SNPs are strongly concentrated in ZNF genes (OR CI: [15.1,25.8]; Figure 1G).

These findings indicate that zinc-finger proteins (ZFPs) play a role in the long-range control of DNA methylation, in line with their role in the regulation of heterochromatin and ZNF gene expression [START_REF] Marchal | Emerging concept in DNA methylation: role of transcription factors in shaping DNA methylation patterns[END_REF][START_REF] O'geen | Genome-wide analysis of KAP1 binding suggests autoregulation of KRAB-ZNFs[END_REF][START_REF] Quenneville | The KRAB-ZFP/KAP1 system contributes to the early embryonic establishment of sitespecific DNA methylation patterns maintained during development[END_REF].

Genetic variants controlling the expression of nearby transcription factors (TF) have been found

to have widespread, long-range effects on the DNA methylome [START_REF] Bonder | Disease variants alter transcription factor levels and methylation of their binding sites[END_REF][START_REF] Hop | Genome-wide identification of genes regulating DNA methylation using genetic anchors for causal inference[END_REF]. Of the 3,643 genes with a nearby variant associated with a remote CpG site, 33% have its expression altered by the same variant in the eQTLGen database [START_REF] Võsa | Unraveling the polygenic architecture of complex traits using blood eQTL meta-analysis[END_REF]. In total, 200 protein-coding genes have local variants that are associated with 5mC levels at ≥ 10 remote CpG sites (Figure 1H). Variants that affect the largest number of remote CpG sites are located nearby well-known structural factors and TFs such as SENP7, BCLAF1, CTCF, NFKB1 and NFE2, and, consistently, CpG sites remotely associated with a TF local eQTL are strongly enriched in binding sites of the corresponding TF, or a TF related to it (Tables 1 andS3). For example, the rs10417143 variant alters ZNF257 mRNA levels (P adj = 2.8×10 -70 ; [START_REF] Võsa | Unraveling the polygenic architecture of complex traits using blood eQTL meta-analysis[END_REF]) and 5mC levels at 16

CpG sites, which are enriched in binding sites for the ZNF534 TF (Figure 1I). Likewise, rs12491955 is associated with increased SENP7 mRNA levels (P adj = 4.2×10 -302 ; [START_REF] Võsa | Unraveling the polygenic architecture of complex traits using blood eQTL meta-analysis[END_REF] and 5mC levels at 35 CpG sites, 30 of which are located in two clusters of ZNF genes on chromosome 19 [START_REF] Lemire | Long-range epigenetic regulation is conferred by genetic variation located at thousands of independent loci[END_REF]. Of these 35 CpG sites, 23 are located in binding sites for KAP1 (encoded by TRIM28, Figure 1J), a chromatin remodeler regulated by SENP7 [START_REF] Garvin | The deSUMOylase SENP7 promotes chromatin relaxation for homologous recombination DNA repair[END_REF]. Of note, 20 out of the 50 most enriched TFs are ZFPs (Table S3). These results collectively support the notion that transcriptional variation of TFs results in DNA methylome-wide changes due to differential occupancy of their binding sites and highlight the role of chromatin remodelers and ZFPs in the regulation of DNA methylation.

Aging Elicits DNA Hypermethylation Related to Polycomb Repressive Complexes and

Increased Epigenetic Dispersion

The link between DNA methylation and aging is well established [START_REF] Hannum | Genome-wide methylation profiles reveal quantitative views of human aging rates[END_REF][START_REF] Heyn | Distinct DNA methylomes of newborns and centenarians[END_REF][START_REF] Johansson | Continuous Aging of the Human DNA Methylome Throughout the Human Lifespan[END_REF][START_REF] Jones | DNA methylation and healthy human aging[END_REF][START_REF] Wang | Epigenetic influences on aging: a longitudinal genome-wide methylation study in old Swedish twins[END_REF]; however, given that blood cell composition is altered with age [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF], it remains unclear how age impacts DNA methylation in a heterogeneous tissue such as blood [START_REF] Jaffe | Accounting for cellular heterogeneity is critical in epigenomewide association studies[END_REF]. We thus investigated how the DNA methylome is shaped by the intertwined processes of cellular aging (i.e., direct, cell-composition-independent effects of age) and age-related changes in blood cellular composition (i.e., effects of age mediated by changes in cellular composition) (STAR Methods). At a false discovery rate (FDR) of 0.05, we found that age has a significant total (i.e., direct and mediated) effect on 5mC levels at 258,830 CpG sites (40.2% of CpG sites; Figure 2A and Table 2). When estimating direct effects by adjusting on measured immune cell proportions, a significant effect was observed at 144,114 CpG sites (22.4% of CpG sites; FDR < 0.05; Table 2), highlighting the widespread, cellcomposition-independent effect of age on 5mC levels in blood. Importantly, when we applied a commonly used deconvolution method to correct for cellular heterogeneity [START_REF] Koestler | Improving cell mixture deconvolution by identifying optimal DNA methylation libraries (IDOL)[END_REF], the number of age-associated CpG sites almost doubled (n = 277,209), suggesting that corrections based on estimated proportions of major cell subsets are incomplete.

We observed that 69% of CpG sites directly associated with age show a decrease in 5mC levels (Figure 2B,C). This pattern predominates in quiescent chromatin, actively transcribed genic regions and enhancers. In contrast, DNA hypermethylation was observed in 87% of age-associated CpGs within CpG islands (CGIs; Figure S3A,B). Consistently, CpG sites exhibiting increasing 5mC levels with age are predominantly found in Polycomb-repressed regions, bivalent TSSs and bivalent enhancers, which are CGI-rich regions (Figures 2B,C and S1K,L). Furthermore, these CpG sites are the most enriched in binding sites for RING1B, JARID2, RYBP, PCGF1, PCGF2 and SUZ12 TFs (OR > 10.0; Figure 2D and Table S3), which are all part of the Polycomb repressive complexes (PRC) 1 and 2. PRC1 and PRC2 mediate cellular senescence and modulate longevity in invertebrates [START_REF] Bracken | The Polycomb group proteins bind throughout the INK4A-ARF locus and are disassociated in senescent cells[END_REF][START_REF] Siebold | Polycomb Repressive Complex 2 and Trithorax modulate Drosophila longevity and stress resistance[END_REF]. Importantly, when we restricted the analysis to CpG sites outside of CpG islands, we found similar enrichments in Polycomb-repressed regions (OR 95% CI [10.4,14.0]) and PRC TF binding sites (RING1B OR 95% CI: [11.4, 14.1]; PCGF2 OR 95% CI [9.89,13.5]). Finally, genes with age-increasing methylation levels are strongly enriched in developmental genes (P adj = 2.2×10 -42 ; Table S4), which are regulated by PRCs [START_REF] Boyer | Polycomb complexes repress developmental regulators in murine embryonic stem cells[END_REF]. These results support a key regulatory role of Polycomb proteins in age-related hypermethylation [START_REF] Dozmorov | Polycomb repressive complex 2 epigenomic signature defines ageassociated hypermethylation and gene expression changes[END_REF].

Global hypomethylation of the genome and CGI-associated hypermethylation are both hallmarks of cancer [START_REF] Timp | Cancer as a dysregulated epigenome allowing cellular growth advantage at the expense of the host[END_REF]. We found that genes with a TSS that is increasingly methylated with age are significantly enriched in tumor suppressor genes (OR = 1.55,CI: [1.28,1.85]; Fisher's exact test P = 4.0×10 -6 ) [START_REF] Zhao | TSGene 2.0: an updated literature-based knowledgebase for tumor suppressor genes[END_REF]. For example, 5mC levels linearly increase by 0.2% per year of age nearby the TSS of BCL6B (P adj = 4.6×10 -66 ; Figure 2E), a tumor suppressor gene that is hypermethylated in cancer (Xu et al., 2012a). In addition, 5mC levels increase by 0.1% per year of age nearby the TSS of DNMT3A (P adj = 6.6×10 -17 ; Figure S3E), which encodes a DNA methyltransferase that plays a key role in tumorigenesis [START_REF] Jost | Epimutations mimic genomic mutations of DNMT3A in acute myeloid leukemia[END_REF]. Our results indicate that genomic hypomethylation and CGI-associated hypermethylation are lifelong progressive processes, possibly due to an altered maintenance of DNA methylation after multiple mitotic cell divisions [START_REF] Teschendorff | Age-dependent DNA methylation of genes that are suppressed in stem cells is a hallmark of cancer[END_REF][START_REF] Yang | Correlation of an epigenetic mitotic clock with cancer risk[END_REF][START_REF] Zhou | DNA methylation loss in late-replicating domains is linked to mitotic cell division[END_REF], and support an intricate link between aging and oncogenesis.

Finally, we determined whether the variance of 5mC levels among individuals depends on age, a phenomenon known as "epigenetic drift" (i.e., the divergence of the DNA methylome as a function of age due to stochastic changes) [START_REF] Fraga | Epigenetic differences arise during the lifetime of monozygotic twins[END_REF][START_REF] Jones | DNA methylation and healthy human aging[END_REF], by fitting models parameterizing the residual variance with a linear age term, and adjusting for cellular composition, aging, CMV serostatus and sex in the mean function (STAR methods). We observed a significant dispersion with age for 16.3% of all CpG sites. Strikingly, 90% of these CpGs show an increase in the variance of 5mC levels with age (Figure 2F), suggesting a decrease in the fidelity of epigenetic maintenance associated with aging. Examples of CpG sites with large, age-increasing dispersion are found in the TSS of the MAFA and CBLN1 genes (P adj = 2.08×10 -49 and 4.63×10 -45 ; Figure 2G,H).

Similar results were obtained when adjusting the variance function for cellular composition (STAR methods). In addition, we found that, out of 104,786 CpG sites with age-related dispersion, 63%

show no significant changes in 5mC levels with age (Figure S3F), which suggests that these results are not driven by relationships between the average and variance of 5mC levels. Collectively, these findings indicate that aging elicits numerous DNA methylation changes in a cell-compositionindependent manner, including global epigenome-wide demethylation, hypermethylation of PRCassociated regions, and increased variance, highlighting the occurrence of different mechanisms involved in epigenetic aging.

Immunosenescence-Related Changes in Cellular Composition Mediate DNA Methylation

Variation with Age

We detected a significant cell-composition-mediated effect of age at ~12% of CpG sites (n = 75,301; Table 2), indicating that a substantial fraction of age-associated changes in DNA methylation are due to age-related changes in immune cell proportions. However, mediated effects are typically weaker than direct effects (Figure 2C,K) and CpG sites with the strongest direct age effects show no mediated effect, suggesting that their changes in 5mC levels are shared across cell types (Figure 2A). In contrast to direct effects, mediated effects are most often associated with demethylation, regardless of their localization (Figures 2K andS3C, D). Yet, enhancers and TSSflanking regions (but not TSS themselves) are enriched in CpG sites with a significant cellcomposition-mediated, positive effect of age (Figure 2J), possibly because these regions tend to be active in a cell-type specific manner [START_REF] Roadmap Epigenomics | Integrative analysis of 111 reference human epigenomes[END_REF]. In addition, CpG sites 10 with a cell-composition-mediated increase in DNA methylation are enriched in TF binding sites for RUNX1 and RUNX3, two key regulators of hematopoiesis (Figure 2L and Table S3). Genes with CpG sites showing a mediated increase or decrease in DNA methylation with age are respectively enriched in genes involved in lymphoid or myeloid cell activation (P adj < 2.0×10 -14 ; Table S4), indicating that mediated age effects on DNA methylation are related to progressive differences in the composition of the lymphoid and myeloid compartments.

We then determined if age elicits 5mC changes specific to certain immune cell subsets, by deriving and verifying an interaction model capable of dealing with the compositional nature of immune cell proportions in blood, similarly to previous work [START_REF] Zheng | Identification of differentially methylated cell types in epigenome-wide association studies[END_REF] (STAR Methods). Because inference was unstable for rare cell subsets, we restricted the analysis to six major immune cell types (CD4 + and CD8 + T cells, CD4 -CD8 -T cells, B cells, NK cells, and monocytes), and estimated how DNA methylation changes with age within these cell types, compared to the most frequent cell type, i.e., neutrophils. We found that 17% of tested CpG sites (n = 106,899) show a cell-type-dependent association with age in CD8 + T cells (Figure 2I), 69% of which show decreasing 5mC levels (Figure S3G). These results support previous studies reporting that this T cell subset undergoes substantial, lifelong epigenetic changes [START_REF] Goronzy | Epigenetics of T cell aging[END_REF][START_REF] Tserel | Age-related profiling of DNA methylation in CD8+ T cells reveals changes in immune response and transcriptional regulator genes[END_REF]. Together, our findings provide strong statistical evidence that DNA methylation variation with age rely on different, non-mutually exclusive mechanisms: the progressive decline of the epigenetic maintenance system, common to all cell types, and the increased heterogeneity of immune cell subsets that characterizes immunosenescence [START_REF] Nikolich-Zugich | The twilight of immunity: emerging concepts in aging of the immune system[END_REF].

Sex Differences in DNA Methylation are Predominantly Cell-and Age-Independent

Given that substantial differences in immune cell composition between sexes have been observed [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF], we next assessed how cellular heterogeneity contributes to sex differences in DNA methylation [START_REF] Singmann | Characterization of whole-genome autosomal differences of DNA methylation between men and women[END_REF][START_REF] Yousefi | Sex differences in DNA methylation assessed by 450 K BeadChip in newborns[END_REF]. We found ~29% of CpG sites (n = 186,545) with a significant total effect of sex, ~20% (n = 126,904) with a significant direct effect, and ~7% (n = 44,667) with a significant cell-composition-mediated effect (FDR < 0.05; Table 2 andFigure S4A). The largest direct effects of sex were observed at DYRK2, DNM1, RFTN1, HYDIN, and NAB1 genes (P adj < 1.0×10 -285 ). For example, the DYRK2 promoter is 11% and 45% methylated in men and women, respectively, at the CpG site bound by the X-linked PHF8 histone demethylase (Figure S4B,C). DYRK2 phosphorylates amino acids and plays a key role in breast and ovarian cancer development [START_REF] Correa-Saez | Updating dual-specificity tyrosine-phosphorylation-regulated kinase 2 (DYRK2): molecular basis, functions and role in diseases[END_REF].

DNA methylation levels were higher in women at 78% of sex-associated autosomal CpG sites (Figure S4D,E), a pattern also observed in newborns [START_REF] Yousefi | Sex differences in DNA methylation assessed by 450 K BeadChip in newborns[END_REF]. This proportion is similar across different genomic regions, based on either chromatin states or CpG density (Figure S4E,I). When quantifying how sex differences in DNA methylation vary during adulthood, by adding a sex-by-age interaction term to our models (STAR Methods), we found only 23 CpG sites with a significant, sex-dependent effect of age (FDR < 0.05; Table S5), confirming previous findings [START_REF] Mccartney | An epigenome-wide association study of sex-specific chronological ageing[END_REF][START_REF] Yusipov | Age-related DNA methylation changes are sex-specific: a comprehensive assessment[END_REF]. The most associated genes are FIGN, associated with risk-taking behaviors [START_REF] Karlsson Linner | Genome-wide association analyses of risk tolerance and risky behaviors in over 1 million individuals identify hundreds of loci and shared genetic influences[END_REF] and educational attainment [START_REF] Lee | Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals[END_REF], and PRR4, associated with the dry eye syndrome, a hormone-dependent, late-onset disorder [START_REF] Perumal | Proteomics analysis of human tears from aqueous-deficient and evaporative dry eye patients[END_REF]. Overall, our findings indicate that the blood DNA methylome is widely affected by sex, but its effects are typically not mediated by cellular composition and do not change during adulthood.

Cytomegalovirus Infection Alters the Blood DNA Methylome through Regulation of Host Transcription Factors

We next leveraged the extensive questionnaire and phenotyping conducted in the Milieu Intérieur cohort to identify environmental factors that elicit cell-composition-independent changes in the blood methylome. Specifically, we estimated how 5mC levels are influenced by 139 variables (Table S1), including factors related to upbringing (e.g., birth weight, delivery route, rural or urban childhood), socio-economic status (e.g., educational attainment, monthly income, work-hours), dietary habits (e.g., eating frequency of various foods), health-related habits (e.g., smoking, BMI, physical exercise), lipid metabolism (e.g. low-and high-density lipoproteins, cholesterol and triglyceride levels), mental health and sleeping habits (e.g., self-reported depression, hours of sleep), exposure to pollutants (e.g., asbestos, benzene, silica), reproductive life cycle and contraception in women (e.g., contraceptive use, age at menopause), past and present exposure to infectious agents (e.g., cytomegalovirus, Epstein-Barr virus), total serum antibody concentrations (e.g. levels of IgG, IgE and IgM), and vaccination history (e.g., MMR, hepatitis A vaccine). Tests between each of the variables and 5mC levels at all measured CpG sites were considered a separate family and were adjusted to control the FDR at 0.05 (STAR Methods). All models were adjusted for associated meQTLs, genetic ancestry, batch variables and factors that impact cell composition, including sex and a non-linear age term.

The factor that is associated with the largest number of CpG sites is CMV serostatus. CMV is the causative agent of a latent, mainly asymptomatic, infection with prevalence ranging from 40% to 100% [START_REF] Cannon | Review of cytomegalovirus seroprevalence and demographic characteristics associated with infection[END_REF], which drastically alters the composition of the CD8 + T cell compartment in blood [START_REF] Klenerman | T cell responses to cytomegalovirus[END_REF]. CMV seropositivity has a significant total effect on ~36% of CpG sites (n = 233,014; Figure 3A andTable 2). When adjusting for blood cell 12 composition, a significant direct effect was detected for ~10% of CpG sites (n = 64,383; FDR < 0.05). Of note, the 16 cell proportions we adjusted for include central, effector memory and EMRA CD8 + T cells, which we have previously shown to be strongly associated with CMV serostatus [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. When we used the standard deconvolution method to correct for cellular heterogeneity, which does not include estimates of CD8 + sub-compartments [START_REF] Koestler | Improving cell mixture deconvolution by identifying optimal DNA methylation libraries (IDOL)[END_REF], we found twice as many CpG sites directly associated with CMV serostatus (~19%, n = 120,024), indicating again that this standard correction for cellular heterogeneity is not complete.

One of the strongest direct effects of CMV infection was observed nearby the TSS of LTBP3 (β value scale 95% CI: [3.0%, 4.2%], P adj = 2.2×10 -33 ; Figure 3E). LTBP3 is a regulator of latent transforming growth factor β (TGF-β) [START_REF] Morita | Egr2 and Egr3 in regulatory T cells cooperatively control systemic autoimmunity through Ltbp3-mediated TGF-beta3 production[END_REF], which is induced in CMV latently infected cells [START_REF] Mason | Human cytomegalovirus latency alters the cellular secretome, inducing cluster of differentiation (CD)4+ T-cell migration and suppression of effector function[END_REF]. We found that effects of CMV are typically smaller than those of age and sex (Figure S2) and are associated with an increase in 5mC levels in 92% of CpG sites within CGIs and a decrease in 76% of CpGs outside CGIs (Figure S5A,B). As for age, we observed in CMV + donors an overall increase in 5mC levels in Polycomb-repressed regions and binding sites of PRC-related TFs, and a decrease in regions of strong transcription (Figure 3B-D), suggesting dysregulation of the host gene inactivation machinery as a result of latent infection. Similar results were found when restricting the analysis to CpG sites outside of CGIs (Polycomb-repressed regions OR = 2.9, CI: [1.9, 4.24]). Interestingly, CpG sites showing increased 5mC levels in CMV + donors are strongly enriched in binding sites for the BRD4 TF (OR = 13.5, CI: [11.9, 15.3], P adj < 1.0×10 - 320 ; Figure 3D and Table S3), a bromodomain protein that plays a critical role in the regulation of latent and lytic phases of CMV infection [START_REF] Groves | Bromodomain proteins regulate human cytomegalovirus latency and reactivation allowing epigenetic therapeutic intervention[END_REF]. In addition, CpG sites showing a decrease in DNA methylation in CMV + donors are strongly enriched in binding sites for BATF3 (OR = 9.0, CI: [8.4, 9.7], P adj < 1.0×10 -320 ; Figure 3F and Table S3), which is paramount in the priming of CMV-specific CD8 + T cells by cross-presenting dendritic cells [START_REF] Torti | Batf3 transcription factordependent DC subsets in murine CMV infection: differential impact on T-cell priming and memory inflation[END_REF].

We investigated if the large shift in the composition of the CD8 + T cell compartment caused by CMV [START_REF] Klenerman | T cell responses to cytomegalovirus[END_REF] is accompanied by changes in 5mC levels. We found that 33% of CpG sites show a significant cell-composition-mediated effect of CMV serostatus (n = 217,223, FDR < 0.05). Importantly, 93% of CpG sites with a significant direct effect also show a significant mediated effect (n = 60,194; Figure 3A), and we observed a clear correlation between total and mediation effect sizes (R = 0.93; Figure 3A). For example, the CpG site in the TSS of LTBP3 with a large direct effect of CMV has also a large, mediated effect (CI: [2.4%, 3.5%], P adj = 1.8×10 -18 ); CMV + donors show higher proportions of CD8 + T EMRA cells (P = 1.38×10 -49 ), which in turn are associated with higher 5mC levels at LTBP3 (P = 6.9×10 -92 ), supporting mediation by this T cell subset (Figure 3G,H). Collectively, our analyses indicate that CMV infection affects a large fraction of the human blood DNA methylome through the dysregulation of host TFs and finegrained changes in cellular composition.

Strong Effects of Smoking are Reversible and Independent of Blood Cell Composition

The second exposure that is associated with the largest number of CpG sites is cigarette smoking, for which the total effect was significant for 7,257 CpG sites (~1.1%, Figure 3I and Table 2).

Although active smoking is known to elicit reproducible changes in DNA methylation [START_REF] Dugue | Smoking and blood DNA methylation: an epigenomewide association study and assessment of reversibility[END_REF][START_REF] Gao | DNA methylation changes of whole blood cells in response to active smoking exposure in adults: a systematic review of DNA methylation studies[END_REF], we previously showed that smoking also has a broad effect on blood immune cell subsets [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF], suggesting possible mediation by cellular composition.

When adjusting for the 16 immune cell proportions, we found that smoking directly alters 5mC levels at 2,416 CpG sites (~0.4% of CpG sites; FDR < 0.05; Table 2), 62% of which show a decrease in 5mC levels. For example, smokers show strongly decreased 5mC levels in the introns of the dioxin receptor repressor gene AHRR (β value scale 95% CI: [-23%, -20%], P adj = 9.9×10 -128 ), the second exon of F2RL3 (CI: [-10%, -8.6%], P adj = 3.2×10 -78 ) and the first intron of RARA (CI: [-10%, -8.5%], P adj = 5.4×10 -67 ), in agreement with previous studies [START_REF] Dugue | Smoking and blood DNA methylation: an epigenomewide association study and assessment of reversibility[END_REF][START_REF] Gao | DNA methylation changes of whole blood cells in response to active smoking exposure in adults: a systematic review of DNA methylation studies[END_REF]. No clear differences in the distribution of effect sizes were observed over genomic regions (Figure S5E,F). CpG sites that are demethylated in smokers are significantly enriched in binding sites for the hypoxia-related TFs EPAS1 and HIF2A (P adj < 7.0×10 -6 ), as well as AHRR (OR = 7.22, CI: [3.91, 12.2], P adj = 2.4×10 -8 ; Figure 3J and Table S3). This indicates that AHRR up-regulation in smokers elicits decreased 5mC levels at AHRR binding sites. To determine if such direct effects are reversible, we compared, for all smoking-associated CpG sites, the changes in 5mC levels with years since last smoke for past smokers, to the changes with years of smoking for active smokers. In agreement with previous studies [START_REF] Dugue | Smoking and blood DNA methylation: an epigenomewide association study and assessment of reversibility[END_REF][START_REF] Gao | DNA methylation changes of whole blood cells in response to active smoking exposure in adults: a systematic review of DNA methylation studies[END_REF], we found a strong negative correlation between effect sizes (R = -0.70; slope = -1.12; Figure S5I), supporting the reversibility of the direct effects of smoking on DNA methylation.

We estimated that 5mC levels are significantly altered by smoking due to changes in cell composition at ~3.2% of CpG sites (n = 20,381 CpG sites with a mediated effect; FDR < 0.05).

Among the most strongly affected CpG sites, we found IL18RAP (CI: [0.85%, 1.48%], P adj = 1.23×10 -9 ), a subunit of the receptor for IL18 that is differentially expressed by NK cells [START_REF] Crinier | High-Dimensional Single-Cell Analysis Identifies Organ-Specific Signatures and Conserved NK Cell Subsets in Humans and Mice[END_REF]. We observed that active smoking induces a reduction in the proportion of NK cells (P = 1.17×10 -10 ), which is in turn associated with lower 5mC levels at IL18RAP (P = 2.77×10 -29 ; Figure 3K, L), in line with an effect of smoking mediated by NK cells. Of note, mediated effects of smoking on 5mC levels were also reversible, to a degree similar to that of direct effects (R = -0.69; slope=-0.84, Figure S5J). Importantly, only 5.3% of CpG sites with a direct effect of smoking status also have a significant mediation effect, and, on average, mediated effects of smoking are weaker than direct effects (Figure 3I). Out of the 50 CpG sites with the largest total effect of smoking status (P adj < 2.4×10 -18 ), 49 CpG sites, including those in AHRR, showed no significant mediation effect (FDR < 0.05). Collectively, these findings indicate that the largest effects of cigarette smoking on the blood DNA methylome are reversible and independent of blood cell composition.

Other Environmental Exposures do not Trigger Strong, Widespread Changes in the Adult

DNA Methylome

The third environmental exposure that we identified as affecting DNA methylation variation is circulating levels of C-reactive protein (CRP), a marker of chronic, low-grade inflammation in healthy adults. Associations between CRP levels and hundreds of 5mC marks have been detected [START_REF] Ligthart | DNA methylation signatures of chronic lowgrade inflammation are associated with complex diseases[END_REF], but the strong relationship of CRP levels with the immune system [START_REF] Sproston | Role of C-Reactive Protein at Sites of Inflammation and Infection[END_REF] and genetic variation [START_REF] Ligthart | Genome Analyses of >200,000 Individuals Identify 58 Loci for Chronic Inflammation and Highlight Pathways that Link Inflammation and Complex Disorders[END_REF] suggests that that these factors could confound associations. Specifically, changes in blood cell composition may be the cause of changes in CRP levels, and this could induce spurious associations, instead of mediated effects, at CpG sites associated with immune cell proportions.

We found an association between CRP and 5mC levels at 20,043 CpG sites (~3.1% of CpG sites; FDR < 0.05; Table 2), a figure that, when adjusting for cellular composition, dropped to only 480, of which 80% (n = 386) showed decreased 5mC levels with increased CRP levels. We detected a CpG site within an enhancer nearby BCL2, a key regulator of apoptosis and inflammation [START_REF] Chong | Noncanonical Cell Fate Regulation by Bcl-2 Proteins[END_REF], where 5mC levels increase with increasing CRP levels (β value scale 95% CI: [0.6%, 1%], P adj = 1.06×10 -5 ; Figure S5K). Another example is a CpG site within an enhancer nearby ABCG1 (CI: [0.4%, 0.8%], P adj = 1.20×10 -5 ; Figure S5L). In our cohort, 5mC levels at the same site are also associated to triglyceride (CI: [1.2%, 2.1%], P adj = 2.42×10 -8 ) and HDL (CI: [-3.9, -2.1], P adj = 9.35×10 -9 ) levels. The associations were retained in a model including CRP, HDL and triglyceride levels, indicating that they affect ABCG1 5mC levels independently. CRP is known to inhibit cellular cholesterol efflux by downregulating ABCG1 mRNA levels, which are impaired in patients with type 2 diabetes, obesity, and hypertension [START_REF] Li | The transcription levels of ABCA1, ABCG1 and SR-BI are negatively associated with plasma CRP in Chinese populations with various risk factors for atherosclerosis[END_REF]. These results indicate that associations between CRP levels and DNA methylation are mainly, but not exclusively, due to changes in blood cell composition, and generate new hypotheses on the epigenetic mechanisms relating subclinical inflammation to metabolic conditions.

Besides CMV infection, smoking status and chronic inflammation, we found limited evidence of a direct effect of non-heritable factors on the DNA methylome of healthy adults (Table 2). We found a significant total effect of heart rate, ear temperature and hour of blood draw on 5mC levels at 76,018, 59,728 and 38,884 CpG sites, respectively (FDR < 0.05), but no associations remained significant when adjusting for cellular heterogeneity. In total, we found 59 significant cellcomposition-independent associations between the remaining non-heritable factors and 5mC levels (Table 2), the majority of which relate ABCG1, DHCR24 and CPT1A genes with lipid-related traits and BMI [START_REF] Braun | Epigenome-wide association study (EWAS) on lipids: the Rotterdam Study[END_REF]. In addition, we found a single association between log protein levels and 5mC levels at a CpG site close to DAO, a gene encoding D-amino acid oxidase involved in protein catabolism (95% CI: [3.4%, 7%], P adj = 0.015). We detected an association close to TCERG1L and educational attainment (95% CI: [-0.028%, -0.013%], P adj = 0.019). Genetic variation in TCERG1L is associated with years of education [START_REF] Lee | Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals[END_REF], but not with 5mC levels in our cohort. We also found a significant association between log uric acid levels and 5mC levels at the SLC2A9 gene (cg00071950; P adj = 0.0034), which is no longer significant when adjusting on the local meQTL SNP (P adj = 1.0), illustrating how DNA sequence variation can confound EWAS results. Nutritional habits, assessed based on 20 dietary frequency variables, have no detectable effects on the blood DNA methylome, except the frequency of raw fruit consumption at GLI2 (95% CI: [-2.5%, -1.1%], P adj = 0.0022). Of note, we did not replicate previously reported associations between DNA methylation and serum IgE levels [START_REF] Ek | Epigenome-wide DNA methylation study of IgE concentration in relation to selfreported allergies[END_REF][START_REF] Liang | An epigenome-wide association study of total serum immunoglobulin E concentration[END_REF] and did not detect any association with current socio-economic status. Collectively, these results indicate that environmental exposures related to upbringing, socio-economic status, nutrition or vaccination do not induce strong changes of the blood DNA methylome in our cohort of healthy adults.

Gene × Environment and Gene × Cell Type Interactions Affect DNA Methylation Variation

Gene × environment interactions are thought to underlie adaptable human responses to environmental exposures through epigenetic changes [START_REF] Feinberg | The Key Role of Epigenetics in Human Disease Prevention and Mitigation[END_REF]. Having established that age, sex, CMV serostatus, smoking status and chronic low-grade inflammation (CRP levels) are the main non-heritable determinants of DNA methylation variation, we evaluated whether their effects are genotype-dependent. We thus tested for genotype × age, genotype × sex or genotype × exposure interactions, adjusting for 16 measured cell proportions (STAR Methods). We found evidence of genotype-dependent effects at 175, 41, 4, 29 and 0 CpG sites for age, sex, smoking status, CMV serostatus and CRP levels, respectively (P adj < 0.05, MAF > 0.10; Figure 4A; Table S5), the interacting SNP being local in all except 7 cases. We detected a strong genotype × age interaction for three CpG sites located in the BACE2 gene, the 5mC levels of which decrease with age only in donors carrying the nearby rs2837990 G>A allele (β value scale 95% CI: [11%, 13%], P adj = 2.83×10 -10 ; Figure 4B; Table S5). BACE2 encodes beta-secretase 2, one of two proteases involved in the generation of amyloid beta peptide, a critical component in the etiology of Alzheimer's disease [START_REF] Holler | BACE2 expression increases in human neurodegenerative disease[END_REF].

We then explored whether genetic variants affect 5mC levels specifically in different immune cell types, i.e., cell-type-dependent meQTLs. Because inferences were unstable for rare immune cell subsets, we estimated the effects of associated variants within six major cell types, compared to the effect of the variants within neutrophils -the most frequent blood cell subset. We found that genotypes affected DNA methylation differently according to cellular composition at 695 CpG sites.

We found 264, 157, 62, 56, 32, and 19 significant interaction effects for CD4 + T cells, CD8 + T cells, NK cells, B cells, CD4 -CD8 -T cells and monocytes, respectively (P adj < 0.05; Figure 4A). One of the strongest signals was found between 5mC levels at the TSS of CD300A and the nearby rs12939435 variant, the effects of which depend on the proportion of CD8 + T cells (CI: [-0.50%, -0.29%], P adj = 2.19×10 -14 ; Figure 4C; Table S5). CD300A is an immunomodulatory molecule that is expressed in various immune cell types and is associated with a cytotoxic molecular signature in CD8 + T cells (Xu et al., 2012b). Overall, our analyses identify several environment-and cell-typedependent meQTLs, supporting a strong, but limited impact of gene × environment and gene × cell type interactions on the blood DNA methylome.

Genetics and Cellular Heterogeneity Drive DNA Methylation Variation in Human Blood

Having established how genetic variation, cellular composition, intrinsic factors and a broad selection of non-heritable factors shape the blood DNA methylome, we next sought to compare the relative impact of these factors on DNA methylation. We classified the factors into four groups: (i) the cellular heterogeneity group, which consists of the 16 measured cell proportions; (ii) the intrinsic group, which consists of age and sex; (iii) the genetic group, which consists of the most associated local-meQTL variant around each CpG site; and (iv) the exposure group, which consists of smoking status, CMV serostatus and chronic low-grade inflammation. Since these groups vary in their degrees of freedom, we measured the relative predictive strength for each CpG site by the outof-sample prediction accuracy, estimated by cross-validation (STAR methods). To ensure unbiased estimates, we mapped local meQTLs anew within each training set.

The model explains < 5% of out-of-sample variance for 51% of CpG sites (Figure 5A), which are typically characterized by low total 5mC variance (Figure S6A). This suggests that these sites are constrained in the population and that small fluctuations in 5mC levels determine their variation, possibly due to measurement errors or biological noise. Nevertheless, the model explains > 25% of DNA methylation variance for 21% of CpG sites (n = 133,180). The strongest predictor for these CpGs is cellular composition, genetics, intrinsic factors and exposures in 74%, 22%, 4% and 0% of cases, respectively. Cellular composition explains > 25% of out-of-sample variance for 13% of CpG sites (n = 86,046; Figure 5A, C and Table S6), with the highest variance explained by cellular composition for one CpG site being 68.5%. The 16,034 CpG sites for which > 50% of variance is explained by cellular composition are typically located in genes related to the immune system (Top 3 gene ontology terms: leukocyte activation, cell activation, cell activation involved in immune response, P adj < 1.0×10 -28 ; Table S4). These CpG sites are concentrated in enhancer regions (95% CI: [3.41 3.69]; Figure S6B), and largely depleted from TSS (95% CI: [0.0761 0.106]), reflecting the importance of enhancer DNA methylation in cell-type identity.

For the 2,521 CpG sites where the model explains > 75% of variance, local genetic variation is the strongest predictor in 99% of cases (Figure 5C and Table S6). Local genetic variation explains > 25% of DNA methylation variance at 23,796 CpG sites, and almost as many when adjusting for cellular composition (n = 23,062) (Figure 5A,B), indicating that genetic effects on 5mC levels are cell-composition-independent. Intrinsic factors explain > 25% of out-of-sample variance at 3,621 CpG sites, and > 75% at 17 sites (Figure 5C). When conditioning on cell composition, these numbers dropped to 379 and 7 CpG sites, respectively, suggesting that the predictive ability of age and sex is partly mediated by immune cell composition (Figure 5B). Interestingly, environmental exposures are the weakest predictor of 5mC levels, explaining > 25% of the variance at only 23 CpG sites and with a maximum variance explained for a CpG site of 51%.

Finally, we estimated the proportion of variance explained by genotype × age, genotype × sex and genotype × exposure interactions, by considering the difference of the out-of-sample variance explained by models including interaction terms and models with only main effects (STAR Methods). We found a significant increase in predictive ability when including interaction terms for 1,984 CpG sites (ANOVA P adj < 0.05). However, the effects were typically modest: only 35 CpG sites showed an increase in the proportion of variance explained larger than 4% (Figure 5B). The largest difference was found for a CpG site in the TSS of the ENOSF1 gene, where the interaction model explained an additional 11.1% of DNA methylation variance (Table S6). Collectively, these results show that cellular composition and local genetic variation are the main drivers of DNA methylation variation in the blood of adults, reinforcing the critical need to study epigenetic risk factors and biomarkers of disease in the context of these factors.

DISCUSSION

Here, we present a rich data resource that delineates the contribution of genetics, age, sex, environmental factors, cellular composition and their interactions to variation in the DNA methylome. All the results can be explored via a web-based browser (MIMETH browser), to facilitate the exploration of the estimated effects of these factors on DNA methylation variation. We show that genetic variation controlling 5mC levels is likely to affect phenotype variation and disease risk, and often controls the expression of TFs. Furthermore, the remote genetic control of DNA methylation is driven by variants nearby ZNF genes, consistent with a role of ZFPs as direct regulators of 5mC levels [START_REF] Marchal | Emerging concept in DNA methylation: role of transcription factors in shaping DNA methylation patterns[END_REF]. Furthermore, we show that remote meQTLs preferentially affect 5mC levels of ZNF genes, supporting the view that the major targets of ZFPsmediated regulation are ZNF genes themselves [START_REF] O'geen | Genome-wide analysis of KAP1 binding suggests autoregulation of KRAB-ZNFs[END_REF]. Most ZFPs possess a Krüppelassociated box (KRAB) domain, a DNA-binding domain that elicits KAP1-mediated transcriptional repression and induce heterochromatin by recruiting chromatin remodelers and DNA methyltransferases [START_REF] Quenneville | The KRAB-ZFP/KAP1 system contributes to the early embryonic establishment of sitespecific DNA methylation patterns maintained during development[END_REF][START_REF] Vogel | Human heterochromatin proteins form large domains containing KRAB-ZNF genes[END_REF][START_REF] Zuo | Zinc finger protein ZFP57 requires its co-factor to recruit DNA methyltransferases and maintains DNA methylation imprint in embryonic stem cells via its transcriptional repression domain[END_REF], providing a putative mechanism for the direct regulation of DNA methylation by KRAB-ZFPs. This is also supported by the widespread effect of a SENP7 regulatory variant on 5mC levels of a KRAB-ZNF gene cluster on chromosome 19; SENP7 is a SUMO protease involved in the deSUMOylation of KAP1 that allows its chromatin remodelling activity [START_REF] Garvin | The deSUMOylase SENP7 promotes chromatin relaxation for homologous recombination DNA repair[END_REF].

Our study reveals three different biological mechanisms underlying age-related changes in DNA methylation. The first elicits increased 5mC variance with age and is related to epigenetic drift [START_REF] Fraga | Epigenetic differences arise during the lifetime of monozygotic twins[END_REF][START_REF] Jones | DNA methylation and healthy human aging[END_REF], likely caused by the progressive decline in fidelity of the DNA methylation maintenance machinery. The second produces cell-composition-independent, global DNA demethylation and CGI-associated hypermethylation. Age-associated DNA demethylation could be related to the downregulation of DNMT3A/B de novo methyltransferases, whereas CGI-associated hypermethylation may result from the downregulation of the Polycomb repressive complexes 1 and 2 and/or TET proteins, coupled with a loss of H3K27me3 marks [START_REF] Beerman | Proliferation-dependent alterations of the DNA methylation landscape underlie hematopoietic stem cell aging[END_REF][START_REF] Li | Genome-wide analyses reveal a role of Polycomb in promoting hypomethylation of DNA methylation valleys[END_REF][START_REF] Williams | DNA methylation: TET proteins-guardians of CpG islands?[END_REF]. Alternatively, these changes may be related to the mitotic clock, which assumes a progressive accumulation of DNA methylation changes with mitotic divisions, including loss of methylation at partially methylated domains (PMD) and gain of methylation at PRC2-marked CpG-rich regions [START_REF] Kim | Counting human somatic cell replications: methylation mirrors endometrial stem cell divisions[END_REF][START_REF] Yang | Correlation of an epigenetic mitotic clock with cancer risk[END_REF]Zhou et al., 2018). Both scenarios are supported by the enrichment of Polycomb-repressed regions in ageassociated CpG sites, and of binding sites of PRC-related TFs in CpG sites methylated with age.

The third mechanism elicits cell-composition-mediated demethylation at all compartments of the epigenome, particularly at enhancers of myeloid activation genes. This process likely reflects an increased degree of differentiation in the lymphoid compartment. Single-cell methylomes of differentiating and dividing white blood cells will help determine the role of mitotic and post-mitotic 5mC changes during epigenetic aging.

Latent infections are known to profoundly alter the number, activation status and transcriptional profiles of immune cell populations, yet their epigenetic consequences have attracted little attention.

We found that CMV infection elicits widespread changes in the blood DNA methylome, in contrast with other herpesviruses such as EBV, HSV-1, HSV-2 and VZV. We observe that most CMV effects are mediated by the profound changes in blood cell composition caused by CMV [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF], including the inflation of CMV-specific memory CD8 + T cells [START_REF] Klenerman | T cell responses to cytomegalovirus[END_REF]. However, we also detected cell-composition-independent effects of CMV infection, suggesting that the herpesvirus can directly regulate the host epigenome. Methylated CpG sites in CMV + donors are targeted by BRD4, a key host regulator of CMV gene expression and latency [START_REF] Groves | Bromodomain proteins regulate human cytomegalovirus latency and reactivation allowing epigenetic therapeutic intervention[END_REF], suggesting that this TF, when upregulated during latent CMV infection, binds both viral and host genomes. Furthermore, CMV + donors are characterized by a strong increase in 5mC levels at LTBP3, the product of which is involved in TGF-β secretion. TGF-β is a well-known immunosuppressive cytokine induced by CMV infection [START_REF] Mason | Human cytomegalovirus latency alters the cellular secretome, inducing cluster of differentiation (CD)4+ T-cell migration and suppression of effector function[END_REF], which represents a possible strategy of the virus to escape host immunity. These results suggest that the capacity of CMV to manipulate the host epigenetic machinery results in epigenetic changes of latently infected cells.

Another interesting finding of our study is that environmental exposures explain a small fraction of the variance of DNA methylation in healthy adults, at odds with the common view that the epigenome is strongly affected by the environment [START_REF] Feil | Epigenetics and the environment: emerging patterns and implications[END_REF]. Twin studies have estimated the heritability of DNA methylation to range from ~20-40% [START_REF] Bell | Epigenome-wide scans identify differentially methylated regions for age and age-related phenotypes in a healthy ageing population[END_REF][START_REF] Grundberg | Mapping cis-and trans-regulatory effects across multiple tissues in twins[END_REF][START_REF] Van Dongen | Genetic and environmental influences interact with age and sex in shaping the human methylome[END_REF], suggesting that environmental effects, along with gene × environment interactions, account for the remaining 60-80% [START_REF] Teschendorff | Statistical and integrative system-level analysis of DNA methylation data[END_REF]. However, other factors, including cellular composition and measurement error, may account for most of the unexplained variance. Consistently, we estimated that cellular composition explains >25% of the variance for ~13% of the DNA methylome, and it has been estimated that measurement error may explain >50% [START_REF] Li | Causes of blood methylomic variation for middle-aged women measured by the HumanMethylation450 array[END_REF]. Nevertheless, a limitation of our study is that perinatal and early life exposures, which are thought to contribute extensively to epigenetic variation in adulthood [START_REF] Feil | Epigenetics and the environment: emerging patterns and implications[END_REF], have not been extensively assessed in the Milieu Intérieur cohort. In addition, it has been hypothesized that gene × environment interactions are central to understand the role of epigenetics in development [START_REF] Boyce | Development and the epigenome: the 'synapse' of geneenvironment interplay[END_REF], but statistical evidence for interaction effects requires larger cohorts [START_REF] Fleiss | Design and analysis of clinical experiments[END_REF], suggesting that our results might represent the small, perceptible fraction of a large number of weak effects [START_REF] Czamara | Integrated analysis of environmental and genetic influences on cord blood DNA methylation in new-borns[END_REF][START_REF] Teh | The effect of genotype and in utero environment on interindividual variation in neonate DNA methylomes[END_REF]. Large, longi-tudinal cohorts addressing the developmental origins of disease are needed to shed new light on the role of DNA methylation in the interplay between genes and the environment.

Collectively, our findings have broad consequences for the study and interpretation of epigenetic factors involved in disease risk. First, a third of the DNA methylome is affected by genetic variants, some of which are associated with disease risk. Epigenetic associations with a given disease or trait may thus result from the pleiotropic effects of genetic variants on DNA methylation, which may confound interpretation. Second, because age, sex, CMV infection, smoking and chronic lowgrade inflammation influence disease risk [START_REF] Furman | Chronic inflammation in the etiology of disease across the life span[END_REF][START_REF] Mauvais-Jarvis | Sex and gender: modifiers of health, disease, and medicine[END_REF][START_REF] Niccoli | Ageing as a risk factor for disease[END_REF][START_REF] Samet | Tobacco smoking: the leading cause of preventable disease worldwide[END_REF][START_REF] Savva | Cytomegalovirus infection is associated with increased mortality in the older population[END_REF], our results highlight the critical need to consider such factors in EWAS. Third, our analyses clearly show that the effects of age, CMV serostatus and CRP levels are largely mediated by fine-grained changes in immune cell proportions. This reinforces the view that EWAS must be interpreted with caution, particularly when standard corrections using estimated cell proportions [START_REF] Houseman | DNA methylation arrays as surrogate measures of cell mixture distribution[END_REF][START_REF] Koestler | Improving cell mixture deconvolution by identifying optimal DNA methylation libraries (IDOL)[END_REF][START_REF] Teschendorff | A comparison of referencebased algorithms for correcting cell-type heterogeneity in Epigenome-Wide Association Studies[END_REF] are incomplete. The integration of DNA methylation profiling and finegrained measurements of immune cell subsets, such as the data used here, could also help improving the estimation of blood cell composition from DNA methylation and corrections for cellular heterogeneity. Finally, our findings highlight the major epigenetic impact of aging, persistent viral infections and inflammation through fine-grained changes in blood cell proportions, highlighting the need to assess the respective role of DNA methylation and altered cellular composition in the etiology of disease [START_REF] Lappalainen | Associating cellular epigenetic models with human phenotypes[END_REF]. Large-scale studies using single-cell approaches will help overcome these challenges, and are anticipated to further decode the epigenetic mechanisms underlying healthy aging and the environmental causes of human disease. Further information and requests for resources and information should be directed to and will be fulfilled by the Lead Contact, Pr. Lluís Quintana-Murci (quintana@pasteur.fr).

FIGURES AND LEGENDS

• EXPERIMENTAL MODEL AND SUBJECT DETAILS

The Milieu Intérieur cohort

The Milieu Intérieur cohort was established with the goal to identify genetic variation and environmental exposures that affect phenotypes related to the immune system in the adult, healthy population. The 1,000 healthy donors of the Milieu Intérieur cohort were recruited by BioTrial (Rennes, France), and included 500 women and 500 men. Donors included 100 women and 100 men from each decade of life, between 20 and 69 years of age. Donors were selected based on various inclusion and exclusion criteria that are detailed elsewhere [START_REF] Thomas | The Milieu Interieur study -an integrative approach for study of human immunological variance[END_REF]. Briefly, donors were required to have no history or evidence of severe/chronic/recurrent pathological conditions, neurological or psychiatric disorders, alcohol abuse, recent use of illicit drugs, recent vaccine administration, and recent use of immune modulatory agents. To avoid the influence of hormonal fluctuations in women, pregnant and peri-menopausal women were not included. To avoid genetic stratification in the study population, the recruitment of donors was restricted to individuals whose parents and grandparents were born in Metropolitan France.

Ethical approvals

The • METHOD DETAILS

DNA sampling and extraction

Whole blood samples were collected from the 1,000 Milieu Intérieur healthy, fasting donors on Liheparin, every working day from 8AM to 11AM, from September 2012 to August 2013, in Rennes,

France. Tracking procedures were established in order to ensure delivery to Institut Pasteur (Paris) within 6 hours of blood draw, at a temperature between 18°C and 25°C. Upon receipt, samples were kept at room temperature until DNA extraction. DNA was extracted using the Nucleon BACC3 genomic DNA extraction kit (GE Healthcare, Illinois, USA). High-quality genomic DNA was obtained for 978 out of the 1,000 donors.

DNA methylation profiling and data quality controls

Extracted genomic DNA was treated with sodium bisulfite (Zymo Research, California, USA).

Bisulfite-converted DNA was applied to the Infinium MethylationEPIC BeadChip (Illumina, California, USA), using the manufacturer's standard conditions. The MethylationEPIC BeadChip measures 5mC levels at 866,836 CpG sites in the human genome. Raw IDAT files were processed with the minfi R package [START_REF] Fortin | Preprocessing, normalization and integration of the Illumina HumanMethylationEPIC array with minfi[END_REF]. All samples showed average detection P-values lower than 0.005. No sample showed a mean of methylated intensity signals lower than 3 standard deviations from the cohort average. Thus, no samples were excluded based on detection P-values or methylated intensity signals. The sex predicted from 5mC signals on sex chromosomes matched the declared sex for all samples (Figure S1A). Using the 59 control SNPs included in the MethylationEPIC array, a single sample showed high genotype discordance with the genome-wide SNP array data (see 'Genome-wide DNA genotyping' section) and was thus excluded (Figure S1B).

Unmethylated and methylated intensity signals were converted to M-values. A total of 2,930 probes with >1% missingness (i.e., detection P-value > 0.05 for more than 1% of donors) were excluded and remaining missing data (missingness = 0.0038%) were imputed by mean substitution. Using the irlba R package, Principal Component Analysis (PCA) of M values identified nine outlier samples, including eight that were processed on the same array (Figure S1C), which were also excluded. The "noob" background subtraction method [START_REF] Triche | Low-level processing of Illumina Infinium DNA Methylation BeadArrays[END_REF] was applied on M values for the remaining 969 samples, which showed highly consistent epigenome-wide 5mC profiles (Figure S1D,E).

To identify batch effects on the DNA methylation data, we searched for the factors that were the most associated with the top 20 PCs of the PCA of noob-corrected M values. We used a linear mixed model that included the proportion of lymphocytes, age, sex and cytomegalovirus (CMV) serostatus as fixed effects, and slide position and sample plate as random effects. Strong associations were observed between the first four PCs and slide position and sample plate (Figure S1F,G). M values were thus corrected for these two batch effects using ComBat [START_REF] Johnson | Adjusting batch effects in microarray expression data using empirical Bayes methods[END_REF]. After ComBat correction, the ten first PCs of a PCA of M values were associated with factors known to affect DNA methylation, including blood cell composition, age and sex (Figure S1H-J), indicating no other, strong batch effect on the data. M-values were converted to β values, considering that β = 2 M / (2 M + 1). Because outlier 5mC values due to measurement error could inflate the type I error rate of regression models, we excluded, for each CpG site, M or β values that were greater than 5 × standard deviations (SD) from the population average, corresponding to <0.1% of all measures. We also excluded (i) 83,380 non-specific probes that share >90% sequence identity with several genomic regions (see details in [START_REF] Price | Additional annotation enhances potential for biologically-relevant analysis of the Illumina Infinium HumanMethylation450 BeadChip array[END_REF]), (ii) 118,575 probes that overlap a SNP with MAF>1% in the Milieu Intérieur cohort or in European populations from the 1,000 Genomes project [START_REF] Auton | A global reference for human genetic variation[END_REF], (iii) 558 probes that were absent from the Illumina annotations version 1.0 B4 and (iv) 16,876 probes located on sex chromosomes. As a result, the final, quality-controlled data was composed of 968 donors profiled at 644,517 CpG sites.

Flow cytometry

Protocols, panels, staining antibodies and quality control filters used for flow cytometry analyses are detailed elsewhere [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. Briefly, immune cell proportions were measured using ten eight-color flow-cytometry panels. The acquisition of cells was performed using two MACSQuant analyzers, which were calibrated using MacsQuant calibration beads (Miltenyi, Germany). Flow cytometry data were generated using MACSQuantify software version 2.4.1229.1.

The mqd files were converted to FCS compatible format and analyzed by FlowJo software version 9.5.3. A total of 110 cell proportions were exported from FlowJo. Abnormal lysis or staining were systematically flagged by trained experimenters. We removed outliers by using a scheme detailed previously [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. Briefly, we used a distance-based approach that, for each cell-type, removes observations in the right tail if the distance to the closest observation in the direction of the mean is larger than 20% of the range of the observations. Similarly, observations in the left tail were removed if the distance to the closest observation in the direction of the mean is more than 15% than the range the observations. We removed 22 observations in total, including a maximum of 8 observations for a single cell type (i.e., for the proportion of neutrophils). Finally, missing data were imputed using the random forest-based missForest R package [START_REF] Stekhoven | MissForest--non-parametric missing value imputation for mixed-type data[END_REF].

Genome-wide DNA genotyping

Protocols and quality control filters for genome-wide SNP genotyping are detailed elsewhere [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. Briefly, all the 1,000 Milieu Intérieur donors were genotyped on both the HumanOmniExpress-24 and the HumanExome-12 BeadChips (Illumina, California, USA), which include 719,665 SNPs and 245,766 exonic SNPs, respectively. Average concordance rate between the two genotyping arrays was 99.9925%. The final data set included 732,341 high-quality polymorphic SNPs. After genotype imputation and quality-control filters, a total of 11,395,554

SNPs was further filtered for minor allele frequencies > 5%, yielding a data set composed of 1,000 donors and 5,699,237 SNPs for meQTL mapping. Ten pairs of first to third-degree related donors were detected with KING 1.9 [START_REF] Manichaikul | Robust relationship inference in genome-wide association studies[END_REF]. Out of the 968 donors whose blood methylome was profiled, 958 unrelated donors were kept for subsequent analyses.

• QUANTIFICATION AND STATISTICAL ANALYSIS

Circulating immune cells

One of the key questions in this study is whether differences in 5mC levels observed with respect to different factors are due to epigenetic changes occurring within cells or if they in fact reflect changes in cell composition. To answer this question, we adjusted models on measured proportions of 16 major subsets of blood: naïve, central memory (CM), effector memory (EM) and terminally differentiated effector memory (EMRA) subsets of CD4 + and CD8 + T cells, CD4 -CD8 -T cells, B cells, dendritic cells, natural killer (NK) cells, monocytes, neutrophils, basophils and eosinophils [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. We also investigated whether some factors affect 5mC levels differently within cellular subsets. To answer this question, we derived interaction models, where measured cell proportions interacted with the factor. Since these models showed inflated variance for small subsets, we used a reduced set of 7 immune cell types for this analysis: CD4 + and CD8 + T cells, CD4 -CD8 -T cells, B cells, NK cells, monocytes and neutrophils.

Local meQTL mapping

Local meQTL mapping was performed using the MatrixEQTL R package [START_REF] Shabalin | Matrix eQTL: ultra fast eQTL analysis via large matrix operations[END_REF].

Association was tested for each CpG site and each SNP in a 100-Kb window around the CpG site, by fitting a linear regression model assuming an additive allele effect. Models included the set of 16 immune cell proportions (see above) as predictors. They also included factors we have previously identified to have a large impact on blood and its molecular characteristics: a nonlinear age term encoded by 3 degrees-of-freedom (DoF) natural splines, sex, smoker status, ex-smoker status and CMV serostatus [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. We also adjusted for the top two PCs of a PCA of the genotype data. We did not include more PCs because of the low population substructure observed in the cohort [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. For the ݅:th individual and the :th CpG site, let ݕ be the measured 5mC levels on the M value scale, SNܲ the number of minor alleles of the m:th tested SNP for the CpG site and ݂ ఉ ಲ ሺ݁݃ܣ ሻ a nonlinear age term of natural splines with corresponding parameter vector ߚ . Moreover, let the vector ܿ be measurements of the 16 immune cell types for the i:th individual and ߚ be the corresponding parameter vector. The additive allele effect of the SNP was estimated by the parameter ߚ in the model,

ݕ ൌ ߤ ܵܰܲ ߚ ܿ ் ߚ ԛ ԛ݂ ఉ ಲ ሺ݁݃ܣ ሻ ܹ݊ܽ݉ ߚ ௐ ݎ݁݇݉ݏݔܧ ߚ ா௫௦ ݎ݁݇݉ܵ ߚ ௌ ܸܯܥ ߚ ெ 1ܥܲ ߚ ଵ 2ܥܲ ߚ ଶ ߝ ,
Eq. 1

Where ߝ is a symmetrical zero-mean distribution with constant variance.

Long-range meQTL mapping

Testing all possible associations between 644,517 CpG sites and 5,699,237 SNPs would require performing 3,769 billion statistical tests. To reduce the number of tests, long-range meQTL mapping was conducted on a selection of 50,000 CpG sites with the highest residual variance in the model described in Eq. 1, but with m indexing in this case the most associated local SNP for each site. For each of the 50,000 selected CpG sites, we then fitted one model per SNP located outside of a 1-Mb window around the CpG site. For each SNP-CpG pair, we estimated the additive allele effect of the remote SNP using the model specified in Eq. 1.

Local and long-range meQTL mapping were adjusted for multiple testing by employing a two-stage hierarchical procedure designed for the structure of tested hypotheses [START_REF] Peterson | Many Phenotypes Without Many False Discoveries: Error Controlling Strategies for Multitrait Association Studies[END_REF].

Consider all performed hypothesis tests, ܪ , ,  ൌ 1, … , ܰ, ݉ ൌ 1, … , ܯ , where ܰ is the number of CpG sites and ܯ is the number of SNPs considered for the :th CpG site, with corresponding P-values, ܲ , . Define the family of all hypothesis tests performed for the :th CpG site, ܪ ൌ ቄܪ ,ଵ , … , ܪ ,ெ ቅ.

For each such family, we tested the intersection hypothesis of no genetic control of 5mC levels at the CpG site. The P-value for this hypothesis was computed as the smallest Bonferroni-adjusted Pvalue in the set of P-values for the family, ܲ ൌ min ቄܯ ܲ , ቅ.

To adjust for multiplicity due to the number of CpG sites, we adjusted the P-value collection ܲ ,  ൌ 1, . . , ܰ by the Benjamini-Hochberg procedure. We considered an intersection hypothesis to be rejected if its Benjamini-Hochberg adjusted P-value was below 0.05. Let ܵ be the number of rejected intersection hypotheses. In the final stage, we performed hypothesis tests for association of SNPs with 5mC levels of CpG sites under genetic control, i.e., within families with a rejected intersection hypothesis. We considered a hypothesis ܪ , within a selected family to be rejected if,

ܯ ܲ , ܰ ܵ ൏ 0.05.
This procedure controls the false discovery rate (FDR) for discovery of CpG sites under genetic control, the global FDR over all tests and the average family-wise error rate (FWER) over selected families [START_REF] Peterson | Many Phenotypes Without Many False Discoveries: Error Controlling Strategies for Multitrait Association Studies[END_REF].

Detection of independent remote meQTLs

We designed the following scheme to compute a set Φ of independently associated remote SNPs for each CpG site, where all such SNPs are associated with 5mC levels ݕ at the :th CpG site conditional on the most associated local SNP and other SNPs in Φ. Define ܺ ଵ to be the set of SNPs with a long-range association to ݕ and let ݔ be the most associated significant local SNP, if it exists. The set ܺ ଵ includes many SNPs that are in linkage disequilibrium (LD). The algorithm uses an iterative procedure to build sets ܯ of SNPs, where in the ݆:th iteration, SNPs that are not associated with 5mC levels at the CpG site conditional on SNPs included in ܯ ିଵ are discarded, while the most associated is retained in ܯ ୨ . In the final step, the set Φ is constructed by elements of ܯ ୨ that are associated with 5mC levels at the CpG site conditional on all the other elements.

Intuitively, Φ consists of the most associated SNP in each LD block. The algorithm is given in pseudocode in Algorithm 1, where the condition ߚ ് 0 is determined by an F test on the level ߙ ൌ 10 ି .

Algorithm 1: Forming a set of long-range independently associated SNPs with a CpG site If the CpG site is under local genetic control then let ܯ ଵ ൌ ݔ , otherwise let ܯ ଵ ԛ ൌ ԛ

Repeat for

݆ ൌ 1, 2, … ܲ ൌ ሼݔ א ܺ ך ܯ : ԛߚ ௫ ് 0 in ݕ ൌ ߤ ߚݔ ௫ ∑ ߚݖ ௭אெ ೕ ε, ε  ሺ0, σ ଶ ሻሽ If ܲ ൌ  Exit ܺ ାଵ ൌ ܲ M ୨ାଵ ԛ ൌ ԛM ୨ ԛ  ሼ:ݔ ݔ SNP with the smallest P-value in Pሽ End Φ ൌ ሼݔ א ܯ ାଵ ך ݔ : ԛߚ ௫ ് 0 in ݕ ൌ ߤ ߚݔ ௫ ∑ ݖ ௭אெ ೕశభ ךሼ௫ሽ ߚ ௭ ߝ, ߝ  ሺ0, σ ଶ ሻ

Epigenome-wide association studies of non-genetic factors

We assessed the effect of 141 non-genetic variables (Table S1) on the blood methylome of adults.

The measured 5mC levels at a CpG site is an average of the methylation state of this CpG site in all cells in the blood sample. Cell composition is unlikely to have a strong causal effect on most of the investigated variables, with few exceptions, such as C-reactive protein (CRP) levels. However, many of the 141 candidate variables are likely to influence cell composition, which will cause a corresponding change in 5mC levels. We denote this effect the "(cell-type-)mediated effect". In addition, the variable might alter 5mC levels within individual cells, or within cell types. We denote this effect the "direct effect" (See Figure S3H for a schematic directed acyclic graph of the system).

Several important factors are known to have a large effect on blood cell composition in healthy donors, the most important being age, sex, CMV serostatus and smoking [START_REF] Patin | Natural variation in the parameters of innate immune cells is preferentially driven by genetic factors[END_REF]. As an added complexity, these factors are also associated with most of the other variables in the study.

Based on this framework, we investigated four questions, each one targeted by a separate statistical model.

The total effect

The total effect includes both changes in 5mC levels induced by changes in cellular composition and those induced within cell types. For each variable of interest ݔ and CpG site pair, the total effect was estimated in a regression model including 5mC levels of the CpG site on the M value scale as response variable and ,ݔ a nonlinear age term of 3 DoF natural splines, sex, CMV serostatus, smoking status, the most associated significant local SNP, independently associated remote SNPs and the two first PCs of the genotype matrix as predictors. In addition, since we noticed variability in 5mC levels across days of blood draw, we included date of blood draw as a random effect. Let ݆ be the day of blood draw for the ݅:th individual. For the p:th CpG site, let ݕ be the 5mC levels of the ݅:th individual on the M value scale, ݂ ఉ ಲ ሺ݁݃ܣ ሻ a nonlinear age term of 3 DoF natural splines and SNP a vector of the number of minor alleles of independently associated SNPs with corresponding parameter vector ߚ ௌே . The total effect of the variable ݔ was estimated by the corresponding parameter β ௫ in the model,

ݕ ൌ ߤ ݔ ߚ ௫ ݂ ఉ ಲ ሺAge ୧ ሻ Woma݊ ߚ ୭୫ୟ୬ Exsmokeݎ ߚ ୶ୱ୫୭୩ୣ୰ Smokeݎ ߚ ୗ୫୭୩ୣ୰ CMܸ ߚ େ PC1 ߚ େଵ PC2 ߚ େଶ SNܲ ் β ୗ DayOfSamplin݃ ሺሻ ε ,
Eq. 2

where DayOfSamplin݃ ሺሻ  ࣨሺ0, σ ௗ ଶ ሻ and ߝ  ሺ0, ߪ ଶ ሻ. Aging was tested by removing ݔ and replacing the non-linear age term with a linear one in Eq. 2. The effects of sex, smoking status and CMV serostatus were tested by removing ݔ in Eq. 2. For variables concerning women only (e.g., age of menarche), we excluded men from the analysis and removed the ܹ݊ܽ݉ ߚ ௐ term.

Hypothesis tests were performed by the Kenward-Roger approximation of the F-test for linear mixed models, implemented in the pbkrtest R package [START_REF] Halekoh | A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models -The R Package pbkrtest[END_REF].

The direct effect

Let the vector ܿ be measurements of the 16 immune cell types for the i:th individual and ߚ be the corresponding parameter vector. Using the same notation as for the total effect, the direct effect of the variable ݔ was estimated by ߚ ௫ in the model,

ݕ ൌ ߤ ݔ ߚ ௫ ܿ ் ߚ ݂ ఉ ಲ ሺ݁݃ܣ ሻ ܹ݊ܽ݉ ߚ ௐ ݎ݁݇݉ݏݔܧ ߚ ா௫௦ ݎ݁݇݉ܵ ߚ ௌ ܸܯܥ ߚ ெ 1ܥܲ ߚ ଵ 2ܥܲ ߚ ଶ ܵܰܲ ் ߚ ௌே ݈݂ܱ݃݊݅݉ܽܵݕܽܦ ߝ .
Eq. 3

For age and sex, age and CMV serostatus, and age and smoking status, we also estimated their interaction effect by including one interaction term at a time in the model specified in Eq. 3.

Hypothesis tests were performed by the Kenward-Roger approximation of the F-test for linear mixed models, implemented in the pbkrtest R package [START_REF] Halekoh | A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models -The R Package pbkrtest[END_REF].

The mediated effect

We estimated the mediated effect of aging, sex, variables related to smoking and CMV serostatus. It was estimated as the effect on 5mC levels mediated by changes in proportions of the 16 cell subsets due to the given factor. Estimates were computed by a two-stage procedure. Introduce the vector ݇ of covariates: age (an entry for each spline term), sex, smoking, CMV serostatus and ancestry (2 PCs), but excluding the variable of interest, ݔ (mediated effect of aging was estimated with a linear term), and let ܿ be a vector of measured proportions of the 16 blood subsets. We fitted two different groups of models. In the first, measured proportions of immune cells were response variables. For the model of the ݊:th cell type, let ߚ be the parameter vector for covariates ݇ and ߚ ௫ the parameter for the variable of interest. Let ܿ denote the n:th entry of the vector ܿ , the measured proportion of the n:th cell type for the i:th individual. For the model of 5mC levels in the M value scale at the :th CpG site, ݕ , let ߠ ௫ be a parameter for the variable of interest and ߠ and ߠ parameter vectors for the effects of cell proportions and covariates. To estimate mediated effects of the variable of interest ,ݔ we fit the models,

ܧሼܿ פ ݔ , ݇ ሽ ൌ ߚ ݔ ߚ ௫ ݇ ் ߚ , for ݊ ൌ 1, . . ,16,
and

ݕሼܧ פ ݔ , ܿ , ݇ ሽ ൌ ߠ ݔ ߠ ௫ ܿ ் ߠ ݇ ் ߠ .
The mediated effect of ݔ on DNA methylation was estimated by ߚ ௫ ் ߠ [START_REF] Vanderweele | Explanation in Causal Inference: Methods for Mediation and Interaction[END_REF].

Inference was done by the parametric bootstrap.

IDOL-adjusted effect

To compute the IDOL-adjusted effect, we estimated proportions of CD4 + and CD8 + T cells, B cells, NK cells, monocytes and neutrophils by the estimateCellCounts2 function in the FlowSorted.Blood.EPIC package with IDOL optimized CpG sites [START_REF] Salas | An optimized library for reference-based deconvolution of wholeblood biospecimens assayed using the Illumina HumanMethylationEPIC BeadArray[END_REF]. For age, sex, smoking status and CMV serostatus, we estimated the IDOL-adjusted effect by adjusting for these estimated 6 proportions in the model specified by Eq. 3, instead of the 16 measured proportions.

Detection of the dispersion of DNA methylation with age

To estimate the change in dispersion of 5mC levels with age, we fit regression models where the residual variance depends on age. Let ݕ be methylation levels on the M value scale for the :th

CpG site and the ݅:th individual. Using a similar notation as above, we estimated the dispersion effect of age by the parameter θ in the model,

ݕ ൌ ߤ ݔ ߚ ௫ ܿ ் ߚ ܵܰܲ ் ߚ ௌே ݂ ఉ ಲ ሺ݁݃ܣ ሻ ܹ݊ܽ݉ ߚ ௐ ݎ݁݇݉ݏݔܧ ߚ ா௫௦ ݎ݁݇݉ܵ ߚ ௌ ܸܯܥ ߚ ெ 1ܥܲ ߚ ଵ 2ܥܲ ߚ ଶ ߝ ,
Eq. 4

where ߝ  ࣨሺ0, ߪ ଶ ሻ, logሼ ߪሽ ൌ ߬ Age ߠ. We devised a hypothesis test for ߠ by a likelihood ratio test comparing that model to a model with ߝ  ࣨሺ0, ߪ ଶ ሻ, logሼ ߪሽ ൌ ߬ in Eq. 4. As a sensitivity analysis, we also fitted a model with ߝ  ࣨሺ0, ߪ ଶ ሻ, logሼ ߪሽ ൌ ߬ ܿ ் ߚ Age ߠ in Eq. 4.

Hypothesis test for ߠ in this case was done by comparing to a model with ߝ  ࣨሺ0, ߪ ଶ ሻ, logሼ ߪሽ ൌ ߬ ܿ ் ߚ in Eq. 4. In this analysis, 77,708 CpG sites showed significant dispersion with age, 10%

of which showed an increase in dispersion.

Cell-type specific methylation changes

Let ݕ be 5mC levels measured at a CpG site on the β value scale for the ݅:th individual. Further, let ݉ be 5mC levels and ܿ measured proportions of the ݈:th cell type and ݔ a vector of variables of interest. Expected average 5mC levels over all cells can be decomposed into,

ݕሼܧ פ ܿ ଵ , … , ܿ , ݔ ሽ ൌ ܿ ܧ ሼ݉ פ ݔ ሽ.
Eq. 5

Now assume that the expected value of ݉ depends linearly on covariates of interest ݔ , ܧሼ݉ פ ݔ ሽ ൌ ߤ ݔ ் ߠ .

Eq. 6

Inserting Eq. 6 in Eq. 5 yields

ݕሼܧ פ ܿ ଵ , … , ܿ , ݔ ሽ ൌ ܿ ߤ ܿ ݔ ் ߠ
Eq. 7

Now, ܿ are proportions, so

ܿ ଵ ൌ 1 െ ܿ . ஷଵ Eq. 8
Inserting Eq. 8 into Eq. 7 and rearranging gives,

ݕሼܧ פ ܿ ଵ , … , ܿ , ݔ ሽ ൌ ߤ ଵ ݔ ் ߠ ଵ ܿ ሺߤ െ ߤ ଵ ሻ ஷଵ ܿ ݔ ் ሺߠ െ ߠ ଵ ሻ ஷଵ ,
Eq. 9

which, by a change of parameters

ߤ ൌ ߤ ଵ , ߚ ൌ ߠ ଵ , ߬ ൌ ߤ െ ߤ ଵ , ߚ ൌ ߠ െ ߠ, gives the interaction model, ݕሼܧ פ ܿ ଵ , … , ܿ , ݔ ሽ ൌ ߤ ݔ ் ߚ ܿ ሺ߬ ሻ ஷଵ ܿ ݔ ் ߚ ஷଵ .
Eq. 10

We can now interpret the parameters of this model. The intercept ߤ corresponds to the intercept level of 5mC levels in a baseline cell type. The main effect terms ߚ and ߬ are the effects of covariates ݔ on 5mC levels in the baseline cell type and the difference in intercept term in the ݈:th cell type compared to the baseline cell type. Finally, the interaction parameter vector of interest ߚ is the difference in the effect of the covariates on 5mC levels in the ݈:th cell type compared to the baseline cell type. Because it is the largest immune cell subset in blood, we used neutrophils as baseline cell type. This model had inflated variance for very small blood subsets. We therefore used the reduced set of 7 immune cells in this analysis. We estimated cell-specific effects on 5mC levels of age, sex and CMV serostatus by using,

ݔ ் ൌ ሺ݁݃ܣ ݔ݁ܵ ܸܯܥ ሻ.
For cell-specific effects of genetic variants, we used a model that additionally included main effect terms for aging, sex and CMV serostatus,

ݕሼܧ פ ܿ ଵ , … , ܿ , ݔ ሽ ൌ ߤ ݔ ் ߚ ௫ ݇ ் ߚ ܿ ሺ߬ ሻ ஷଵ ܿ ݔ ் ߚ ஷଵ ,
Eq. 11

where the vector ݇ contains age (3 DoF spline term), sex and CMV serostatus, and ݔ ൌ ܵܰܲ , the minor allele dosage of the ݊:th SNP associated with 5mC levels at the CpG site for the ݅:th individual. Inference was done by Wald tests with heteroscedasticity-consistent standard errors estimated by the sandwich R package [START_REF] Zeileis | Various Versatile Variances: An Object-Oriented Implementation of Clustered Covariances in R[END_REF]. To test the interaction models in Eq. 10

and Eq. 11, we performed simulations of a system like the one specified in Eq. 5 and Eq. 6, except that the output of Eq. 6 was logit-transformed to ensure that it was a proportion. We used observed 44 cell proportions for the simulations. When we simulated cell-type specific effects of aging in CD8 + T cells and NK cells, drawn from a normal distribution with mean and standard deviation taken from those estimated for a moderate signal in the main effect age EWAS (P adj ≈ 10 -5 ), and zero effects of aging in the other cells, our model correctly detected cell-specific age effects in CD8 + T cells and NK cells, but not in the other cells (Figure S3I).

Detection of gene × environment interactions

We tested whether age, sex, CMV serostatus, smoking status and CRP levels could have a different effect on the methylome depending on genotypes. For the i:th individual, let ܵܰܲ be minor allele dosages of SNPs associated with 5mC levels at the :th CpG site in the M value scale, ݕ and let ܿ be a vector of measured proportions of blood subsets with corresponding parameter vector ߚ ୡ .

Interaction effects for each variable of interest and each associated SNP were estimated for each

CpG site in the model,

Eq. 12

We investigated CRP in a separate model that simply added corresponding log-transformed CRP terms to Eq. 12. Inference was done by Wald tests with heteroscedasticity-consistent standard errors estimated by the sandwich R package [START_REF] Zeileis | Various Versatile Variances: An Object-Oriented Implementation of Clustered Covariances in R[END_REF].

Estimation of proportions of explained 5mC variance

According to our analyses, 5mC levels in the population are mainly associated with local genetic variation, blood cell composition, age, sex, smoking, CMV infection and CRP levels. We grouped these variables into 4 groups: genetic, cell composition, intrinsic and exposures. For a particular CpG site and the ݅:th individual, we collected observations of the minor allele dosage for the most associated local SNP in ݔ , proportions of the 16 major cell types in the vector ݔ , intrinsic factors (sex and natural spline expanded values of age) in the vector ݔ and exposures (smoking status, CMV serostatus and log-transformed CRP levels) in the vector ݔ , with corresponding parameter vectors ߚ , ߚ , ߚ and ߚ . We interpret log-transformed CRP levels as a proxy measure of the exposure of chronic low-grade inflammation. For each group, we define the linear predictor terms:

݂ ൫ݔ ൯ ൌ ݔ ߚ ,
Eq. 13

݂ ሺݔ ሻ ൌ ሺݔ ሻ ் ߚ ,
Eq. 14

݂ ൫ݔ ൯ ൌ ൫ݔ ൯ ் ߚ ,
Eq. 15

݂ ሺݔ ሻ ൌ ሺݔ ሻ ் ߚ .
Eq. 16

These functions vary in complexity, so to get a fair comparison between them, we estimated group effect sizes as the out-of-sample proportion of variance explained by each group predictor. This estimation is done by indexing samples into two disjoint index groups ܫ ଵ and ܫ ଶ , fitting the model on samples from ܫ ଵ , and evaluating the prediction accuracy on samples from ܫ ଶ .

Let ݕ be 5mC levels at a CpG site on the M value scale.

Eq. 18

We can then define the total effect size for group ݊ as the squared correlation between observations and the out-of-sample prediction, ሺܴ Tot ሻ ଶ ൌ cor ቀݕ , ݂ ൫ݔ ൯ቁ ଶ , ݆ א ܫ ଶ .

Eq. 19

For groups other than the cell composition group, we also computed a direct effect. For each group, it was computed as the added out-of-sample proportion of variance explained when adding the group predictor term to that of the cell composition group. The effect was computed for group ݊ by

ሺܴ ሻ ଶ ൌ ሺܴ ା Tot ሻ ଶ െ ሺܴ Tot ሻ ଶ Eq. 20
Where ሺܴ ା Tot ሻ ଶ is the total effect of the predictor including both group ݊ and proportions of cell types

݂ ା ሺݔ , ݔ ሻ ൌ ߤ̂ ሺݔ ሻ ் ߚ ሺݔ ሻ ் ߚ .
Eq. 21

To mitigate the impact of sampling on estimates of total and direct effects, we did four independent repeats of fivefold cross-validation and averaged the effect sizes across all 20 drawn samples. To have an unbiased estimation of the out-of-sample explained variance, we redid a local meQTL mapping on the training set in each iteration of the cross-validation scheme. The algorithm for drawing samples of the total effect is detailed in Algorithm 2.

Algorithm 2: Cross-validation for estimating out-of-sample group total effect size Repeat 4 times:

For ݇ ൌ 1, … , 5
Index a fifth of individuals as ܫ , the others are indexed as Eq. 23

Biological annotations

Information about the position, closest gene and CpG density of each CpG site was obtained from the Illumina EPIC array manifest v.1.0 B4. We retrieved the chromatin state of regions around each CpG site, using the 15 chromatin states inferred with ChromHMM for CD4 + naive T cells by the ROADMAP Epigenomics consortium [START_REF] Roadmap Epigenomics | Integrative analysis of 111 reference human epigenomes[END_REF]. We used CD4 + naive T cells as a reference because it is a large, relatively homogeneous subset of cells that are less differentiated than memory cells. We obtained similar results when using other cell subsets as reference (data not shown). The data was downloaded from the consortium webpage (https://egg2.wustl.edu/roadmap/web_portal/chr_state_learning.html). The transcription factor binding site data used was public CHIP-seq data collected and processed for the 2020 release of the ReMap database [START_REF] Cheneby | ReMap 2020: a database of regulatory regions from an integrative analysis of Human and Arabidopsis DNA-binding sequencing experiments[END_REF], including a total of 1,165 TFs. Binding sites include both direct and indirect binding. Enrichment analyses were performed by creating a simple two-way table for each target set and TF, and then performing a Fisher's exact test. Gene ontology enrichments were computed with the gometh function in the missMethyl R package [START_REF] Phipson | missMethyl: an R package for analyzing data from Illumina's HumanMethylation450 platform[END_REF].

We tested if a set of ݔ local or remote meQTL SNPs is enriched in disease-or trait-associated variants, by sampling at random, among all tested SNPs, 15,000 sets of ݔ SNPs with minor allele frequencies matched to those of meQTL SNPs. For each resampled set, we calculated the proportion of variants either known to be associated with a disease or trait, or in linkage disequilibrium (LD; set here to r 2 > 0.6) with a disease/trait-associated variant (P-value < 5×10 -8 ; 48 EBI-NHGRI Catalog of GWAS hits version e100 r2021-01-1). The enrichment P-value was estimated as the percentage of resamples for which this proportion was larger than that observed in meQTL SNPs. LD was precomputed for all 5,699,237 SNPs with PLINK 1.9 (with arguments 'show-tags all-tag-kb 500-tag-r2 0.6') [START_REF] Chang | Secondgeneration PLINK: rising to the challenge of larger and richer datasets[END_REF].

•

Figure 1 .

 1 Figure 1. Local and Remote Genetic Control of DNA Methylation Variation in Whole Blood (A) Distributions of 5mC levels at 644,517 CpG sites averaged over 958 donors, across 15 chromatin states. (B) Distributions of the variance of 5mC levels at 644,517 CpG sites among 958 donors, across 15 chromatin states. (C) Enrichment in CpG sites associated with local meQTL variants, across 15 chromatin states. (D) Average Genomic Evolutionary Rate Profiling (GERP) scores, across 15 chromatin states. (E) Enrichment of local meQTL variants in disease/trait-associated variants, across 15 chromatin states.

Figure 2 .

 2 Figure 2. Direct and Cell-Composition-Mediated Effects of Aging on the Blood DNA Methylome (A) Total effects against cell-composition-mediated effects of age on 5mC levels (50 year effect). Only CpG sites with a significant total and/or cell-composition-mediated effect are shown. Labels denote genes with strong total or cell-composition-mediated effects of age. (B) Enrichment in CpG sites with significant direct effects of age, across 15 chromatin states. (C) Distributions of significant direct effects of age, across 15 chromatin states. Numbers on the right indicate the number of associated CpG sites and proportion of positive effects. (D) Enrichment of CpG sites with a significant positive, direct effect of age in binding sites for TFs. The 15 most enriched TFs are shown, out of 1,165 tested TFs. (E) Genomic distribution of direct age effects at the BCL6B locus. (F) Number of CpG sites with a significant decreased or increased variance with age.(G) Increased variance of 5mC levels with age at the MAFA locus.

Figure 3 .

 3 Figure 3. Effects of Latent Cytomegalovirus Infection and Smoking on the Blood DNA Methylome (A) Total effects against cell-composition-mediated effects of latent CMV infection on 5mC levels. (B) Enrichment in CpG sites with a significant direct effect of CMV infection, across 15 chromatin states. (C) Distributions of significant direct effects of CMV infection across 15 chromatin states. Numbers on the right indicate the number of associated CpG sites and proportion of positive effects. (D) Enrichment of CpG sites with a significant direct, positive effect of CMV infection in binding sites for TFs. The 15 most enriched TFs are shown, out of 1,165 tested TFs. (E) Genomic distribution of direct effects of CMV infection at the LTBP3 locus. (F) Enrichment of CpG sites with a significant direct, negative effect of CMV infection in binding sites for TFs. The 15 most enriched TFs are shown, out of 1,165 tested TFs.

Figure 4 .

 4 Figure 4. Effects of Gene × Environment Interactions on the Blood DNA Methylome (A) P-value distributions for significant effects of genotype × age, genotype × sex, genotype × exposures and genotype × cell type interactions. The number of significant associations is indicated on the left. Labels denote genes with strong interaction effects.(B) Genotype-dependent effect (rs2837990 variant) of age on 5mC levels at the BACE2 locus.(C) CD8 + T cell-dependent effect of the rs12939435 variant on 5mC levels at the CD300A locus.

Figure 5 .

 5 Figure 5. Best Predictors of the Blood DNA Methylome of Adults (A) Complementary cumulative distribution function of the out-of-sample variance explained by the full model, blood cell composition, genetic factors, intrinsic factors and exposures, for 644,517 CpG sites. (B) Complementary cumulative distribution function of the out-of-sample variance explained by genetic factors, intrinsic factors, exposures and gene × environment (G × E) interactions, when conditioning on blood cell composition, for 644,517 CpG sites. (C) Proportion of the explained out-of-sample variance of 5mC levels for the 20,000 CpG sites with the variance most explained by blood cell composition, genetic factors, intrinsic factors and exposures, respectively.

  study is sponsored by the Institut Pasteur (Pasteur ID-RCB Number: 2012-A00238-35) and was conducted as a single center study without any investigational product. The Milieu Intérieur clinical study was approved by the Comité de Protection des Personnes -Ouest 6 (Committee for the protection of persons) on June 13, 2012 and by the French Agence Nationale de Sécurité du Médicament (ANSM) on June 22, 2012. The samples and data used in this study were formally established as the Milieu Intérieur biocollection (study# NCT03905993), with approvals by the Comité de Protection des Personnes -Sud Méditerranée and the Commission nationale de l'informatique et des libertés (CNIL) on April 11, 2018.

Table 1 . Genetic variants with multiple, remote effects on the blood DNA methylome

 1

	TABLES
	DNA Sequence variant	Chr. Position	Closest gene	eQTL P-value	eQTL direction	#CpG sites	Positive effects (%)	TF	TFBS enrichment [95% CI]
	rs77081633	6	136589425 BCLAF1	-	-	25	48%	BCLAF1	7.4 [3.04, 18]
	rs11850055	14 105754532	BRF1	1.1×10 -22	Positive	11	0%	BRF1	1011.3 [162.91, 4300]
	rs60626639	16 67625797	CTCF	-	-	22	96%	CTCF	16.8 [2.71, 694]
	rs11986122	8	10009949	MSRA	1.3×10 -177	Negative	16	56%	-	-
	rs79755767	12 54698408	NFE2	-	-	14	100%	NFE2	141 [21.16, 5765]
	rs1585215	4	103444474 NFKB1	8.9×10 -45	Positive	14	0%	NFKB1	59.8 [15.79, 337]
	rs12491955	3	101146597	SENP7	4.2×10 -302	Positive	35	9%	KAP1	16.3 [6.66, 48]
	rs10889104	1	59046496 TACSTD2 4.2×10 -302	Positive	10	0%	-	-
	rs3809627	16 30103160	TBX6	3.5×10 -135	Negative	26	31%	-	-
	rs1005278	10 38218748	ZNF25	8.9×10 -171	Positive	10	100%	-	-
	rs10417143	19 22373303	ZNF257	2.8×10 -70	Positive	19	16%	ZNF534	5.5 [0.62, 23]

. DNA sequence variants that affect 5mC levels at more than 10 remote CpG sites can control mRNA levels of a nearby TF or chromatin remodeler, and remotely associated CpG sites can be enriched in binding sites for the corresponding TF, or a TF related to it. Chr. denotes the chromosome where the DNA sequence variant is located. eQTL Pvalue is the P-value of association between the master variant and mRNA levels of a close gene

[START_REF] Võsa | Unraveling the polygenic architecture of complex traits using blood eQTL meta-analysis[END_REF]

. #CpG sites denote the number of CpG sites remotely associated with the DNA sequence variant. TFBS stands for TF binding sites.

Table 2 . Number of CpG sites significantly associated with intrinsic factors and exposures

 2

	pbkrtest R package	(Halekoh and Højsgaard, 2014)	https://cran.r-project.org/web/packages/pbkrtest/index.html
	missMethyl R package	PMID: 26424855	https://doi.org/doi:10.18129/B9.bioc.missMet hyl
	• CONTACT FOR REAGENT AND RESOURCE SHARING

. Out of 141 tested factors, 20 and 16 have a significant total effect (i.e., cell-composition-independent and cell-composition-mediated effects) or a direct effect (i.e., cell-composition-independent effect) on the 5mC levels of more than one CpG site, respectively. "log" denotes log-transformation, CMV: cytomegalovirus, CRP: C-reactive protein, and HDL: high-density lipoprotein.

  To compute the total effect of each group ݊ on CpG methylation, we first fit the predictor function in individuals indexed to ܫ ଵ ,

		݂	ሺݔ	ሻ ൌ ߤ̂ ሺݔ	ሻ	் ߚ	, ݅ א ܫ ଵ	Eq. 17
	with parameters estimated by least squares,
	ሺ ߤ̂ߚ	ሻ ൌ argmin μ,β n ൫ݕ െ ߤ െ ݂ ሺݔ ሻ൯ ଶ	.
				אூ భ

  ܫ ך Select SNP for the predictor ݂ by performing a local meQTL mapping on individuals in ܫ ך For predictor ݂ א ሼ݂ , ݂ , ݂ , ݂ ሽ Estimate ݂ by Eq. 17 and Eq. 18 with ܫ ଵ ൌ ܫ ך Compute ሺܴ Tot ሻ ଶ by Eq. 19 with ܫ ଶ ൌ ܫ Finally, we computed an effect size for interactions between genetic and non-genetic factors. It was computed, similar to Eq. 20, as the added out-of-sample proportion of variance explained by the regression function, ݂ Int ሺ݁݃ܣ , ܹ݊ܽ݉ , ܸܯܥ , ݎ݁݇݉ܵݔܧ , ݎ݁݇݉ܵ , ܴܲܥ ሻ ൌ ߤ ܵܰܲ ߚ ௌே ݁݃ܣ ߚ ܹ݊ܽ݉ ߚ ௐ ܸܯܥ ߚ ெ ݎ݁݇݉ܵݔܧ ߚ ா௫ௌ ݎ݁݇݉ܵ ߚ ௌ ܴܲܥ݃ሺ݈ ሻߚ ோ Main ሺ݁݃ܣ , ܹ݊ܽ݉ , ܸܯܥ , ݎ݁݇݉ܵݔܧ , ݎ݁݇݉ܵ , ܴܲܥ ሻ ൌ ߤ ܵܰܲ ߚ ௌே ݁݃ܣ ߚ ܹ݊ܽ݉ ߚ ௐ ܸܯܥ ߚ ெ ݎ݁݇݉ܵݔܧ ߚ ா௫ௌ ݎ݁݇݉ܵ ߚ ௌ ܴܲܥ݃ሺ݈ ሻߚ ோ .

				Eq. 22
	ܵܰܲ ൫݁݃ܣ ߚ	ௌே ܹ݊ܽ݉ ߚ ௐ ௌே	ܸܯܥ ߚ ெ ௌே
	ݎ݁݇݉ܵݔܧ ߚ ா௫ௌ ௌே	ݎ݁݇݉ܵ ߚ ௌ ௌே	ܴܲܥ݃ሺ݈ ሻߚ ோ ௌே ൯,
	compared to the same regression function without interaction terms:

݂

  DATA AND SOFTWARE AVAILABILITYThe Infinium MethylationEPIC raw intensity data have been deposited at the European Genomephenome Archive (EGA, https://www.ebi.ac.uk/ega/) under accession number EGAS0000XXXXXX. Data access applications are reviewed by a data access committee (DAC) and access is granted if the request is consistent with the consent provided by Milieu Intérieur participants. All association statistics obtained in this study (i.e., local meQTL mapping, the 141 EWAS and interaction models) can be explored and downloaded from the web browser hub05.hosting.pasteur.fr/MIMETH_browser/. All the code supporting the current study has been uploaded to GitHub: https://github.com/JacobBergstedt/mimeth.

ACKNOWLEDGMENTS

We thank Sarah Merrill, Nicole Gladish, Violaine Saint-André, Lucas Husquin and the Milieu Intérieur scientific advisory board for helpful discussions. We acknowledge the help of the HPC Core Facility of Institut Pasteur for this work. This work benefited from support of the French government's program 'Investissement d'Avenir', managed by the Agence Nationale de la Recherche (reference 10-LABX-69-01).

The Milieu Intérieur Consortium ¶ is composed of the following team leaders: Laurent Abel (Hôpital

STAR METHODS

• KEY RESOURCES

SUPPLEMENTAL INFORMATION

The supplemental PDF includes 6 supplemental figures. Supplemental tables are provided as 6 separate excel files:

Table S1. Candidate intrinsic and environmental factors tested for association with the blood DNA methylome of adults.

Table S2. Summary statistics for significant remote-effect meQTLs.

Table S3. Significant enrichments of variable-associated CpG sites in binding sites for transcription factors (TFs).

Table S4. Significant gene ontology enrichments of genes close to variable-associated CpG sites.

Table S5. CpG sites significantly associated with two interacting variables.

Table S6. Proportions of variance explained by intrinsic factors, exposures, cell composition and local SNPs for the 10,000 CpG sites with the most explained variance.