Epigenome-wide analysis identifies DNA methylation mediators of treatment-related cardiometabolic risk in survivors of childhood cancer

Domain: neuroscience
Published: 2026-01-23
Editorial outcome: major_revision

Abstract
Childhood cancer survivors face increased cardiometabolic risks from cancer treatment exposures, yet mechanisms remain unclear. Here, epigenome-wide analysis identifies 1893 DNA methylation (DNAm) sites in peripheral-blood-mononuclear-cells (PBMCs) associated with at least one cardiometabolic risk factor (CMRF), including obesity (n = 1720), abnormal glucose (n = 201), hypertriglyceridemia (n = 145), hypercholesterolemia (n = 38) and hypertension (n = 34) in 2938 survivors from the St. Jude Lifetime Cohort. A core set of five DNAm sites near CPT1A and LMNA is associated with all CMRFs. Mediation analyses identify 24 sites mediating associations between treatments and CMRFs, implicating inflammatory and metabolic pathways. Notably, cg20370568, a cis-expression quantitative trait methylation site for ANTXR2, mediates 20% of the effect of body-trunk-radiotherapy on abnormal glucose. These findings suggest that prior genotoxic cancer treatments may become biologically embedded through DNAm variations that could contribute to cardiometabolic dysfunction and highlight candidate biomarkers for refining risk stratification and guiding intervention strategies in survivorship care.

Introduction
Cardiovascular disease (CVD)-related death is the leading cause of non-cancer late mortality (i.e., 5 or more years after diagnosis of primary cancer) in childhood cancer survivors, accounting for ~26% of deaths within 45 years post diagnosis1. Survivors face a substantially higher burden of CVD than individuals without a history of childhood cancer, with an earlier age of onset and a well-established association with prior exposures to genotoxic cancer treatments, including anthracycline chemotherapy, and chest-directed radiotherapy2. Importantly, cardiometabolic risk factors (CMRFs)—including abnormal glucose, obesity, hypertension, hypertriglyceridemia, and hypercholesterolemia—potentiate risks of future CVD and late mortality3,4. While treatment-related CMRFs are highly prevalent, the underlying biological mechanisms linking cancer treatments to cardiometabolic dysfunctions remain poorly understood. Elucidation of these molecular mechanisms is critical to advance the precision medicine paradigm of personalized risk stratifications and targeted interventions, ultimately improving long-term survivorship outcomes.

DNA methylation (DNAm) is a dynamic and environmentally responsive epigenetic mechanism that integrates signals from internal and external stressors, including genotoxic cancer therapies such as chemotherapy and radiotherapy (RT). Although inter-individual variability in DNAm is common, individuals who received specific cancer therapies can exhibit persistent shifts relative to those without exposures5. Because these DNAm changes may endure for years, they represent a plausible link in the causal pathway between prior cancer treatment and subsequent long-term adverse health outcomes. Our previous study in the well-characterized St. Jude Lifetime Cohort (SJLIFE) identified 935 5′-cytosine-phosphate-guanine-3′ (CpG) sites exhibiting persistent DNAm variations associated with various cancer treatments5. Moreover, a subset of these CpGs mediated the association between specific cancer treatments and CMRFs, suggesting that DNAm may serve as a molecular intermediary linking cancer therapy exposures to long-term cardiometabolic dysfunctions. Notably, DNAm at four independent CpGs (cg06963130, cg21922478, cg22976567, cg07403981) mediated 70.3% of the association between abdominal-RT and hypercholesterolemia. Two of these CpGs reside near genes with well-established roles in lipid metabolism: cg21922478 (ITGA1) and cg22976567 (LMNA), suggesting the biological plausibility of treatment-induced DNAm variations contributing to CMRFs.

Beyond its role in the biological embedding of cancer treatment effects, DNAm is a putative regulator of gene expression, influencing metabolic pathways implicated in cardiometabolic dysfunction6. DNAm-regulated genes are involved in fatty acid synthesis, lipid transport, and storage7,8, as well as inflammatory pathways9,10 that contribute to insulin resistance, abnormal glucose, and hypertension in the general population. Among long-term survivors of childhood cancer, we previously conducted another outcome-centric epigenome-wide association study (EWAS) that complemented our exposure-centric5 EWAS approach, identifying CpGs associated with circulating lipid levels11. This analysis identified 106 lipid-associated CpGs, reinforcing the functional relevance of DNAm in metabolic dysfunction and highlighting its potential as a biomarker for identifying survivors with existing lipid abnormalities or those at higher risk for developing metabolic disorders.

Despite these insights, the biological mechanisms linking cancer treatment exposures to cardiometabolic toxicities remain incompletely understood. Survivors exhibit substantial variability in the timing of onset, clinical presentation, and long-term outcomes of CMRFs, highlighting the need for mechanistic studies that move beyond the statistical association. While our prior work demonstrated that DNAm may mediate the effects of specific treatments on select CMRFs, those analyses were designed around treatment-associated CpGs, assessed a limited subset of treatment-CMRF pairs, and did not examine transcriptional consequences.

Here, we adopt an outcome-first design to assess whether peripheral-blood DNAm mediates associations between prior treatments and CMRFs (Fig. 1). We first map CpGs associated with each CMRF via an epigenome-wide analysis, then evaluate their links to treatment exposures, and estimate mediation for treatment-DNAm-CMRF trios. Additionally, we consolidate uncorrelated CpGs into composite DNAm scores for composite mediation analysis, and we integrate cis-eQTM analyses to connect mediator CpGs to nearby gene expression. By integrating EWAS, causal mediation analysis, and cis-eQTM mapping, we identify CpGs that are not only statistically associated with CMRFs but also biologically informative. This comprehensive framework advances our understanding of DNAm as one type of molecular intermediary between cancer treatments and late effects and highlight candidate biomarkers with potential utility for refining risk stratification and guiding future strategies to mitigate CMRFs and further CVD risk as well as premature mortality in survivors of childhood cancer.

Results
Characteristics of the study population
The study included 2938 childhood cancer survivors (53% male) from the SJLIFE cohort (Fig. 1a). The median age at primary cancer diagnosis was 6.1 (interquartile range [IQR], 2.8–12.6) years, with blood sampling for DNAm conducted at a median age of 30.1 (IQR, 21.6–38) years, and the most recent clinical follow-up occurring at a median age of 33.4 (IQR, 25.4–41.8) years (Table 1). Survivors represented a diverse range of primary cancer diagnoses, including leukemia (34.0%), lymphoma (18.6%), embryonal tumors (13.4%), central nervous system (CNS) tumors (12.7%), sarcoma (12.3%), and other malignancies (9.0%).

Survivors underwent multimodality therapy. Most received vinca alkaloids (71.4%), anthracyclines (57.0%), and alkylating agents (56.7%), while antimetabolites (48.3%), corticosteroids (45.5%), and epipodophyllotoxins (35.2%) were also common. Given the strong correlations among chest, abdominal, and pelvic RT exposures (Phi coefficient > 0.65) and their likely shared systemic effects on peripheral-blood-mononuclear-cells (PBMCs) derived DNAm, we defined a composite body-trunk-RT variable as irradiation to any of these body regions. Overall, 26.7% of survivors received body-trunk-RT and 26.5% received brain-RT. Treatment exposures were frequently concurrent, with strong correlations between corticosteroids and antimetabolites (Phi = 0.77), asparaginase enzymes and antimetabolites (0.71), corticosteroids and asparaginase enzymes (0.69), and corticosteroids and vinca alkaloids (0.57; see Supplementary Fig. 1).

