Epigenome-wide analysis reveals potential biomarkers for radiation-induced toxicity risk in prostate cancer

Abstract

Background

Prostate cancer is the second most common cancer globally, with radiation therapy (RT) being a key treatment for clinically localized and locally advanced cases. Given high survival rates, addressing long-term side effects of RT is crucial for preserving quality-of-life. Radiogenomics, the study of genetic variations affecting response to radiation, has primarily focussed on genomic biomarkers, while DNA methylation studies offer insights into RT responses. Although most research has centred on tumours, no epigenome-wide association studies have explored peripheral blood biomarkers of RT-induced toxicities in prostate cancer patients. Identifying such biomarkers could reveal molecular mechanisms underlying RT response and enable personalized treatment.

Methods

We analysed 105 prostate cancer patients (52 cases and 53 controls). Cases developed grade ≥ 2 genitourinary and/or gastrointestinal late toxicity after 12 months of starting RT, whereas controls did not. An epigenome-wide association study of post-RT toxicities was performed using the Illumina MethylationEPIC BeadChip, adjusting for age and cell type composition. We constructed two methylation risk scores—one using differentially methylated positions (MRSsites) and another using differentially methylated regions (MRSregions)—as well as a Support Vector Machine-based methylation signature (SVMsites). We evaluated RT effects on biological age and stochastic epigenetic mutations within established radiation response pathways. Gene Ontology and pathway enrichment analyses were also performed.

Results

Pre-RT methylation analysis identified 56 differentially methylated positions (adjusted p-value ≤ 0.05), and 6 differentially methylated regions (p-value ≤ 0.05) associated with the genes NTM, ACAP1, IL1RL2, VOOP1, AKR1E2, and an intergenic region on chromosome 13 related to Short/Long Interspersed Nuclear Elements. Both Methylation Risk Scores (MRSsites AUC = 0.87; MRSregions AUC = 0.89) and the 8-CpG Support Vector Machine signature (SVMsites AUC = 0.98) exhibited strong discriminatory accuracy in classifying patients in the discovery cohort. Gene ontology analysis revealed significant enrichment (adjusted p-value ≤ 0.05) of genes involved in DNA repair, inflammatory response, tissue repair, and oxidative stress response pathways.

Conclusions

Epigenetic biomarkers show potential for predicting severe long-term adverse effects of RT in prostate cancer patients. The identified methylation patterns provide valuable insights into toxicity mechanisms and may aid personalized treatment strategies. However, validation in independent cohorts is essential to confirm their predictive value and clinical applicability.

Background

Prostate cancer is the second most common cancer worldwide, with an incidence of 14.1% (n = 1,414,259) in 2022, only 0.2 points behind lung cancer (14.3%, n = 1,435,943). It is estimated to affect over two million men by 2035 [1]. Radiation therapy (RT) is one of the primary treatment protocols for the management of clinically localized and locally advanced prostate cancer with curative intent [2]. Advances in RT are occurring rapidly, making it increasingly precise and allowing for a reduction in the radiation dose delivered to surrounding health structures. Nevertheless, the five-year incidence of grade ≥ 2 gastrointestinal (GI) and genitourinary (GU) toxicities are still around 20% and 40%, respectively [3].

In 2023, Hamdy et al. revealed that only 3% of RT patients had died of prostate cancer after a 15-year follow-up [4]. Considering these survival rates, it is crucial to weigh the long-term effects of treatments to ensure a quality-of-life after undergoing RT. Therefore, optimal patient management and precise personalisation of RT is an urgent need. Radiogenomics, which involves the study of genetic variation associated with the response to RT [5], is emerging as a key area of research. Since its inception, research in this field has primarily focussed on the identification of genomic biomarkers. Initially, this was conducted through candidate gene studies [6,7,8] and later through Genome-Wide Association Studies (GWAS) [9,10,11]. However, despite these efforts, the results obtained to date are still limited and lack clinical applicability for accurately personalizing patients’ treatment based on their post-RT toxicities and future quality-of-life.

With the emergence and consolidation of recent omics technologies, radiogenomics has found new and promising approaches to elucidate the response to RT. DNA methylation, which involves the chemical modification of cytosines and cytosine-guanidine dinucleotides, is an inheritable epigenetic mark known for its profound repressive effects. In 2010, Kim et al. compared the epigenetic profiles of two lung cancer cell lines, one radioresistant and one radiosensitive, identifying 1,091 genes with differential methylation rates [12].

