Study population
For this study, we used data from a population-based German birth cohort on the Influence of Life-style factors on Development of the Immune System and Allergies in East and West Germany (LISA). From 1997 to 1999 a total of 3,097 full-term healthy newborns were recruited at four study centers (Munich, Wesel, Leipzig and Bad Honnef). The study was approved by local ethics committees (Bavarian Board of Physicians, Board of Physicians of North-Rhine-Westphalia and Medical Faculty of the University of Leipzig) and written, informed consent was obtained from the parents or legal guardians. In the present study, only data from participants enrolled in the Munich study center with parental consent for genetic analyses at both six and ten years is included (Nmax=240).
Aeroallergen sensitization
Positive aeroallergen sensitization was defined as a specific IgE threshold of >0.35kU/L (at least Radio-Allergo-Sorbent-Test (RAST) class 1), measured for a mix of common aeroallergens (SX1 mix: Dermatophagoides pteronyssinus, cat, dog, rye, timothy grass, Cladosporium herbarum, birch and mugwort). Serum at six, ten and 15 years was analyzed using the CAP-RAST FEIA system (Pharmacia Diagnostics, Freiburg, Germany) according to the manufacturer’s instructions.
Risk factors for aeroallergen sensitization
Genome-wide data in the LISA study was measured using the Affymetrix Chip 5.0 and 6.0 (Thermo Fisher Scientific, USA). More information on genetic data can be found in the supplementary material of Grosche et al. 2021(22). We calculated a PRS for any allergic disease based on the genome-wide significant hits reported in Ferreira et al., 2017 & 2019(15,16). Single nucleotide polymorphisms (SNPs) were extracted for each participant and weighted with the reported effect size. Multiallelic SNPs, highly correlated variants (Linkage disequilibrium R2 > 0.7), those with a low imputation quality (< 0.4) or a minor allele frequency of less than 1% were excluded. Further information on quality control and PRS calculation can be found elsewhere(7).
Information on family history of allergic diseases was collected at birth and defined as a binary factor indicating no family history or at least one biological parent reporting ever experiencing asthma, atopic dermatitis or hay fever.
Maternal smoking during pregnancy was defined as smoking in the second and/or third trimester of pregnancy, with controls defined as either stopped smoking before the second trimester or never smoking. Potential confounders after literature research are sex, age, season at blood withdrawal, cell type proportions, Body-Mass-Index (BMI), socio-economic status (SES), and air pollution, defined as nitrogen dioxide (NO2) at birth address (Table S1).
DNAm data
DNAm was measured for 256 participants from blood clots taken at six and ten years using the Methylation EPIC BeadChip (Illumina, Inc., San Diego, CA). We applied functional normalization(23) and ComBat(24) to normalize the data and remove technical variation. Probes were removed if they were located on the sex chromosomes, had missing values, or failed the detection p-value of 0.01 in more than 1% of samples. Samples were removed if they were outliers, sex mismatches, or did not fulfill the bad-sample threshold of methylated and unmethylated intensities. Cell type proportions were estimated using the EpiDISH package(25). Further information on quality control and data processing can be found elsewhere(7).
Methylation Risk Scores
MRS were calculated for six allergy-related EWAS, namely high IgE(2), aeroallergen sensitization(3), asthma(5), any allergic disease(6) and two on atopy, defined as high total IgE(1) or positive specific IgE as well as a positive skin-prick test(4). Details on the calculation and evaluation of these allergy-related MRS have been published previously(7). In short, we calculated each MRS by weighting the CpG beta-values with the respective effect size identified by the EWAS and transformed to z-scores. The selection of CpG sites was conducted using a pruning and thresholding approach(26). As described previously(7), the MRS that reached the highest prediction accuracy for allergic sensitization at six years of age across all p-value thresholds, was used in the downstream analyses.
Statistical analysis
To evaluate whether changes in DNAm are predictors or consequences of allergic diseases, we tested the following two hypotheses: (H1, DNAm as predictor) The association between exposure (maternal smoking during pregnancy; family history of allergic disease; PRS for any allergies) and allergic sensitization is mediated by prior changes in DNAm (measured by MRS or methylation in individual CpG sites); (H2, DNAm as consequence) The association between exposure (maternal smoking during pregnancy; family history of allergic disease; PRS for any allergies) and changes in DNAm (measured by MRS or methylation in individual CpG sites) is mediated by prior allergic sensitization. In our main analyses, mediators were measured at six years and outcomes at ten years, both for hypothesis (H1) and (H2). In addition, we conducted a secondary analysis for hypothesis (H1), in which mediators (DNAm) were measured at ten years and outcome (aeroallergen sensitization) at 15 years (Fig. 1 & S1).
Mediation analyses rely on the following three assumptions(27): 1) no exposure-mediator confounding, 2) no mediator-outcome confounding and 3) no exposure-outcome confounding. To fulfill these assumptions to the best of our knowledge, we constructed directed acyclic graphs (DAGs) to visualize each of these paths using dagitty(28) (Figures S2-S9). A minimal sufficient adjustment set was identified for each pathway via the tracing of association directions and elimination of any potential confounders already associated with a precursory confounder. Exposure-mediator models were adjusted for SES (Exposure: maternal smoking during pregnancy), SES and NO2 exposure at birth (family history) and sex (PRS) for both hypotheses. Mediator-Outcome models were adjusted for all potential confounders according to the DAGs (Figures S2-S9). A detailed description of the definition and assessment of these covariates is provided in Table S1.
Associations with continuous outcomes (MRS or DNAm in individual CpG sites) were analyzed using linear regression and associations with binary outcomes (allergic sensitization) were analyzed using logistic regression.
Causal mediation analysis of MRS
Causal mediation analysis, using the R package mediation(29), was applied to test the two hypotheses (H1) and (H2) for allergy-related MRS. Results were adjusted for multiple testing using the Benjamini-Hochberg procedure(30) for false-discovery rate (FDR) together within each H1 and H2.
High-dimensional mediation analysis of individual CpGs
High-dimensional mediation analyses (HMA) were used to test the two hypotheses for individual CpGs. H1 was tested using the Divide-Aggregate Composite-Null test (DACT), HIMA, and gene-based HMA (gHMA). H2 was tested using only DACT, because HIMA and gHMA are only applicable for high-dimensional mediators but not for high-dimensional outcomes.
1) Previous knowledge + Divide-Aggregate Composite-Null test (DACT)
Based on previously published EWAS of total IgE(1,2), aeroallergen sensitization(3,4), childhood asthma(5) and any allergic disease(6) we used existing knowledge on allergy-relevant CpGs to reduce the multiple testing burden. Of the 1673 previously reported CpGs, 1501 were available in the LISA cohort and 583 were significantly validated (False discovery rate <= 0.05) in the LISA cohort at six years(7), which were further taken as testing-set of potential mediators. We used DACT for the composite null hypothesis of no mediation effect as suggested by Liu et al.(31) to improve the multiple testing burden. In short, DACT takes the p-values from the exposure-mediator and the mediator-outcome model to compute a new joint list of p-values, which will be used to determine significance (p-value < 0.05). This is done by aggregating the weighted p-values of the three possible null-hypotheses leading to no mediation effect and calibrating this using Efron’s empirical null framework(32).
2) HIMA
Whereas the previous approach relied on existing knowledge as a baseline selection of mediators, HIMA as proposed by Zhang et al.(33) uses a three step procedure to identify significant CpGs throughout the whole epigenome. First, the top CpGs with the largest effect sizes (beta of standardized inputs) for the response variable are identified using sure independence screening (SIS)(34). The total number of top hits (N) varies per model and is calculated by N = 2*n/log(n), with n being the input sample size. To capture relevant CpGs with our smaller sample size, we applied a looser threshold than the original publication. In the final steps, HIMA estimates the mediation effect using minimax concave penalty and performs joint significance tests.
3) Gene-based HMA (gHMA)
We further applied gene-based high-dimensional mediation analysis (gHMA) as proposed by Fang et al.(35). The idea behind this approach is that not single CpGs but genes act as biological units and should therefore be analyzed together. The functions further provide different modeling options for linear or non-linear relationships and an omnibus-test to combine both, which outperformed the single models in their simulation study. First, we annotated every CpG to their nearest gene within 20,000 base pairs as done previously(36), resulting in 40,916 different genes. We then applied gHMA to each of these 40,916 genes, each covering between one and 1758 CpGs, performing the linear, non-linear and omnibus-test for significance. We used differing kernel-thresholds of 0.7, 0.8 and 0.9 as values for explained variance by the kernel principal components. Results of the omnibus-test were corrected using the Benjamini-Hochberg procedure(30).
Validation of CpG sites using causal mediation analysis
All significant CpG sites identified with the three described methods above are followed up using a causal mediation analysis to determine the direct, indirect, and total effects as well as the proportion mediated. Multiple testing correction followed the one applied for the MRS evaluation by calculating the FDR for all H1 CpGs together, the same correction was applied for H2 CpGs. Models and adjustment are the same as for MRS analyses and we conducted a sensitivity analysis additionally adjusting all exposure-mediator associations for estimated cell types.
All analyses were performed in R V.4.1.2(37).