CMRFs were highly prevalent, affecting 58% of survivors. Prevalent CMRFs were defined using the highest Common Terminology Criteria for Adverse Events (CTCAE) grade recorded before or after the blood draw for DNAm. Obesity was most common (41.0%) and defined as CTCAE grade ≥ 3, (BMI ≥ 30 kg/m² in adults or ≥ 95th percentile for their age and sex in children). Hypertension (29.2%), hypertriglyceridemia (20.4%), hypercholesterolemia (16.4%), and abnormal glucose (10.4%) were defined as CTCAE grade ≥ 2, reflecting moderate or more severediseases. Despite the high prevalence of CMRFs, co-occurrence was modest (pairwise Phi < 0.35), except for hypertriglyceridemia and hypercholesterolemia (Phi = 0.78; Supplementary Fig. 2). Survivors had a median of two CMRFs, with 1242 exhibiting none, 827 one, 406 two, 235 three, 178 four, and 78 all five CMRFs.

To support causal inference in downstream mediation analyses, we separately defined incident cases as those developing or worsening after DNAm sampling. Incidence was lower than prevalence given short duration of follow-up since blood draw for DNAm; obesity (8.7%), hypertension (7%), hypertriglyceridemia (2.7%), hypercholesterolemia (1.7%), and abnormal glucose (2.5%). These incident cases provided a temporal alignment for investigating DNAm mediation of treatment-related cardiometabolic risk.

DNAm signatures of CMRFs reveal distinct and shared epigenetic associations
To delineate the epigenetic landscape underlying post-treatment cardiometabolic risks, we performed separate EWAS analyses for each of the five CMRFs using prevalent case status (Fig. 1b). Across 752,683 CpGs, we identified 1894 unique CpGs significantly associated with one or more CMRF (p < 9 × 10−8), yielding a total of 1912 CpG-CMRF associations. Obesity, the most prevalent CMRF, had the largest number of associated CpGs (n = 1720), followed by abnormal glucose (n = 201), and hypertriglyceridemia (n = 145), whereas hypercholesterolemia (n = 38) and hypertension (n = 34) yielded fewer associations (Fig. 2a). These CpGs showed distinct distributions of effect sizes by CMRF (Fig. 2b) and were enriched in epigenetically dynamic regions, particularly open sea regions and 5′ untranslated regions (UTRs), while being significantly less common in CpG islands and near transcription start sites (TSSs) (Supplementary Fig. 3).

Most CMRF-associated CpGs (91%; 1715/1893) were specific to a single condition, while 9% overlapped across two or more CMRFs (Fig. 2c). The largest overlaps occurred between obesity and abnormal glucose (n = 74 CpGs), obesity and hypertriglyceridemia (n = 27), and obesity, abnormal glucose, and hypertriglyceridemia (n = 16). Genes mapped to shared CpGs were enriched in metabolic and inflammatory pathways, including the FTO (fat mass and obesity associated) Obesity Variant Mechanism12, JAK-STAT13, and IL-2/STAT514 signaling. Additionally, hypertriglyceridemia and obesity showed co-enrichment in fatty acid beta-oxidation and integrin signaling, suggesting a unified metabolic–inflammatory axis (Supplementary Data 1).

Strikingly, five CpGs (cg00574958, cg03725309, cg05325763, cg17058475, cg22976567) were associated with all five CMRFs, suggesting potential regulatory hubs. Three (cg00574958, cg05325763, cg17058475) mapped to CPT1A, a gene central to fatty acid metabolism15,16, while cg22976567 mapped to LMNA, a gene implicated in heart disease and dyslipidemia17,18. These shared loci suggest an interconnection across key biological pathways underlying cardiometabolic dysfunction in survivors. The presence of these shared CpGs highlight possible convergence points across metabolic and cardiovascular pathways, underscoring the biological interconnectedness of post-treatment cardiometabolic dysfunction in survivors.

To separate total from lifestyle-independent associations, we fit EWAS with prespecified additions of smoking and physical activity, then body mass index (BMI) (non-obesity outcomes), and diet quality (Health Eating Index [HEI]) in a HEI-complete subset of survivors. Across conditions, smoking/physical activity introduced little changes, whereas most loss of statistical significance occurred when BMI was added. Retention of significant results from the baseline model ranged from high for obesity (no BMI; 81% retained) to low for hypertension (8%), with intermediate retention for hypertriglyceridemia (48%), hypercholesterolemia (37%), and abnormal glucose (16%; Supplementary Fig. 4). Effect sizes for retained loci were highly concordant across models (Spearman ρ = 0.75–0.96), indicating mainly uniform shrinkage rather than sign reversals (Supplementary Fig. 5); median absolute attenuation of the regression β coefficient was modest (2–11%; up to 18% for the single hypertension locus; Supplementary Fig. 6). In the HEI-complete subset, adding HEI had negligible additional impact (ρ = 0.99 vs. the BMI+smoking+physical activity model; Supplementary Fig. 7). Notably, multi-phenotype “hub” CpGs, including three in CPT1A, remained robust for obesity and hypertriglyceridemia but were particularly weakened for hypertension after BMI adjustment (Supplementary Fig. 8).

Relevance of CMRF-associated CpGs in survivors relative to the general population
To assess whether CMRF-associated CpGs identified in survivors reflect broader cardiometabolic biology, we cross-referenced them with the EWAS Catalog19, which aggregates findings from primarily non-cancer populations. We observed significant enrichment for CpGs associated with C-reactive protein (CRP)—a biomarker of systemic inflammation20 (one-sided hypergeometric test, Benjamini–Hochberg [BH] p < 9.29 × 10−81, Supplementary Data 2). CRP-associated CpGs were highly represented among those linked to hypertension (90.9%), hypercholesterolemia (81.6%), hypertriglyceridemia (81.4%), and abnormal glucose (85.1%), yet notably lower in obesity (57.3%).

While many CMRF-associated CpGs overlapped with loci reported in general population EWAS, a substantial subset, particularly among obesity-associated CpGs, appeared unique to survivors. Obesity had the highest proportion of CpGs not reported for any trait in the EWAS Catalog (14%), compared to abnormal glucose and hypercholesterolemia (8%), hypertriglyceridemia (6%), and hypertension (3%), suggesting survivor-specific epigenetic signatures. Among the 1720 obesity-associated CpGs, 155 (9%) had been linked to BMI, waist circumference, or obesity-related traits (Supplementary Data 3). While some CpGs align with established obesity phenotypes, a substantial fraction appears unique to this cohort, potentially reflecting survivor-specific epigenetic signatures (Supplementary Data 4). These findings underscore both shared and unique components of cardiometabolic risk in survivors.

Cancer treatment exposure is associated with widespread DNAm variations
We next examined the extent to which cancer treatments were associated with CMRF-related CpGs, focusing on seven chemotherapy agents (asparaginase enzymes, alkylating agents, antimetabolites, anthracyclines, corticosteroids, epipodophyllotoxins, and vinca alkaloids) and two site-specific radiotherapies (body-trunk-RT and brain-RT; Fig. 1c). We did not examine associations between primary cancer diagnosis and DNAm, as cancer diagnoses are highly collinear with treatment exposures in this cohort (See “Methods”; Supplementary Data 5). Among 1893 CpGs associated with prevalent CMRFs, 64% (1218 CpGs) were also associated with at least one treatment (BH-adjusted p < 0.05), resulting in 2574 significant treatment-CpG associations (Fig. 3).

RT exposures had the strongest epigenetic signatures, with body-trunk-RT (n = 603) and brain-RT (n = 499), accounting for most associations. Among chemotherapies, alkylating agents (n = 346) and asparaginase enzymes (n = 319) had the most CpG associations, while vinca alkaloids (n = 67) had the fewest. Over half (55%) of treatment-associated CpGs were linked to multiple treatments; 12 CpGs were associated with seven or more treatments. One CpG, cg04880464, was associated with eight therapies (all but vinca alkaloids) and also associated with prevalent abnormal glucose, hypertriglyceridemia, and obesity; this CpG has previously been linked to CRP levels in the EWAS Catalog.