Subsequent efforts have primarily focussed on characterizing tumours in their response to RT through epigenetic patterns, with a particular emphasis on methylation. Most findings have emerged from candidate gene studies aimed at identifying differentially methylated promoters [13]. For instance, Smrdel et al. described the hypermethylation of the methylguanine-methyltransferase promoter as a favourable prognostic factor for glioblastoma, as it reduces high methylation rates among tumour cells, thereby enhancing radiosensitivity [14]. Consequently, hypomethylating drugs, such as 5-aza-dC, are being developed for co-administration with RT to increase treatment efficacy [15].

While most published epigenetic studies have followed a candidate gene approach, epigenome-wide association studies (EWAS) represent a different but applicable approach for identifying differential methylation associated with response to RT. Analogous to GWAS, EWAS aims to establish associations between epigenetic markers and phenotypes, potentially elucidating the causes of diseases and promoting the development of new therapies and diagnostic methods [16]. Despite its potential and feasibility, only a few EWAS studies in the field of radiogenomics are found in the literature. Among these, only two studies have been performed on peripheral blood after RT, identifying methylation biomarkers associated with cardiometabolic diseases [17] or affecting survival and time to relapse in head and neck squamous cell carcinoma patients [18].

In this study, we present the first EWAS of adverse effects from RT in prostate cancer patients. Our findings offer valuable insights into the molecular mechanisms underlying genetic variation in radiation sensitivity, which is crucial for ensuring effective and safe treatment. Furthermore, the identification of candidate epigenetic biomarkers has the potential to stratify patients according to their susceptibility to these adverse effects.

Methods

Aim

The objective of this study was to identify epigenetic markers in peripheral blood associated with toxicity in prostate cancer patients undergoing RT. These biomarkers could potentially aid in predicting late-onset GI and GU toxicities in these patients.

Study subject and design

Samples for this study were obtained from prospective selected Galician prostate cancer patients who were treated at the Radiation Oncology Department of the Hospital Clínico Universitario de Santiago de Compostela (northwest Spain). Patients treated with radical, salvage or adjuvant RT received a total dose for the planning target volume (PTV) in the ranges of 70–76 Gy, 66–70 Gy and 60–66 Gy, respectively. All individuals were recruited into the RADIOGEN-Prostate study, initiated in 2005, which collects blood samples, dosimetry data, clinical information and toxicity data [7, 8]. Written informed consent was obtained for each subject in accordance with the protocols approved by the ethics review board of the Galician Ethical Committee for Clinical Research.

Toxicity evaluation was conducted before the initiation of treatment, at 6- and 12-months post-completion of RT, and annually thereafter, following the Common Terminology Criteria for Adverse Events v4.0 (CTCAE) [19]. From the 1673 patients in the RADIOGEN-Prostate cohort, a subset of 294 individuals irradiated between 2014 and 2018 was included in the initial selection. Of these, 105 individuals meeting the established criteria were selected for analysis, with 52 serving as cases and 53 as controls. Cases were defined as patients who experienced an increase in the maximum CTCAE toxicity score ≥ 2 compared to the pretreatment symptom assessment, from months 12 to 72 after completing RT. Controls were individuals who reported a baseline-adjusted CTCAE toxicity score = 0 throughout the five-year follow-up period. The following GI and GU toxicities were assessed: proctitis, bowel perforation, bowel obstruction, bowel fistula, bowel stenosis, bowel ulceration, diarrhoea, flatus, rectal bleeding, management of sphincter control, haematuria, urinary tract obstruction, urinary incontinence, urinary frequency, urinary urgency, urinary retention, urinary fistula, and urinary obstruction.

DNA isolation and methylation profiling

Blood samples were collected in 10 ml EDTA tubes until processing. DNA isolation was performed using the Chemagic™ Magnetic Separation Module (MSM) I robot (Chemagen Biopolymer-Technologie AG) with magnetic bead technology. The Chemagic DNA Blood 200 Kit (PerkinElmer Inc.) was used following the manufacturer’s protocol.

Bisulphite conversion was performed using the EZ DNA Methylation-Gold Kit (Zymo Research). Processed samples were hybridized to the Illumina Infinium MethylationEPIC BeadChip and scanned using the Illumina iScan platform. The intensity of the images was extracted with the GenomeStudio Methylation Software Module (v1.9.0). The samples' location in the Illumina Infinium Methylation EPIC array was randomized by age and phenotype to minimize batch effects.

Statistical analysis

Clinical characteristics between cases and controls were compared using a generalized linear model. Adjustment for multiple comparisons was performed using the False Discovery Rate (FDR) method.