Overlap between treatment- and CMRF-associated CpGs varied by condition (Supplementary Fig. 9). Body-trunk-RT showed the highest median overlap across CMRFs (25.1%), followed by alkylating agents (16.8%), while vinca alkaloids showed minimal overlap (1.5%). Of note, obesity-associated CpGs were less frequently associated with treatments; 23.4% of obesity-associated CpGs showed no treatment associations, compared to 7.5% for hypertension, 4.3% for hypertriglyceridemia, 2.4% for abnormal glucose, and 1.9% for hypercholesterolemia, indicating that a subset of obesity-related DNAm variations may be independent of treatment exposures.

DNAm partially mediates the effects of cancer treatments on CMRFs
To investigate whether DNAm acts as an intermediary between cancer treatments and CMRFs, we examined the three-part axis of treatment, DNAm, and CMRF (Fig. 1d). To preserve the temporal sequence necessary for causal inference, we limited this analysis to incident CMRF events arising after DNAm blood sampling (Table 1). CpGs were selected for mediation testing if they were significantly associated with both a treatment and a CMRF in incident-only models; some prevalent associations did not replicate, likely reflecting reduced statistical power (Supplementary Fig. 10). We applied a counterfactual-based mediation framework21 to each of eleven treatment-CMRF pairs.

Following prior methodology5, we defined “strong mediators” as CpGs with BH-adjusted average causal mediation effect (ACME) p < 0.05 that explained more than 10% of the total effect. Using this criterion, we identified 24 CpGs mediating five treatment-CMRF associations (Fig. 4a and Supplementary Data 6). Mediation proportions ranged up to 24%, suggesting that DNAm variation may partially explain the cardiometabolic toxicity associated with cancer treatment exposure. Several CpGs mediated more than one treatment-CMRF association, suggesting shared epigenetic pathways. Five CpGs (cg00574958, cg05325763, cg09737197, cg17058475, cg26731839) mediated the effects of both brain-RT and body-trunk-RT on abnormal glucose. Another CpG (cg09921385) mediated the effect of both body-trunk-RT and anthracyclines on hypertriglyceridemia, highlighting potential convergence of chemotherapy and radiotherapy pathways in driving metabolic dysregulation.

Of particular note, cg20370568 mediated 20% of the effect of body-trunk-RT on abnormal glucose (ACME p = 1.0 × 10−3). This CpG is located near ANTXR2, a gene implicated in extracellular matrix (ECM) remodeling, inflammatory regulation, and vascular homeostasis22, supporting its potential relevance in therapy-related metabolic dysregulation.

Core DNAm mediators preserved despite lifestyle and outcome adjustments
Restricting the mediation analysis to incident-only CMRFs (i.e., those developed strictly after the DNAm blood draw) produced modest reductions in proportion mediated (<5% attenuation) for most condition-treatment pairs, with the largest decrease for abnormal glucose (4.2%) and a small increase for hypercholesterolemia (+2.0%; Supplementary Fig. 11). By treatment, attenuation was 4.9% for brain-RT, 1.8% for anthracyclines, and 1.6% for body-trunk-RT. ACME estimates were highly consistent across definitions (Spearman ρ = 0.89, p = 3.3 × 10−7, Supplementary Fig. 12), and effect directions were preserved for all 30 trios. Retention was full for abnormal glucose with body-trunk-RT (7/7 retained) and hypercholesterolemia with brain-RT (1/1 retained), partial for hypertriglyceridemia with anthracyclines (6/8), and attenuated for abnormal glucose with brain-RT (5/12). All five CpGs mediating radiotherapy effects on abnormal glucose in the primary analysis (cg00574958, cg05325763, cg09737197, cg17058475, cg26731839) remained significant with at least one exposure, and cg20370568 near ANTXR2 retained robust mediation (20.6%, p < 1 × 10−3, Supplementary Data 6).

We next evaluated robustness to lifestyle covariates adjustment, sequentially adding smoking/physical activity and then BMI. For abnormal glucose and hypercholesterolemia, proportion mediated was largely unchanged after smoking/physical activity (<0.2% attenuation) but decreased after BMI (2.7 and 3.5%, respectively, Supplementary Fig. 13), with similar trends for ACME (Supplementary Fig. 14). In contrast, hypertriglyceridemia showed greater reduction at the smoking/physical activity step (4.78%) with little further change after BMI (0.44%). ACME rankings were strongly preserved across models (Spearman ρ = 0.92, p < 1 × 10−6; ρ = 0.72, p < 1 × 10−5 after sequential adjustment), and effect directions were unchanged. Retention patterns mirrored these shifts, with some loss for brain-RT mediators after BMI (abnormal glucose 3/12 retained; hypercholesterolemia 0/1 retained; Supplementary Fig. 15), while trunk-RT mediators were fully retained. Notably, CpGs at CPT1A and ANTXR2 (cg20370568) remained significant, including cg20370568 mediating 15% after BMI adjustment (5% attenuation, p = 0.006; Supplementary Data 6).

Composite methylation scores capture cumulative mediation effects
Among the five treatment-CMRF pairs with strong mediators, four involved multiple CpGs that met mediation criteria, with a median of seven mediators per pair (range, 2–12), representing 23 unique CpGs in total (Supplementary Data 6). Using the correlation-pruned mediator sets (defined in “Methods”), we constructed a composite DNAm score for each treatment-CMRF pair and re-estimated mediation. The composite analyses yielded a single DNAm signature-level indirect effect (Fig. 1e).

Two treatment-CMRF pairs—body-trunk-RT and brain-RT with abnormal glucose—each had four independent mediators based on pairwise correlation. Mediation analysis using the resulting composite scores revealed significant effects: 12% of the association between brain-RT and abnormal glucose (ACME p < 1 × 10−3; Fig. 4b) and 22% for body-trunk-RT and abnormal glucose (ACME p < 1 × 10−3; Fig. 4c) were mediated by DNAm.

Sensitivity analyses using the minimum (rather than the maximum) ACME estimate to select representative CpGs yielded similar results, with mediation proportions of 13% for brain-RT and 20% for body-trunk-RT, both significant at ACME p < 1 × 10−3. These results support the robustness of the composite methylation score approach in capturing multi-CpG mediation effects and provide further evidence that DNAm may play a role in the biological embedding of treatment exposures contributing to CMRFs.

DNAm associated with expression of genes in CMRF-related pathways
To examine whether CMRF-associated CpGs influence gene expression, we conducted a cis-eQTM analysis in 199 survivors with matched blood-based transcriptomic data (Fig. 1f). Using linear regression, we tested CpG-gene pairs within ±1 megabase (Mb) of a gene’s TSS, identifying 79 CpGs (4.2%) significantly associated with expression of 76 genes (eGenes, BH-adjusted p < 0.05), yielding 83 unique eQTM pairs. Most eGenes (91%) were associated with a single CpG (eCpG). Mapping these eQTMs to their respective CMRFs revealed 100 CMRF-eQTM associations, predominantly from obesity (n = 72 eQTMs), followed by abnormal glucose (n = 18), hypertriglyceridemia (n = 5), and hypercholesterolemia (n = 3), and hypertension (n = 2). Many eCpGs had prior associations with metabolic and inflammatory traits in the EWAS Catalog (Supplementary Data 7), suggesting functional relevance and possible causal contribution to CMRFs.

We also identified 144 treatment-eQTM associations (Fig. 5a), with body-trunk-RT (n = 31) and brain-RT (n = 25) showing the most associations, followed by alkylating agents (n = 24), asparaginase enzymes (n = 21), and anthracyclines (n = 14). Several CpGs were associated with multiple treatments. Notably, the eQTM cg06096184-LRIG1, was associated with seven treatments, implicating LRIG1—a negative regulator of epidermal growth factor receptor (EGFR) signaling23—in shared epigenetic response to therapy-induced cellular stress. This CpG has been associated with smoking, aging, COPD, CRP levels, type 2 diabetes, and nitrogen dioxide exposure in the EWAS Catalog.

Finally, cg20370568, identified as a mediator of body-trunk-RT on abnormal glucose, was significantly associated with expression of ANTXR2 (β = −0.179, BH p = 0.048; Fig. 5b, c). The inverse DNAm-expression association was substantially stronger in individuals exposed to body-trunk-RT (Pearson r = −0.48, p = 2.49 × 10−5) compared to unexposed individuals (Pearson r = −0.11, p = 0.21), and a formal interaction test confirmed a significant DNAm-by-RT interaction on gene expression (p = 0.0035; Fig. 5d). This CpG has also been associated with CRP levels, type 2 diabetes, and COPD in the EWAS Catalog (Supplementary Data 7). Its dual identification as both a treatment-related mediator and an eQTM highlights it as a compelling candidate for future mechanistic investigations of therapy-related cardiometabolic risk.

Methods
Participants
Given the SJLIFE cohort’s demographic composition (>80% are individuals of European ancestry), the current availability of DNAm data, and the need to minimize population stratification, analyses were restricted to survivors of European ancestry. Among 3059 survivors with DNAm data, 3044 had full genotype data (for genetic ancestry), and 2966 had CMRFs clinically ascertained. We further excluded 28 participants without radiotherapy details, yielding a final sample of 2938. Potential selection bias arising from these exclusions was evaluated as described in the “Methods: Statistical Analysis”; no systematic associations between exclusion and measured pre-exposure characteristics were detected (Supplementary Data 8).

CMRFs
CMRFs were ascertained within the SJLIFE cohort using a standardized protocol that integrates retrospective medical record abstraction, patient- or proxy-reported health questionnaires, and comprehensive in-person clinical evaluations conducted at baseline and approximately every 5 years thereafter. These evaluations included physical examinations, laboratory testing, and diagnostic assessments, supplemented by external medical records to verify interim medical events. This ensured consistent ascertainment across baseline and follow-up.

Outcomes were graded using a modified version of the National Cancer Institute CTCAE (version 4.03)38, with grades ranging from 0 (no problem) to 4 (life-threatening). A multidisciplinary clinical review panel developed a standardized rubric to harmonize definitions and ensure clinical relevance. As previously described39, grading thresholds were adapted from established guidelines, including the National High Blood Pressure Education Program (hypertension) and Centers for Disease Control and Prevention definitions for BMI percentiles (obesity), with pediatric-specific thresholds applied where appropriate. Abnormal glucose metabolism, hypertriglyceridemia, and hypercholesterolemia were defined using CTCAE laboratory and grading criteria (see Supplementary Methods).

To account for potential bias from differential follow-up time, we defined follow-up time as the interval from blood draw for DNAm to each participant’s last clinical assessment and included this variable as a covariate in all models (see “Statistical analysis: Mediation analysis subsection”). Prior work in SJLIFE has shown that survivors who return for follow-up are largely comparable to those who do not across demographic and clinical characteristics, with only modest differences in diagnostic group distribution35. Attrition in follow-up based on the baseline cardiometabolic risk was evaluated and incorporated as detailed in the “Methods: Statistical Analysis” section, comparing follow-up by incident status and baseline CMRF and adjusting all models for follow-up duration.

For analysis, each CMRF was coded as a binary variable, indicating the presence or absence of the condition. The highest recorded CTCAE grade was used to determine case status for each participant. A CTCAE grade of ≥2 indicated clinically significant abnormal glucose, hypercholesterolemia, hypertriglyceridemia, and hypertension. A more stringent threshold of CTCAE grade ≥3 was applied for obesity (consistent with the clinical definition of BMI ≥ 30 kg/m2 in adults and BMI ≥ 95th percentile for their age and sex in children). These thresholds were selected based on clinical relevance, ensuring that cases captured reflect clinically significant disease burden (see Supplementary Methods for grading criteria).

For EWAS analyses, CMRFs were classified as prevalent events, meaning they were considered present if they had occurred at any time before or after blood collection. This approach maximized statistical power and allowed for comprehensive assessment of DNAm correlates of disease status. For mediation analyses, in order to maintain temporal alignment with causal mediation assumptions, we restricted our analysis to incident CMRF events, meaning they were considered present only if they developed or worsened after the blood collection (i.e., an increase in CTCAE grade following baseline reaching grade ≥ 2; Table 1).

Cancer treatments
Cancer treatment histories were abstracted from medical records, capturing both chemotherapy and radiotherapy (RT) exposures with potential relevance to cardiometabolic risk. Chemotherapeutic agents included alkylating agents (classic), anthracyclines, antimetabolites, asparaginase enzymes, corticosteroids, epipodophyllotoxins, and vinca alkaloids. RT exposures were classified based on the primary exposed body region, including brain-, chest-, abdominal-, and pelvic-RT. Due to strong correlations among abdominal-, pelvic-, and chest-RT (Phi coefficient > 0.65), and their systemic effects on peripheral blood-derived DNAm, we generated a composite variable body-trunk-RT to indicate exposure to any of these three regions. This approach reduced collinearity while preserving biological relevance in analyses of RT-associated DNAm variations. All chemotherapy classes and RT exposures were defined as binary variables indicating whether a given treatment was administered within the first 5 years following cancer diagnosis, consistent with the categorizations shown in Table 1. Treatments administered beyond this period were uncommon (<6% of survivors), indicating that nearly all exposures occurred during the initial treatment phase. A sensitivity analysis excluding a small proportion of survivors who received any treatment beyond 5 years post diagnosis yielded results consistent with the primary analysis. Each of the eleven treatment classes were considered separately for all statistical analysss.

Covariates
To account for potential confounders, all analyses were adjusted for key demographic, genetic, and leukocyte composition variables. Age at blood draw was included for models where DNAm was the dependent variable, while age at the most recent clinical visit was used for models assessing CMRFs to ensure alignment with the most updated health status and minimize biases due to variation in follow-up time. Sex was included as a binary covariate. To adjust for population stratification, we included the top ten principal components (PCs) derived from principal component analysis of common genetic variants in whole-genome sequencing (WGS) data. Global DNAm variability was accounted for by including the top four DNAm PCs, selected based on a scree plot and the proportion of variance explained, ensuring that key axes of variation were captured while minimizing overfitting. Leukocyte subtype compositions were estimated from the DNAm data using the Houseman method40, including proportions for B cells, CD4+ T cells, CD8+ T cells, monocytes, neutrophils, and natural killer (NK) cells. Additionally, experimental batch effects were assessed and adjusted in models incorporating gene expression data.

Lifestyle variables (smoking, physical activity, BMI, and diet quality by the Healthy Eating Index [HEI-2010]) were incorporated in prespecified sensitivity models. Smoking was ascertained by survey and coded as ever vs. never (≥100 lifetime cigarettes). Physical activity was coded as ≥150 vs. <150 min/week of at least moderate activity. BMI was modeled categorically; adults (underweight < 18.5 kg/m2; normal 18.5–24.9 kg/m2; overweight 25.0–29.9 kg/m2; obese ≥ 30 kg/m2) and youth (< 20 years old) using CDC growth-chart41 percentiles (underweight < 5th; normal weight 5–85th; overweight 85–95th; obese ≥ 95th). HEI-2010 was derived from a nutrition questionnaire using standard scoring which assesses adherence to the 2010 Dietary Guidelines for Americans42. To reflect potential confounding at the time of DNAm blood draw, we used lifestyle measurements obtained at or nearest to the DNAm visit, preferring pre-draw values when available; measurements within 12 months prior to and up to 3 months after the DNAm draw were eligible, with the value closest in time selected when multiple observations were available. Due to missing HEI, two analysis cohorts were used; a larger cohort with smoking, physical activity, and BMI (N = 2037) and an HEI-complete subset with all four lifestyle variables (N = 1382).