Methylation data underwent quality control, normalization, and subsequent statistical analyses using the statistical R software and Bioconductor packages. The Chip Analysis Methylation Pipeline (ChAMP), as described by Campagna et al. [20], was employed. Raw intensity files were imported into the R environment using the Bioconductor package ChAMP [21] and subsequently transformed into β and M-values matrices. M-values were employed in our statistical analysis and were calculated as the log₂ ratio of intensities between methylated and unmethylated probes per CpG, yielding values ranging from −1 to 1.

Quality control and normalization

Patients with CpG island detection rates of 96% or higher were included in the study. Methylation data underwent filtering to remove probes that could introduce noise and uncertainty. Low-quality probes (< 3 beads in > 5% of samples) or with detection p-values > 0.01 were filtered out. Additionally, probes that were non-CpG, non-autosomal, multi-hitting, or containing Single Nucleotide Polymorphisms (SNPs) were also removed.

The Illumina Infinium MethylationEPIC BeadChip utilizes two types of hybridisation chemistries: Type I and Type II. We applied the beta-mixture quantile normalization method in ChAMP [22] to normalize the distribution of type I and II probe densities, thereby avoiding biased detection of Differentially Methylated Positions (DMPs). Batch effects were identified using the Singular Value Decomposition method [23]. The methylation matrix was corrected for sequencing round, sample frozen state, and sequencing array using the ComBat method to eliminate the batch effect [24].

Bioinformatic analysis

Differential methylation analyses

B-cells, NK-cells, CD4T- and CD8T-cells, Monocytes, Neutrophils and Eosinophils cell type proportions (CTp) were estimated with the Robust Partial Correlations algorithm of the EpiDISH tool [25] using the centDHSbloodDMC.m reference matrix. To identify DMPs between cases and controls, we employed a logistic regression model in the limma package [26] with M-values as the dependent variable and including age and CTp as covariates. Standard errors were moderated using an empirical Bayes method (eBayes). DMPs were filtered based on a significance threshold of adjusted p-value < 0.05 (Benjamini–Hochberg method [27]). Results were further annotated using the DMRcate package [28]. A Fisher’s exact test following a smoothing Gaussian kernel was applied through the DMRcate package to identify Differentially Methylated Regions (DMRs), defined as regions > 1000 bp comprising at least two CpGs and p-value ≤ 0.05. Principal Component Analysis (PCA) was subsequently performed over the significant DMPs to assess their potential stratification ability.

To perform gene set enrichment analysis (GSEA) on DMPs, we utilized the empirical Bayes GSEA method implemented in ChAMP [29]. Additionally, we employed the R package methylGSA [30] to extract enriched biological functions and pathways from the Gene Ontology (GO) and Reactome databases [31,32,33]. For all three analyses, p-values were adjusted using the Benjamini–Hochberg procedure. We considered the enrichment significant if the FDR was < 0.05.

A detailed analysis was performed to assess whether the enriched biological functions or pathways were predominantly hypo- or hypermethylated in cases. For each gene involved in these functions or pathways, the average methylation change was calculated as the difference between the mean methylation levels in cases and controls.

Assessment of biomarkers’ classification ability

The classification ability of DMPs and DMRs was evaluated through the generation of three models: (1) methylation risk score based on DMPs (MRSsites), (2) methylation risk score based on DMRs (MRSregions), and (3) support vector machine model incorporating the most discriminative DMPs (SVMsites). MRSsites was calculated for each patient as the weighted sum of all DMPs, with weights corresponding to the effect sizes derived from the EWAS. The MRSregions was generated by fitting a logistic regression to the weighted sum of all CpGs within each DMR, as detailed in Eq. 1.

$$MRSgenes = \mathop \sum \limits_{i} \left( {w_{i} \cdot \mathop \sum \limits_{j} \left( {\beta_{ij} \cdot CpG_{ij} } \right)} \right)$$

(1)

where w_i represents the weight assigned to DMR i, β_ij is the effect size derived from the EWAS for CpG j in DMR i, and CpG_ij is the methylation level of each CpG site j in DMR i.

To identify the most relevant CpG sites for the SVMsites model, the ReliefF feature selection algorithm was applied using a 100-repetition bootstrap approach on the identified DMPs. CpGs with a mean importance weight greater than 0.11 were selected for further analysis. To minimize redundancy among the selected features, pairwise correlations were calculated using Pearson’s correlation coefficient, and CpGs with r² > 0.6 were excluded. The resulting set of informative, non-redundant CpGs was then used to build the SVMsites model.