DNAm bioinformatic analysis
Genome-wide DNAm data were generated using Illumina Infinium MethylationEPIC BeadChip v1 (Illumina, San Diego, CA, USA), with DNA extracted from PBMCs. Preprocessing was conducted using the minfi v1.46.0 package43 in R, following established quality control procedures including a detection p-value ≤ 0.01. CpGs were excluded if they were located on sex chromosomes, contained single nucleotide polymorphisms, or were previously identified as cross-reactive or non-specific. Quantile normalization was applied at the probe level to ensure comparability across samples. After these steps, 752,683 CpGs remained.

To adjust for potential confounding, we next performed a CpG-specific regression for each of the 752,683 CpGs. Specifically, the M-value at each CpG was regressed on age at blood draw for DNAm, sex, leukocyte subtype proportions, the top ten genetic PCs, and the top four DNAm PCs (as described above). The resulting residuals were then used in subsequent analyses. Age at DNAm assessment was included because it is highly correlated with both time since primary cancer diagnosis and time since treatment completion (Pearson r = 0.88 and r = 0.85, respectively; both p < 2.2 × 10−16), thereby serving as an effective proxy for latency since diagnosis and treatment. Including age at DNAm in the residualization step helps account for interindividual differences in latency, ensuring subsequent models capture variation not attributable to heterogeneity in time at blood draw for DNAm since diagnosis.

Gene expression bioinformatic analysis
Gene expression data were generated from PBMC RNA sequencing, with raw counts quantified using htseq-count44 and aligned to the hg38 human reference genome. Gene expression levels were converted to transcripts per million (TPM), and genes with a median TPM ≤ 1.0 across samples were excluded (considered to have very low or negligible expression). Matched DNAm and RNA sequencing data from the same blood draw were available for 199 individuals. To ensure consistency in downstream analyses, log2-transformed quantile-normalized TPM values were residualized using the same covariates applied in DNAm analyses, including age at blood draw, sex, leukocyte subtype proportions, and genetic ancestry PCs, with an additional adjustment for RNA sequencing batches, and top five gene expression PCs.

CpG-gene mapping
For EWAS functional analyses, CpGs were mapped to their nearest genes using annotations from the Illumina EPIC v1 manifest. For expression quantitative trait methylation (eQTM) analysis, CpGs were mapped to genes based on proximity to the TSS within a ±1 megabase (Mb) window, a threshold commonly used in eQTM studies to capture potential cis-regulatory interactions45,46. Of 12,231 genes analyzed, 7099 contained at least one CpG within the 1 Mb TSS window. In total, 14,513 CpG-gene pairs were identified, with each CpG mapped to a median of two genes (mean 2.05, range [1-9]). Conversely, each gene was associated with a median of six CpGs (mean 8.79, range [1-42]).

Statistical analysis
Selection bias assessment. We evaluated potential selection bias arising from exclusions by modeling the probability of exclusion (1 = no WGS or missing RT details; 0 = included) using logistic regression with pre-exposure covariates (sex, age at diagnosis, diagnosis group). In a secondary analysis, we additionally included chemotherapy class indicators to test whether exclusion varied by treatment exposure or follow-up time. No measured covariate was associated with exclusion (all two-sided p > 0.10, Supplementary Data 8), and estimates were essentially unchanged after adding treatment indicators (Supplementary Data 9).

CMRF EWAS. To identify differentially methylated CpGs associated with CMRFs, we performed an EWAS for each condition using logistic regression models. Analyses included all 752,683 CpGs that passed quality control. Each model treated the binary CMRF outcome as the dependent variable, with the residualized M-value for each CpG standardized (i.e., centered and scaled across samples) as the independent variable, ensuring effect sizes are expressed per one standard deviation change in the methylation level for each CpG. Covariate adjustments included age at the most recent follow up, sex, and relevant treatment exposures. Each treatment class was considered individually and was modeled as a binary variable (1 = exposed, 0 = unexposed), with survivors unexposed to that treatment serving as the reference group. Because many survivors received multiple treatments, other treatment exposures were included as covariates in the same model to account for co-exposure and ensure that associations reflect the independent effect of each treatment class. To minimize collinearity and over-adjustment, co-exposures were included only if (i) they were associated with the respective CMRF in the baseline clinical model (see below) and (ii) showed low correlation with the treatment of interest (Phi < 0.4). Childhood cancer treatments are often administered in standardized combinations determined by diagnosis and protocol and including all treatment classes simultaneously would include high collinearity and obscure interpretable associations. Genome-wide significance was set at p < 9 × 10−8, consistent with established threshold in EWAS to account for multiple testing while preserving statistical power. A summarized statistical analysis workflow can be found in Supplementary Fig. 16.

A sensitivity analysis with lifestyle variables was performed for each CMRF. Because lifestyle data were incomplete, we evaluated the impact of lifestyle adjustments in two cohorts: a larger lifestyle cohort with smoking, physical activity, and BMI available (N = 2037) and an HEI cohort with all four variables including the HEI (N = 1382). Building on top of the same base model for the primary EWAS, we fit a prespecified sequence of nested models that added lifestyle covariates sequentially—smoking, then physical activity, then BMI (BMI omitted when the outcome was obesity); in the HEI cohort, HEI was added last. Lifestyle measurements followed the timing rules described above (closest to DNAm; ≤ 12 months prior and ≤ 3 months after). We summarized the impact of lifestyle by (i) retention of genome-wide significant CpG-CMRF associations at each step (same threshold as primary EWAS: p < 9 × 10−8), and (ii) concordance of effect estimates between adjacent sequential models using Spearman correlations.

CMRF-treatment associations. To assess associations between cancer treatments and CMRFs, we performed logistic regression analyses for each treatment-CMRF pair, with the CMRF as the dependent variable and the treatment of interest as the independent variable. A model adjusting for sex and age at the most recent follow-up, and for other treatments whose Phi coefficient was below 0.4 with the exposure treatment, thereby minimizing collinearity while preserving distinct treatment effects5. Because the variables were coded as binary, the Phi coefficients here are numerically equivalent to the Pearson correlation coefficients (r). Significant treatment-CMRF associations were used for subsequent analyses (Supplementary Fig. 10).

Pathway enrichment analysis. Pathways enrichment analysis was performed to identify biological pathways associated with CMRF-associated CpGs and eQTMs using the enrichR 3.2 R package47. To maximize biological interpretability, we incorporated multiple curated pathway databases covering diverse cellular functions and disease mechanisms, including Reactome Pathways 2024, WikiPathways 2024 Human, BioPlanet 209, KEGG 2021 Human, Elsevier Pathway Collection, MSigDB Hallmark 202, BioCarta 2016, HumanCyc 2016, GO Biological Process 2023, GO Cellular Component 2023, GO Molecular Function 2023, MSigDB Oncogenic Signatures, and dbGaP. These databases were selected for their relevance to molecular signaling, metabolic processes, transcriptional regulation, and disease-related pathways, ensuring comprehensive coverage of potential mechanisms linking DNAm to CMRFs. Pathways were considered significantly enriched if they met a Benjamini–Hochberg (BH)-adjusted p-value threshold of < 0.05.