The classification performance of the MRSsites, MRSregions, and SVMsites model was further evaluated using a fourfold cross-validation approach with 2,000 iterations.

Stochastic epigenetic mutations (SEMs) analysis

SEMs refer to CpGs displaying methylation levels in an individual that fall outside the ranges observed across the entire population for that site [34, 35]. We conducted SEMs analysis to assess epigenetic drift at the individual level.

For each probe in the study, we estimated the methylation ranges for the reference population using the formula:

$$\left[ {Q1\left( {k \cdot IQC} \right), Q3 + \left( {k \cdot 1QR} \right)} \right]$$

where Q1 is the first quartile, Q3 is the third quartile, IQR is the interquartile range (Q3–Q1), and k is a multiplicative factor. We quantified all CpGs that fell outside this range for each individual, using two values for k: a mild one (k = 2) and a restrictive one (k = 3).

We applied a logarithmic transformation on the number of SEMs, expressed as log(SEM) for all analyses. To assess the relationship between SEMs per patient and post-RT toxicities, logistic regression models were fitted with the phenotype, age and CTp as covariates. Initially, all SEMs throughout the genome were included in the analysis. Subsequently, the analysis was repeated, considering the number of SEMs separately within specific genomic regions: TSS200 (Region from the Transcription Start Site to 200 bp upstream), TSS1500 (Region from 1500 to 200 bp upstream of the Transcription Start Site), TSS200 + 1500 (Region that includes both the TSS200 and TSS1500), 1st exon, gene body and intergenic regions.

Furthermore, we investigate the influence of SEMs on three pathways associated with RT response: DNA repair, cell senescence, and the immune system. Genes involved in each pathway were obtained from the Reactome database [33]. Subsequently, we conducted a separate analysis focussing specifically on SEMs located within the genomic regions of each pathway.

Epigenetic clock measures

We estimated the age acceleration (AA) as the difference between the predicted methylation age and the actual individual age. To determine methylation age and AA, we used the most recent iteration of Horvath's epigenetic clock [36]. We then compared the mean values of methylation age and AA between cases and controls using a t-test. Additionally, given that SEMs can influence ageing processes, we explored the association between log(SEMs) and both methylation age and AA.

Results

Sample characteristics

The epigenetic analysis of the RADIOGEN cohort included 105 patients (52 cases and 53 controls) diagnosed with non-metastatic prostate adenocarcinoma. These patients underwent radical treatment with a mean total dose of 68.6 Gy in three-dimensional conformal RT (3D-CRT), at an average age of 70.9 years (71.0 years for cases and 70.7 years for controls). No statistically significant differences were observed between cases and controls in terms of clinical, comorbidities, treatments, dosimetry data, or tumour characteristics (Table 1). Although there was an initial observation of significance in the percentage of the rectum receiving a dose of 65 Gy (V65Gy) with a p-value < 0.05, this significance was not maintained after adjustment for multiple comparisons (FDR-adjusted p-value = 0.94).

Table 1 Characteristics of prostate patients in the RADIOGEN cohort

Differential methylation analysis

Out of 850,000 probes, 681,806 successfully passed the quality control filters, with no patients excluded from the analysis. There were no significant differences in the estimation of per-sample CTp, except for CD8T-cells that were higher in controls (t = 2.58; p-value = 0.011) (Fig. 1).

Fig. 1

Comparison of cell type proportion between cases (red) and controls (light blue). The x-axis represents cell types (CT) and the y-axis count proportions. *p-value < 0.05

Adjusting the linear model for all CTp deconvoluted and patients' age, we identified 56 DMPs associated (Fig. 2A) with 43 unique genes. Among these, 15 DMPs were intergenic, while three DMPs (cg00162348: chr16, 30,773,695 bp| cg11679914: chr7, 99,662,323 bp | cg12593138: chr12, 53,439,997 bp) were situated at positions overlapping two genes. Specifically, 34.7% of the DMPs were located upstream of genes (n = 6 at − 1500 bp and n = 11 at − 200 bp), 16.3% (n = 8) in the 5’-UTR, 6.1% (n = 3) in the first exon, 36.7% (n = 18) in the gene body region, and the remaining 6.1% (n = 3) in the 3’-UTR. The PCA performed on the 56 DMPs shows clear patient stratification, with principal components 1 and 2 explaining 19.6% and 11.3% of the variance, respectively (Fig. 2B). Of these 56 DMPs, 31 exhibited higher methylation levels in cases, while the remaining 25 showed higher methylation in controls (Fig. 2C).