CpG enrichment in genomic features. To evaluate whether CMRF-associated CpGs were over- or underrepresented in specific genomic regions, we performed enrichment analysss using gene annotations from the Illumina EPIC v1 manifest. This analysis aimed to determine whether differential methylation preferentially occurred in regulatory regions, potentially influencing gene expression and cardiometabolic risk. We tested six CpG island annotations (island, north shore, north shelf, south shelf, south shore, open sea) and eight genomic functional regions (TSS1500, TSS200, 5′ UTR, first exon, gene body, 3′ UTR, exon-intron boundaries, intergenic). Enrichment was assessed using Fisher’s exact test, comparing the distribution of CMRF-associated CpGs to all CpGs included on the EPIC array as the background set. Odds ratios and 95% confidence intervals were calculated to quantify effect sizes. Multiple-testing correction was applied using the BH procedure across all 14 tested regions to account for multiple comparisons.

Targeted analysis of treatment-associated DNAm. To assess the effect of cancer treatments on DNAm, we performed an association analysis restricted to the CMRF-associated CpGs. Associations between residualized DNAm and treatment exposures were tested using linear regression models, with each treatment modeled separately as a binary independent variable (treated vs. untreated). Models adjusted for sex and non-correlated treatments (Phi coefficient < 0.4) to account for potential confounding while minimizing collinearity. BH correction was applied to adjust for multiple comparisons across all CpG-treatment tests, with a significance threshold set at FDR < 0.05.

We focused on treatment rather than primary cancer diagnosis because in SJLIFE, treatment regimens are strongly determined by primary diagnosis, yielding high collinearity between diagnosis and treatment exposures. Logistic regression models predicting each treatment from diagnosis yielded highly significant associations (BH-adjusted p < 1 × 10−30; Supplementary Data 5), indicating that diagnosis and treatment are not independently distributed. To avoid multi-collinearity, we focused our analysis on treatment-specific DNAm signatures, consistent with previous literature48.

Mediation analysis. To test whether DNAm mediated the effect of treatment on CMRF risk, mediation analysis was performed on treatment-CMRF-CpG trios that met three criteria: the CpG was significant in both the CMRF EWAS and treatment EWAS, and the treatment was significantly associated with the CMRF in a clinical model (described earlier).

Mediation analyses were performed to evaluate whether DNAm statistically mediated the associations between treatment exposures and subsequent cardiometabolic outcomes. All treatments occurred in childhood or adolescence, DNAm was measured years later in adulthood (median 21.5 years after treatment, IQR 13–30 years), and outcomes were restricted to incident cases that developed after the DNAm assessment (i.e., they were not present prior to or at DNAm assessment and were first detected at a subsequent clinical visit). This temporal structure establishes a clear sequence from treatment to DNAm to outcome.

Confounder choice was guided by a DAG from the treatment (X) to DNAm (M) and further to CMRF (Y) (Supplementary Fig. 17). The DAG defined two sets: exposure-outcome confounders denoted as C(XY) and mediator-outcome confounders denoted as C(MY). C(XY) included sex, age at most recent follow-up, relevant treatment co-exposures (treatments associated with the CMRF but not highly correlated with the exposure of interest; Phi < 0.4; same as in the CMRF EWAS), and genetic ancestry. Primary cancer diagnosis was not included because it is highly collinear with treatment modality. C(MY) included age at DNAm blood draw, sex, relevant treatment co-exposures, leukocyte subtype proportions, and genetic ancestry. Lifestyle factors (smoking, physical activity, BMI) measured before DNAm were treated as potential C(MY); to avoid blocking putative indirect pathways, these were evaluated in additional prespecified sensitivity analyses rather than in the primary mediation models. This DAG-based approach aims to control confounding while avoiding over-adjustment for other variables plausibly on the causal pathway or affected by the exposure, ensuring valid estimation of direct and indirect effects.

Follow-up time was included as a covariate to account for variation in observation windows. We defined follow-up time as the interval (in years) from the date of blood draw for DNAm to each participant’s last clinical assessment. To assess whether differential follow-up time could bias incident-event analyses, we compared follow-up time by incident status (1 vs. 0) within each CMRF using Kruskal–Wallis tests with BH adjustment. Participants who went on to develop a new or worsening CMRF after blood draws were, on average, followed longer than those who have not yet for several conditions (e.g., abnormal glucose metabolism, obesity: BH-adjusted p = 0.041). To mitigate potential bias, models included follow-up duration, defined as the interval between blood draw for DNAm and the most recent clinical visit, ensuring that estimated indirect effects via DNAm were not driven by unequal at-risk periods.

We quantified loss to follow-up (no clinical visit after DNAm) by baseline CMRF status (present vs. absent) using two-sample tests of proportions within each condition and overall; participants with a baseline CMRF were consistently less likely to be lost to follow-up (overall 21.5% vs. 35.4%, Chi-square p = 2.2 × 10−36), with significant differences for abnormal glucose, hypertriglyceridemia, hypercholesterolemia, and hypertension, but not obesity.

Mediation was performed using the mediate function in the mediation v4.5.0 package in R21, with 1000 bootstrap simulations to estimate confidence intervals. The mediation framework combines two interrelated regression models (mediator and outcome) to estimate indirect (ACME) and direct (average direct effect, ADE) effects. The mediator model, capturing how DNAm depends on treatment exposure, was a linear regression with DNAm as the dependent variable and treatment of interest, sex, and relevant co-exposures (as used in the CMRF EWAS) as predictors. The outcome model, capturing how the CMRF depends on both treatment and DNAm, was a logistic regression with the binary incident CMRF outcome as the dependent variable. Predictors included DNAm residuals, sex, age at the most recent clinical follow-up, follow-up duration (included as a covariate to adjust for unequal observation time), and relevant co-exposures. The two models were jointly analyzed using the mediation framework21, which integrates their parameter estimates to quantify how much of the treatment-CMRF association operates indirectly through DNAm versus directly through pathways independent of DNAm.

Multiple testing was controlling using the BH procedure. A CpG was considered a strong mediator if it had a BH-adjusted p-value for the ACME < 0.05 and a proportion mediated > 10%. The 10% threshold was selected to balance sensitivity and biological relevance.

To assess the robustness of mediation findings, we conducted two sensitivity analyses. First, we repeated mediation after restricting outcomes to incident-only events, defined as CMRFs that developed strictly after the DNAm blood draw, thereby excluding worsened cases. Incident-only models included the same covariates as the primary analysis. Second, we evaluated robustness to lifestyle adjustments. Guided by EWAS sensitivity results (which showed modest, similar attenuation for smoking and physical activity and the largest attenuation for BMI), we fit nested mediation models: (i) base model (primary analysis), (ii) base+smoking+physical activity, and (iii) base+smoking+ physical activity+BMI (applied to non-obesity outcomes). Diet quality (HEI) was not included due to minimal impact in EWAS as well as substantial missingness (n = 655 [32%] missing HEI data). For trios significant in the primary analysis, robustness was evaluated by changes in proportion mediated, relative change in ACME (primary-incident)/primary, Spearman correlation of estimates, and concordance and retention of significant trios.

Composite mediation analysis. For treatment-CMRF pairs with multiple strong mediating CpGs, a composite mediation analysis was performed to assess the cumulative effects of DNAm. CpGs were first ranked by the estimated ACME in descending order. To reduce redundancy, CpGs were pruned using pairwise Pearson correlation (r2 < 0.05), retaining the most statistically significant mediator within each correlated cluster to ensure independent signals. Following pruning, composite DNAm scores were generated by summing z-score standardized M-values of the remaining independent CpGs within each treatment-CMRF pair. These composite scores were then used as single mediators in models with the same structure as the individual CpG analyses, adjusting for age, sex, and treatment co-exposures associated with the CMRF but not the exposure treatment. Multiple comparisons were controlled using BH correction.

eQTM analysis. To investigate potential gene regulatory mechanisms, an eQTM analysis was conducted using 199 participants with matched DNAm and RNA sequencing data from the same blood draw. All CMRF-associated CpGs were tested for association with expression levels of nearby genes, defined using a ± 1 Mb window around the gene’s TSS. This window was chosen based on prior studies45,46 demonstrating that the majority of cis-acting methylation-expression relationships occur within this range. Linear models were used to test associations with residualized gene expression as the dependent variable and residualized DNAm as the independent variable. Multiple testing correction was applied across all eQTM tests using the BH procedure, with significance defined as FDR < 0.05.