Fig. 2

EWAS Results. A Volcano plot, B Principal Component Analysis (PCA) plot, C Heatmap, and D Manhattan plot illustrating the association of DNA methylation CpGs with radiation-induced grade ≥ 2 genitourinary and/or gastrointestinal late toxicities in prostate cancer patients (52 cases and 53 controls). In the volcano plot, dark blue and red dots represent significant hypomethylated and hypermethylated CpG sites, respectively, in cases compared to controls (the grey dashed line represents the cutoffs of adjusted p-value < 0.05 and logFC =|0.2|). The PCA plot demonstrates segregation between cases (red dots) and controls (light blue dots), with the explained variance for each principal component (PC1-2) displayed. The heatmap includes unsupervised clustering of the 56 significant differentially methylated positions (DMPs), with a colour gradient where black denotes low methylation and yellow denotes high methylation. In the Manhattan plot, the X-axis represents chromosome positions, while the Y-axis represents the –log10 adjusted p-values. The gray line indicates the significance threshold at an adjusted p-value < 0.05, with red dots highlighting CpGs belonging to significant differential methylation regions

We identified 6 statistically significant DMRs (Fisher's exact test p-value ≤ 0.05) (Fig. 2D), where adjacent CpG probes exhibited varying methylation levels between cases and controls, as shown in Fig. 3. The DMR associated with the neurotrimin gene (NTM) (p-value = 0.002) encompassed a 185 bp region within the first intron, containing three hypermethylated CpG probes in cases (cg09294095, cg18009484, cg15277677). The second most significant DMR (p-value = 0.006) was the only region to display significant hypomethylation in cases compared to controls. This hypomethylation was consistently observed across all 10 CpG probes (cg08324090, cg06721232, cg25900902, cg07925670, cg20217592, cg15713546, cg02448825, cg02676175, cg07148458, cg07520074) within a 1255 bp region covering exons and introns 19, 20, and 21 of the ArfGAP with coiled-coil domain-containing gene (ACAP1). Another DMR with a length of 50 bp, exhibited hypermethylation in cases across three CpGs (cg21181713, cg13779009, cg08023416; p-value = 0.015) within the first intron of interleukin 1 receptor-like 2 gene (IL1RL2). Another DMR, 149 bp in length containing three CpGs (cg03357999, cg25966908, cg03014829), was identified as hypermethylated in cases (p-value 0.023) within the second intron of VOPP1 WW domain-binding protein (VOOP1) gene. The final DMR within a gene region (p-value = 0.032) spanned 626 bp with 11 CpGs (cg00578885, cg25162403, cg22362895, cg09305334, cg19165274, cg20282550, cg19392551, cg13218425, cg00700969, cg12255897), where all but two CpG probes displayed hypermethylation in cases versus controls within the upstream region and first exon of aldo–keto reductase family 1 member E2 gene (AKR1E2). The sixth DMR (p-value = 0.018), situated within an intergenic region, encompassed 10 CpG probes (cg26361286, cg15973954, cg01863042, cg08083251, cg03042692, cg05215994, cg20395040, cg18959621, cg05656990, cg17136073) exhibiting higher methylation rates in cases (Table 3, Fig. 3).

Fig. 3

Boxplots of CpGs Comprising 6 DMRs of Radiation-Induced Severe Adverse Effects. Box plots display the differences in mean DNA methylation levels between cases (red) and controls (light blue). The y-axis represents the methylation rates

In the GO analysis, we identified 18 significantly enriched biological processes. The most prominent GO terms included intrinsic apoptosis, lipid and sugar metabolic processes, and the response to tumour necrosis factor (Fig. 4). Additionally, using the Reactome database, we found five pathways significantly enriched with the DMPs, all displaying an adjusted p-value < 0.05 (Fig. 4). Notably, the GO terms carbohydrate catabolic process, cellular response to tumour necrosis factor, monosaccharide metabolic process, and regulation of small molecule metabolic processes exhibited predominant hypomethylation, while the cellular response to molecule of bacterial origin and intrinsic apoptotic signalling pathway were predominantly hypermethylated (Fig. 5).

Fig. 4

Dot plot of the 18 gene ontology biological processes and the five significantly enriched pathways. Dot colours along the x-axis represent different FDR p-values, and dot size indicates the number of genes involved

Fig. 5

Density plot of methylation changes in genes with significantly enriched biological functions and pathways. Density curves illustrate the distribution of hypo- and hypermethylated genes in cases. A black dashed line marks zero mean methylation change, indicating no net methylation change

Assessment of biomarkers’ classification ability.

The MRSsites model, derived from the weighted sum of the 56 identified DMPs (Table 2), demonstrated a mean AUC of 0.87 (Fig. 6).

Table 2 Significant DMPs associated with late severe GI and GU adverse effects

Fig. 6

Mean ROC curves and AUC violin plots from cross-validation of three methylation classification models. The mean ROC curves A and AUC violin plots B derived from four-fold cross-validation with 2,000 iterations for three methylation classification models: MRSsites (red), MRSregions (light blue), and SVMsites (dark blue). The ROC curves depict the average classification performance of each model, while the violin plots show the distribution of AUC values across cross-validation iterations. The horizontal red dashed line in the plots indicates a zero-effect threshold

The MRSregions model, which integrated the 6 DMRs (Table 3) into a logistic regression framework, achieved a mean AUC of 0.89 (Fig. 6). After fourfold cross-validation with 2,000 iterations of the MRSregions model, the mean effect sizes and directions for each of the 6 DMRs were obtained (Fig. 7, Suppl. Figure 1). DMRs overlapping the ACAP1 (mean effect = 0.62, p-value = 0.004), IL-1RL2 (mean effect = 1.60, p-value = 0.020), NTM (mean effect = 1.20, p-value = 0.025), VOOP1 (mean effect = 1.29, p-value = 0.032), and an intergenic region (mean effect = 0.38, p-value = 0.026) showed significant positive directional effects. In contrast, the DMR within the AKR1E2 gene region displayed inconsistent directional effect across iterations (Fig. 7) and did not reach the significance threshold (mean effect = 0.70, p-value = 0.29).

Table 3 Significant DMRs Associated with Late Severe GI and GU Adverse Effects

Fig. 7

Violin plots showing weight distributions for six significant DMRs in the MRSregions logistic regression model. Each violin plot depicting the weight distributions from the logistic regression model for six significant differentially methylated regions (DMRs) forming the MRSregions model. Each plot represents the effect size and direction of the DMRs associated with the gene regions ACAP1, AKR1E2, IL1RL2, NTM, VOOP1, and an intergenic region (IR)

Using the ReliefF feature selection algorithm with a 100-repetition bootstrap approach, 8 CpG sites with mean importance weight > 0.11 were identified (Suppl. Table 1). These CpG sites also exhibited pairwise correlations r² < 0.6, as determined by Pearson’s correlation coefficient (Suppl. Figure 2). The SVMsites model developed by integrating these features exhibited superior performance, with a mean AUC of 0.98.

Stochastic epigenetic mutations

The cumulative analyses of SEMs did not yield significant results for either k = 2 or k = 3 (p-values 0.78, 0.84 respectively). Similarly, no significant association was found when regressing the phenotype on SEMs grouped by gene regions.

When examining the number of SEMs included within RT pathways, no significant association was observed with the phenotype. This was consistent for genes involved in DNA repair (k2 p-value = 0.95), cell senescence (k2 p-value = 0.95), and the immune system (k2 p-value = 0.64; k3 p-value = 0.343) (Table 4).

Table 4 SEM Regression Analysis Results with k Parameter Set to 2 or 3

Epigenetic clock

No significant differences were observed between the two groups when comparing chronological age (mean cases = 71.0, mean controls = 70.7, t = 0.23, p-value = 0.81), predicted methylation age (mean cases = 72.99, mean controls = 72.39, t = 0.53, p-value = 0.59) or acceleration age (mean cases = 0.21, mean controls = − 0.19, t = 0.56, p-value = 0.57). Additionally, no significant correlation was found between the total number of SEMs (log(SEM)) and patients’ chronological age (k2: r² = − 0.0536, p-value = 0.59 | k3: r² = − 0.0657, p-value = 0.51), epigenetic age (k2: r² = − 0.0030, p-value = 0.98 | k3: r² = 0.00278, p-value = 0.98) or acceleration age (k2: r² = 0.0569, p-value = 0.56 | k3: r² = 0.0797, p-value = 0.42).

Discussion

In this study, we conducted the first epigenome-wide association analysis in prostate cancer patients to identify methylation patterns that could stratify individuals based on their risk of developing toxicities following RT. We discovered 6 significant DMRs and 56 DMPs in peripheral blood, which serve as potential biomarkers for predicting severe long-term GI and GU adverse effects induced by RT. These epigenetic biomarkers are located in genes that provided insights into potential mechanisms underlying these radiation-induced adverse effects. Notably, we observed substantial changes in DNA methylation within genes associated with DNA repair and inflammatory pathways, as well as DMRs linked to tissue repair and oxidative stress response. Additionally, our findings suggest a potential role for Short Interspersed Nuclear Elements (SINEs) and Long Interspersed Nuclear Elements 1 (LINE-1) in radiation-induced adverse effects, warranting further investigation.

Five of the 6 identified DMRs are located within gene regions, while one is situated in an intergenic region. The gene neurotrimine (NTM), which is expressed in multiple tissues including the intestine and urinary bladder [37], promotes neurite outgrowth and cell adhesion through homophilic mechanisms during both normal development and smooth muscle repair [38]. The hypermethylated DMR identified in NTM may contribute to RT-induced GI and GU toxicities by impairing smooth muscle repair through the disruption of neurite outgrowth and cell adhesion.

Our analysis identified a 10-CpG hypomethylated DMR within the body of the ACAP1 gene, which is crucial for integrin β1 recycling. As integrin β1 plays a fundamental role in tissue repair and regeneration [39, 40], the reduced methylation of this DMR in ACAP1 may impair the recycling efficiency of integrin β1. This impairment could lead to defective repair and regeneration of tissues damaged by ionizing radiation, ultimately contributing to GI and GU adverse effects following RT.

A hypermethylated DMR was identified within IL-1RL2, a gene essential for the regulation of proinflammatory cytokine secretion. This differential methylation pattern may be associated with post-RT inflammation [41]. Additionally, we identified two hypermethylated DMRs linked to the oxidative stress response: one within the VOPP1 gene, where reduced expression leads to elevated ROS levels and induced apoptosis [42], and another in the promoter region of the AKR1E2 gene, which is involved in the detoxification of reactive compounds generated during oxidative stress [43].

These findings suggest that epigenetic alterations in genes regulating both inflammatory responses and oxidative stress may contribute to the development of radiation-induced toxicities.

The hypermethylated intergenic DMR on chromosome 13 overlaps with a conserved 70-bp region containing SINE and a LINE-1. Typically, CpG methylation suppresses the retrotransposition of LINE-1s and regulates SINEs retrotransposition [44], preventing exon or intron disruption and mis-splicing. Ionizing radiation can alter this pattern, potentially increasing retrotransposition [45]. While this observation suggests that RT could influence LINE-1 methylation and their expression, the extent of their genomic impact remains unclear. These findings underscore the need for further investigation into the potential role of SINEs and LINE-1s in the genomic response to RT.

Among the 56 identified DMPs, 43 are located within gene regions. Eight of these DMPs are in genes involved in processes relevant to RT response in non-tumoral tissues and their repair, such as apoptosis, fibrosis, epithelial-mesenchymal transition (EMT), inflammation, and DNA repair. Specifically, JHDM1D/KDM7A (cg08657377) modulates TGF-β-mediated transcription of RHOJ, thereby influencing angiogenesis and apoptosis through BCL2 [46,47,48]. FLCN (cg15055494) is involved in TGF-β1-mediated tissue repair and fibrosis [49,50,51,52]. Three DMPs are associated with genes related to EMT, which is essential for wound healing, tissue regeneration, and fibrosis. These include GNB1 (cg25666010) [53, 54], which impacts cell proliferation and migration, CNTN1 (cg26359940), which influences PI3K/Akt signalling [55, 56], and ROR2 (cg14185656), which affects proliferation, migration, and EMT in the context of tumour progression [57]. Additionally, two DMPs are linked to DNA repair and DNA damage response: ANKRD11 (cg08417142), which regulates chromatin accessibility and gene expression [58,59,60,61], and RNF114 (cg22657044), which is involved in the PARylation-dependent DNA damage response [62]. Lastly, CFD (cg02534744) exhibits potential anti-inflammatory effects by modulating IL-6 and IL-4 levels [63].

Additionally, 6 of these DMPs are associated with genes directly involved in the response to ionizing radiation: RNF40, SFN, RNASE6, FLI1, C14orf180, and ADARB2. RNF40, which is phosphorylated by ATM and ATR following ionizing radiation, works alongside RNF20 to facilitate H2B ubiquitination and chromatin reorganization, essential for effective DNA repair [64,65,66] and apoptosis regulation [67]. The FLI1 gene enhances RT resistance by improving DNA repair [68]. Mithramycin A (MithA), an EWS::FLI1 fusion protein inhibitor, has been proposed as a radiosensitizer [69]. The SFN gene encodes a cell cycle checkpoint protein crucial for radiation response [70]. In association studies, RNASE6 has been shown to predict radiation response and survival in glioblastoma [71], C14orf180 is included in a signature related to radiosensitivity in soft tissue sarcoma [72], and a SNP in ADARB2 is associated with radiation cytotoxicity in lymphoblastoid cell lines [73].

We used the identified DMPs and DMRs to develop prognostic prediction models for GI/GU post-RT toxicities. When analysing large-scale datasets (> 850 k CpG sites), the risk of overfitting becomes significant. To mitigate this issue, model performance was evaluated using fourfold cross-validation with 2,000 iterations.

The methylation risk score (MRSregions), derived from the logistic regression model incorporating the 6 DMRs achieved a mean AUC of 0.89, demonstrating strong classification potential [18]. Furthermore, the mean weight coefficients (Fig. 7) and mean p-values (Suppl. Figure 1) for each region within the model indicated that the genes NTM, ACAP1, IL-1RL2, and VOPP1, discussed above, exert the most substantial impact on the progression of GI/GU post-RT toxicities.

The logistic regression model incorporating all DMPs (MRSsites) achieved a mean AUC of 0.87, demonstrating the classification capability of the 56 CpG methylation signature. However, despite efforts to prevent overfitting, using 56 features to classify a cohort of 105 patients may have resulted in the identification of patterns specific to this study population, potentially limiting the model's generalisability.

To address this limitation, we employed the ReliefF feature selection algorithm with a 100-repetition bootstrap approach, identifying 8 key DMPs for building a Support Vector Machine model (SVMsites). This classifier achieved a remarkable mean AUC of 0.98, further highlighting the predictive potential of CpG methylation for post-RT toxicities. As the performance evaluation of the three models was conducted exclusively on the discovery cohort, the high classification accuracy observed may not necessarily translate to independent datasets. Without external validation in replication cohorts, the robustness and clinical utility of the proposed models remain uncertain.

Enrichment analysis of DMPs and DMRs revealed that the intrinsic apoptotic signalling pathway and p53 transcriptional regulation are the most significantly affected processes. These pathways are critical in cellular responses to radiation, causing direct damage like DNA strand breaks and mitochondrial dysfunction, as well as indirect damage through ROS [74]. Such damage can activate intrinsic apoptosis, either through p53 or independently, leading to pre-mitotic cell death [74,75,76]. In high-turnover tissues, like the GI epithelium, p53 promotes apoptosis by interacting with pro-apoptotic Bcl-2 family proteins at the mitochondrial membrane, inducing mitochondrial permeabilization [75, 76].

Certain CpG sites demonstrate consistent changes in DNA methylation with age, allowing for the estimation of biological or epigenetic age. However, our analysis did not reveal any association between these markers and radiation-induced late toxicities in prostate cancer. Additionally, SEM analysis did not identify any statistically significant associations in our study.

This study had several limitations, including the relatively small sample size and recruitment from a single hospital, which might affect the generalisability of our findings, particularly the resulting machine learning classification models. The high AUC values observed must be interpreted with caution, as they were derived from the discovery cohort and may reflect overfitting rather than true predictive power. Validation in independent cohorts from diverse medical centres is essential to assess the robustness and real-world applicability of these models. Additionally, while we utilized blood samples, which may not fully capture the epigenetic landscape of irradiated tissues, blood-based biomarkers remain valuable for non-invasive risk assessment and personalized treatment strategies. Future research should focus on confirming the identified methylation markers, validating the prediction models in larger and more diverse populations, and exploring the potential role of epigenetic modifiers in reducing radiation sensitivity in healthy tissues.

Conclusions

We conducted the first epigenomic study to assess the risk of RT-induced adverse effects in prostate cancer patients. We identified candidate methylation biomarkers and key genes associated with radiation-induced adverse effects, highlighting their involvement in critical biological processes such as fibrosis, inflammation, DNA repair, and tissue regeneration. Additionally, we identified an association potentially related to SINE and LINE-1 elements which may influence genome structure and function, though their precise role requires further investigation.

Our findings suggest that DNA methylation changes could predict RT-induced toxicity in prostate cancer. Incorporating these biomarkers into risk prediction models prior to treatment could facilitate more personalized RT strategies, potentially improving patient outcomes and quality-of-life by minimizing adverse effects. However, as the classification models were developed and evaluated within a single discovery cohort, their generalisability remains uncertain. External validation in larger, independent cohorts is essential to establish their robustness and clinical applicability.