Statistical software. All analyses were conducted in R 4.4.1, with visualizations generated using the ggplot2 v3.5.1 package49. Circle Manhattan plots were creating using CMplot from the rMVP v4.5.1 package50, while UpSet plots and heatmaps were generated using ComplexHeatmap v1.3.351. Logistic regression models were fitted using the glm function in R with a logit link function, and linear regression models were fitted using the lm function in R.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The study by Eulalio et al. investigates the link between genotoxic cancer treatments and the development of cardiovascular disease (CVD) in childhood cancer survivors from the SJLIFE cohort (n=2938), focusing on the role of DNA methylation (DNAm) as a molecular mediator. The authors build on previous research from the same cohort (Song et al. 2021; Dong et al. 2023). Their approach integrates epigenome-wide association studies (EWAS) to identify CpG sites linked to cardiometabolic risk factors (CMRFs) and to determine those CpGs potentially influenced by cancer treatments. They further employ causal mediation analysis in the context of available followup data and cis-eQTM mapping to pinpoint CpGs that are not only significantly associated with CMRFs but also provide insight into underlying biological mechanisms. After validation, these CpGs may have translational potential as biomarkers for risk stratification and as therapeutic targets to improve outcomes for cancer survivors.\nI believe the authors made, generaly, an excellent and methodologically sound work. And the topic is definitely of interest and relevance. My main criticisms are related to fundamental study design aspects that get lost or hiden by the large amount of results being reported, which should be made more visible/available to the readership in the text. More detailed comments follow below.\n1) What are the noteworthy results?\nThe main message of this study is that genotoxic cancer treatment induces DNA methylation (DNAm) changes, which are in turn associated with increased cardiometabolic risk factors (CMRFs) among childhood cancer survivors. The study identifies around 2000 CpGs associated with CMRF prevalence, with only 9% shared among risk factors—primarily in genes involved in metabolic and inflammatory pathways. They also identify an obesity-specific epigenetic signature in survivors and two genes (ANTRX2 and LRIG1).\nDespite these interesting findings, the description of cardiometabolic outcome definitions and the follow-up design is insufficient to assess the potential for bias, specially given the emphasis that authors pose on the prospective causal mediation analysis. Important questions are left unanswered:\n• Were standard CMRF definitions used (e.g., ADA, WHO, ATP III)?\n• Were baseline and follow-up criteria consistent?\nOften, baseline definitions involve laboratory measures, while follow-up may rely on registry data, which affects ascertainment quality.\n• What was the attrition rate across follow-up?\n• Were losses to follow-up differential in terms of baseline cardiometabolic risk?\n• Was follow-up time variable across participants, and how was this handled?\nDifferential attrition and variance in follow-up might eventually introduce bias and should be explicitly addressed using epidemiologic best practices.\n2) Will the work be of significance to the field and related fields?\nThis study contributes meaningfully to survivorship research by linking cancer therapy → epigenetic alterations → CMRFs using prospective mediation analysis—a largely unexplored area. While previous studies have examined treatment impacts on DNAm and separate works have linked DNAm to cardiometabolic health, this is among the first to integrate all three into a coherent causal framework using prospective data.\nThe use of causal mediation analysis with DNAm is innovative, but the authors must clearly discuss key assumptions such as:\n• No unmeasured confounding between exposure–mediator and mediator–outcome,\n• No exposure–mediator interaction,\n• Multiple mediators being independent, an assumption questioned when CpGs are spatially correlated.\nThe authors should discuss these assumptions thoroughly and, if possible, conduct sensitivity analyses. Recent methodological reviews on high-throughput mediation and comparative frameworks of omics mediation stress the importance of addressing these issues and conducting sensitivity analyses to validate causal inference.\n3) Does the work support the conclusions and claims, or is additional evidence needed? In my opinion, the conclusions are appropriately drawn and qualified, with the authors clearly noting the need for external validation, especially in multiethnic cohorts. The use of the EWAS catalog to validate findings externally is a strength, though not equivalent to cohort-level validation.\nNevertheless, smoking and inflammation, both influential on DNAm and cardiovascular outcomes, need more careful evaluation.\n• Recent evidence shows smoking-linked DNAm associations in atherosclerotic tissue and blood, closely related to CVD risk.\n• Similarly, epigenetic trends tied to environmental/behavioral risk factors illustrate smoking-induced CVD epigenetic mechanisms.\nThe authors should add sensitivity analyses, clarify whether these factors were/should be included in adjustment models, and frame their rationale according to established cardiovascular epidemiology (e.g., adjustment sets including age, sex, physical activity, diet, total cholesterol, HDL cholesterol, lipid- and blood pressure-lowering medications, diabetes, systolic blood pressure, BMI, and renal disease).\nThe rationale for including/excluding each variable per CMRF outcome must be clearly provided.\n4) Are there any flaws in data analysis, interpretation and conclusions?\nThere are no fatal flaws; the authors present the limitations appropriately, especially the time lag between cancer diagnosis and DNAm sample collection, and the relatively short period between DNAm measurement and clinical follow-up. This timing affects causal interpretation, but the authors acknowledge it.\nIn high-dimensional omics research, cross-validation or regularization is essential to avoid overfitting. This manuscript does not report using methods such as k-fold cross-validation or bootstrap procedures. Terms suggesting clinical “precision medicine” should be tempered unless cross-validation is available.\n5) Is the methodology sound?\nYes, overall the methodology is strong, but software and package versioning (e.g., minfi, mediation) should be added to ensure full reproducibility. The methods section is quite thorough, with code and GEO data available. With version numbers included, the work is reproducible.\n*Final Recommendation*\nThis study is valuable and addresses an important gap. However, the following moderate revisions are required for rigor and clarity:\n1. Provide detailed follow-up information and harmonized cardiometabolic outcome definitions grounded in accepted guidelines.\n2. Evaluate potential confounding by smoking and inflammation, including sensitivity analyses; clarify adjustment sets per outcome.\n3. Fully articulate causal mediation assumptions, and introduce sensitivity checks.\n4. Implement or acknowledge cross-validation/regularization for high-dimensional DNAm signatures; revise predictive language accordingly.\nOnce addressed, this manuscript will offer a robust, insightful contribution to cancer survivor health research.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "The study by Eulalio et al. investigates DNA methylation CpG sites associated with five cardiometabolic risk factors (CMRFs) and various cancer treatments, including chemotherapy and radiation, in a cohort of 2,938 childhood cancer survivors from the St. Jude Lifetime Cohort. The researchers employed a comprehensive approach to exploring biological relationships, integrating epigenome-wide association studies (EWAS), causal mediation analysis, and cis-expression quantitative methylation (cis-eQTM) mapping. The study provides valuable insights into refining risk stratification and shaping targeted interventions in survivorship care, which could significantly enhance long-term outcomes. Among the many strengths of this research is: 1. The thorough description of the elegant methodologies and applicable statistical analyses employed throughout the study and 2. The limitations of the study are clearly articulated, offering a clear foundation for future research directions. Overall, the soundness of the research reported is outstanding and the results are noteworthy.\nMinor Comments:\nIn line 137, the percentage of survivors with hypertension should be revised to 29.2% according to the table, contrasting with the 29.3% currently indicated in the results section.\nAccording to Figure 2c, the total numbers of unique and overlapping CpGs results in 1,894 CpGs instead of the 1,893 stated in the text in line 154. Please double check.\nClarification is required for the number of 1,912 CpG-CMRF associations mentioned in line 154. The total of associated CpGs with each CMRF does not align with this number in Fig 2c.\nThe section titled “Relevance of CMRF-associated CpGs in survivors relative to the general population” includes paragraphs starting with “Among the 1,720 obesity-associated CpGs, 155 (9%) had been linked to body mass index,…” which belong more appropriately in the preceding results section.\nAuthors are encouraged to name the three CpGs mapped to CPT1A in line 193-194.\nIn the section that starts on line 198, the reader may benefit from text that explicitly states that cancer type was not examined for associations with widespread DNAm variation, like cancer treatment exposure was. This question may come up for some readers so including the related reasoning for this could be helpful.\nThere is a discrepancy between the stated 2,591 significant treatment-CpG associations in line 204 and the total of numbers of N in the bar plot in Fig 3 (2,574).\nSeveral sentences within lines 319 to 338 of the discussion repeat content from the results section. Rewriting these sections to avoid redundancy is recommended.\nThe limitation listed in lines 396-398 is key as the authors acknowledge that even with adjusting for key factors, residual confounding may still influence DNAm and CMRFs. The importance of lifestyle factors (when it comes to future interventions) cannot be overstated so the reader may benefit from some additional text (somewhere in the Discussion) that provides discussion on how this specific factor could potentially influence DNAm and CMRFs and what this could mean for future research.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "This thoroughly analyzed and described study analyzed DNA methylation in PBMCs and cardiometabolic risk factors (CMRFs) among childhood cancer survivors in the St. Jude Lifetime Cohort. The authors identified several methylation sites associated with CMRFs and conducted mediation analyses to estimate direct and indirect causal effects of methylation sites between treatments and CMRFs. The following points need further clarification/justification.\nDiscussion\n• While the authors address generalizability related to race, generalizability related to the cohort source should also be discussed.\n• Discussion of potential misclassification of CMRFs should be included.\nMethods\n• P. 21, Participants: what were the reasons for missing chemotherapy and radiotherapy data? Describe whether this missing data may have introduced a selection bias in the Discussion.\n• P. 21, CMRF data: provide additional detail about the methods for obtaining the CMRF data. Did multiple reviewers classify CMRFs? If so, were reliability and validity evaluated?\n• P. 22, lines 451-452: consider a sensitivity analyses restricting to only those survivors with CMRFs that newly developed after the blood test (excluding those that worsened) to further establish temporality.\n• P. 23, Covariates: because of the focus on causal models and causal mediation analyses, a directed acyclic graph should be used to identify potential confounders for the primary analyses and to fully evaluate the assumptions of mediation analyses.\n• P. 25, Statistical Analysis: how were treatments modeled? Were they categorized as in Table 1 or modeled as binary variables (e.g., chemotherapy yes/no)?\n• P. 27, Mediation Analysis: Details about the methods for conducting the causal mediation analyses, including the required assumptions about confounders, need to be provided. See this reference for discussion of the required assumptions: Naimi AI, Schnitzer ME, Moodie EEM, Bodnar LM. Mediation Analysis for Health Disparities Research. Am J Epidemiol. August 15, 2016 2016;184(4):315-324. doi:http://dx.doi.org/10.1093/aje/kwv329"}, "author_response": "We thank the reviewer for the overall favorable assessment of our study including the methodological rigor of data analyses and translational relevance of the findings. We have now revised the manuscript to address the reviewer’s primary critique regarding the clear delineation of study design from novel findings by enhancing clarity and visibility of the overall study design. Specifically, we have (i) added a concise study-design description in the final paragraph of the Introduction (citing Figure 1) and (ii) consolidated dispersed design descriptions into a dedicated “Study design overview” paragraph in the Methods, with step-by-step alignment to visual illustrations in Figure 1. We appreciate the Reviewer’s feedback and believe these changes make the analytical workflow clearly distinguishable from novel findings and hence substantially improve manuscript readability.\n\n[...]"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made a substantial effort to address my previous comments. Consequently, the manuscript has improved in clarity and rigor.\nHowever, I still have one remaining and critical concern, as several key aspects of the study design remain unclear. The additional explanations provided in this revised version to clarify the design are still ambiguous.\nThe authors state that follow-up was defined as the time from DNA methylation (mediator) assessment to incident outcome. This definition remains conceptually problematic in the mediation analysis. The follow-up should not originate from the mediator, but rather from a time point that precedes the mediator (e.g., treatment) in the (I presume, as not much information on the outcome and mediator models is given) time-to-incident CMRF model. Consequently, the outcome model in the mediation analysis needs to be corrected, with some considerations.\nSpecifically, the authors define treatment as: “All chemotherapy classes and RT exposures were defined as binary variables indicating whether a given treatment was administered within the first five years following cancer diagnosis. Each of the eleven treatment classes was considered separately.”\nThis clarification is helpful, but it raises additional questions:\nWhat is the reference category for each categorical treatment variable? Are comparisons made against individuals who did not receive that specific treatment (i.e. against those who received other treatments than the one coded as 1)? The reference category for interpretation of results must be unambiguous in the text. Importantly, how did the authors handle participants who received multiple treatments?\nAlso, were recurrences or additional treatments after the first five years considered? If participants experienced recurrence, the time of the last completed treatment might be a more appropriate time origin for analyses. A sensitivity analysis considering this possibility would greatly strengthen the study.\nImportantly, there may still be heterogeneity in time since diagnosis across individuals at the time of inclusion, which could introduce bias. Clarifying how time since diagnosis was handled in the analyses is essential.\nWithout a clearly defined temporal sequence—from treatment (exposure) to DNA methylation (mediator) to outcome it is impossible to determine whether the assumptions underlying the mediation analysis are satisfied. Given the uncertainties regarding the study design and the modeling of treatment in the mediation analyses, the authors should refrain from using causal language throughout the manuscript. All references to “mediation” or “causal pathways” should be reframed to describe statistical associations only.\nClarifying these points will ensure that the interpretation of the analyses remains valid and consistent with the underlying study design.\nRelatively minor comment:\nRegarding Figure 1 (“Study design for identifying and characterizing mediator CpGs linking treatment exposures to CMRFs in survivors”): this figure does not fully represent the epidemiologic study design, but rather the statistical analysis strategy. If the authors agree with this interpretation, I recommend revising the title of Figure 1 to replace “study design” with “statistical analysis strategy”.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_3": "The authors have addressed all the comments adequately.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_5": "The study analyzed DNA methylation in PBMCs and cardiometabolic risk factors (CMRFs) among childhood cancer survivors in the St. Jude Lifetime Cohort. The authors identified several methylation sites associated with CMRFs and conducted mediation analyses to estimate direct and indirect causal effects of methylation sites between treatments and CMRFs. This work is significant to the field and will further support future intervention development for high-risk childhood cancer survivors. The work, including the detailed methods and sensitivity analyses in the revised manuscript, support the rigor and reproducibility of the analyses.\nThe revised manuscript thoughtfully considered feedback from the initial reviewers and thoroughly addressed the previous critique. I have no further suggestions for improvement."}, "author_response": "We really appreciate Reviewer 1’s recognition of the improvements in our R1 submission. We also acknowledge that the R1 submission of our manuscript may not have fully conveyed the timing and structure of our data analyses, and we have now clarified these aspects extensively throughout the Methods. In this R2 submission, we carefully revisited each point and provided additional details to make the temporal sequence, modeling framework, and analytic assumptions fully transparent. We believe these revisions have further improved the clarity and overall quality of the manuscript, and we hope that the enhanced descriptions about our study design and corresponding analytical approach address the remaining concerns to the Reviewer’s full satisfaction.\n\n[...]"}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have made a tremendous effort to address my comments. As a result the manuscript has significantly increased in rigor and transparency.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We appreciate the Reviewer’s highly positive feedback on the revised submission."}]

Decision letter