Peripheral blood DNA methylation signatures predict response to vedolizumab and ustekinumab in adult patients with Crohn’s disease: The EPIC-CD study

Floris de Voogd for his assistance in evaluation of ultrasound images. 33 Conflicts of interest: The AmsterdamUMC has a patent pending for the vedolizumab and ustekinumab response 34 prediction models presented in this manuscript. ALY received honoraria from Janssen, Johnson & Johnson, 35 DeciBio and was employed by GSK. WJ received honoraria from Janssen, Johnson & Johnson and is a 36 cofounder of AIBiomics BV. GD received speaker fees from Janssen, Johnson & Johnson. EL is a cofounder of 37 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint NOTE: This preprint reports new research that has not been certified by peer review and should not be used to guide clinical practice. Horaizon BV and AIBiomics BV. The remaining authors disclose no conflicts. TC received honoraria from Janssen 38 and Takeda. The remaining authors have no conflicts of interest to declare. 39 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Graphical abstract 40 41 42 43 44 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint

45 Biological therapeutics are now widely used in Crohn’s disease (CD), with evidence of 46 efficacy from randomized trials and real-world experience. Primary non-response is a 47 common, poorly understood problem. We assessed blood methylation as a predictor of 48 response to vedolizumab (VDZ, anti-a4b7 integrin) or ustekinumab (USTE, anti-IL-12/23p40). 49 We report a two-center, prospective cohort study in which we profiled the peripheral blood 50 DNA methylome of 184 adult male and female CD patients prior to and during treatment with 51 VDZ or USTE in a discovery (n=126) and an external validation cohort (n=58). We defined 52 epigenetic biomarkers that were stable over time and associated with combined clinical and 53 endoscopic response to VDZ or USTE with an area under curve (AUC) of 0.87 and 0.89, 54 respectively. We validated these models in an external cohort yielding an AUC of 0.75 for 55 both VDZ and USTE. These data will now be prospectively tested in a multicenter 56 randomized clinical trial. 57 58 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint

59 Crohn’s disease (CD) is an incurable, chronic, relapsing inflammatory bowel disease (IBD) 60 caused by a complex interplay between the environment, gut microbiome, and a 61 dysregulated immune system in genetically susceptible patients1,2. Accessibility to high-62 throughput “omics” technology has enhanced our understanding of the underlying molecular 63 pathogenesis of CD leading to the development of several monoclonal antibodies or 64 biologicals that target specific inflammatory pathways in an effort to suppress the 65 inflammation and to induce and or maintain a state of clinical and endoscopic remission3. 66 Currently, the repertoire of approved biologicals in CD includes anti-TNF antibodies 67 (infliximab (IFX) and adalimumab (ADA)), the anti-α 4β 7 integrin antibody vedolizumab (VDZ) 68 and the anti-IL12/23p40 antibody ustekinumab (USTE). More recently, specific IL23p19 69 antibodies and JAK inhibitors have also been approved for clinical use4-6. Despite the 70 established efficacy of these biological treatments to induce corticosteroid-free clinical 71 remission in up to 65% of CD patients, sustained endoscopic remission is observed in not 72 more than a third of patients after 1 year of treatment7-9. This creates a clinical challenge 73 since therapeutic guidelines suggest to use endoscopic remission as a target. An 74 increasingly common clinical scenario is the choice between VDZ or USTE as second-line 75 treatment for patients who have not responded to anti-TNF therapies. 76 To date, treatment selection has been based on a trial-and-error approach. Given the limited 77 efficacy, many patients are therefore treated with insufficiently effective treatment, which is 78 associated with an increased risk of complications (stenosis, fistula, abscesses, nutritional 79 deficiencies) and surgery. The development of strategies that allow selection of treatment 80 based on the likelihood of response is an important unmet need. Although various efforts 81 using clinical10, transcriptomic11-13, proteomic14 or microbial15,16 technologies have been 82 investigated and reported, no predictive biomarkers have made their way to clinical 83 application17-19. 84 DNA methylation is one of the most studied epigenetic features characterized by the covalent 85 binding of methyl groups to nucleotides, most often a cytosine in a cytosine-phosphate-86 guanine (CpG) sequence in humans20. DNA methylation is believed to play an essential role 87 in the regulation of gene expression, thereby determining cellular phenotype and behavior 88 without altering the DNA sequence itself21,22. Within the context of IBD, DNA methylation has 89 gained particular interest due to its dynamic interaction with the environment and suggested 90 epigenetic-microbial crosstalk23-25. A number of recent studies demonstrated differential DNA 91 methylation profiles associated with the presence of CD and/or specific CD-phenotypes26,27 in 92 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint peripheral blood leukocytes (PBL)28-30, intestinal mucosa31 and specific cell-types thereof32,33. 93 Most studies proposed a potential role of the DNA methylome in diagnostics and prediction 94 of treatment response. 95 In the current study we performed an epigenome-wide association study (EWAS) in which 96 we identified and validated prognostic DNA methylation signatures in peripheral blood of 97 adult CD patients associated with objective therapeutic response to VDZ and USTE. We 98 subsequently interrogated whether these signatures were stable over time and assessed the 99 influence of common confounding variables. We further performed additional validation 100 analyses of the model against patients with previous non-response to one or both drugs. 101 Finally, we investigated the association of the identified DNA methylation biomarkers with 102 gene expression through transcriptomic analyses. 103 104 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint

105 Study population 106 We prospectively recruited a discovery cohort of 126 adult patients at the IBD Center of 107 Amsterdam University Medical Centers (Amsterdam UMC), Amsterdam, Netherlands. All 108 patients had active symptomatic and endoscopic CD and were scheduled to start VDZ (N = 109 64) or USTE (N = 62). Evidence for active disease was documented with a validated clinical 110 score (Harvey Bradshaw Index (HBI), median 8 (interquartile range (IQR) 4-12)), biochemical 111 tests (serum C-reactive protein (CRP), median 6.1 mg/L (IQR 2.2-14.8) and fecal calprotectin 112 (FCP), median 903 µg/g (IQR 278-1816)) and also endoscopic signs of inflammation 113 measured with the simple endoscopic score for CD (SES-CD), median 9 (IQR 6-15). In 114 addition, an external validation cohort of 58 adult CD patients starting VDZ (N = 25) or USTE 115 (N = 33) biological therapy were recruited at the John Radcliffe Hospital, Oxford, United 116 Kingdom. Peripheral blood leukocyte (PBL) samples were obtained prior to treatment 117 initiation and upon response assessment at a median of 27 (IQR 20-33) weeks into 118 treatment, when patients were classified as responder (R) or non-responder (NR). A detailed 119 overview of all the clinical characteristics across the different cohorts and treatments can be 120 found in the methods section as well as in Table 1. 121 Blood DNA methylation profiling predicts response to vedolizumab and 122 ustekinumab 123 To identify prognostic biomarkers of VDZ and USTE response, we performed supervised 124 machine learning through stability selected gradient boosting34,35 on blood samples obtained 125 shortly before the start of treatment (Fig. 1a). The model was trained on the discovery cohort 126 acquired in Amsterdam and subsequently validated in the validation cohort acquired in 127 Oxford. We were able to generate response-predicting models with an area under the curve 128 (AUC) of 0.87 with a standard deviation of 0.06 and 0.89 with a standard deviation of 0.08 for 129 VDZ and USTE, respectively, when testing against the discovery cohort (Fig. 1b). The 130 models comprised 25 and 68 differentially methylated CpGs for VDZ and USTE, respectively 131 (Fig. 1c, Supplementary Fig. 1-2 and Supplementary Tables 1-2). Validating our models 132 against the independent validation cohort yielded an AUC of 0.75 for both VDZ and USTE 133 (Fig. 1b), indicating reproducible response-associated differences in DNA methylation prior 134 to the start of treatment. 135 While the discrepancy between discovery and validation performances is expected as part of 136 machine learning, we hypothesized that the differences were partially due to differences in 137 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint response assessments. During the COVID-19 pandemic, endoscopic assessments were 138 scheduled less frequently resulting in the validation cohort consisting of patient samples in 139 whom response was defined on both strict (combined clinical-, biochemical- and endoscopic 140 evaluation) or modified criteria (combined clinical- and biochemical evaluation). We therefore 141 investigated whether a difference in performance could be observed between both methods 142 of evaluation. This stratification indicated that accuracy was optimized by using the strictly 143 defined combination of clinical- and endoscopic endpoints for both VDZ (AUCstrict = 0.83 144 versus AUCmodified = 0.66) and USTE (AUCstrict = 0.83 versus AUCmodified = 0.72) 145 (Supplementary Fig. 3a). Application of the strict criteria yielded performances similar to the 146 discovery dataset thereby confirming that our model is likely more capable at predicting 147 response as defined using combined clinical-, biochemical- and endoscopic evaluations. 148 In addition, we were interested whether the performance of our models was affected by prior 149 exposure to anti-TNF medication. We observed a better performance of our models among 150 anti-TNF naïve compared to anti-TNF exposed patients for both VDZ (AUCnon-exposed = 0.85 151 versus AUCexposed = 0.66) and USTE (AUCnon-exposed = 0.97 versus AUCexposed = 0.63) 152 (Supplementary Fig. 3b). 153 Focusing on the practical implications in a clinical setting, we calculated a sensitivity of 0.769 154 and a specificity of 0.67 for VDZ and both a sensitivity and specificity of 0.73 for USTE 155 (Table 2). Next, we computed the likelihood ratio of response and the post-test probability of 156 response to aid clinicians in accurately predicting response following a positive test outcome. 157 From Lowenberg et al. 7, we note that 50 of the 110 VDZ-treated CD patients presented with 158 endoscopic response at week 52 indicating a pre-test probability of response of 0.45. At the 159 calculated sensitivity of 0.77 and a specificity of 0.67, the likelihood ratio of response is 2.31, 160 thereby making the post-test probability of response 0.65 for VDZ. Similarly, for USTE, prior 161 research showed that 75 of the 179 USTE-treated CD patients presented with endoscopic 162 response at week 528, indicating a pre-test probability of response of 0.42. At the calculated 163 sensitivity and specificity of 0.73 the likelihood ratio is 2.67 and the post-test probability of 164 response is 0.66. Taken together, the probability that a patient would actually respond to 165 VDZ or USTE when classified as responder is 0.20 and 0.24 higher than the current standard 166 of care. 167 Methylation status of predictor CpGs remains stable over time 168 As baseline samples were acquired pre-treatment, we sought to understand whether the 169 initiation of treatment affected the methylation status of the predictor CpGs. To this end, we 170 compared samples obtained at the time of response assessment (T2) with samples obtained 171 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint pre-treatment (T1). We could not identify any statistically significant differences in DNA 172 methylation levels for any of the predictor CpGs (Fig. 2a). Indeed, comparing the differences 173 over time suggested that the mean difference between R and NR was similar both pre-174 treatment and at response assessment (Fig. 2b). Furthermore, a two-way, mixed, 175 consistency intra-class correlation (ICC) analysis indicated highly-stable DNA methylation 176 over time with 24 out of 25 VDZ and 62 out of 68 USTE predictor CpGs presenting ICC 177 values ≥ 0.75 (Fig. 2c, Supplementary Fig. 1 and 2). This observation was corroborated by 178 interrogating our previous longitudinal consistency analysis of peripheral blood DNA 179 methylation from 46 adult IBD patients collected at 2 time points with a median of 7 years 180 (range, 2-9 years) in between36. Here, we observed that the majority (16 out of 25 VDZ and 181 52 out of 68 USTE) of the predictor CpGs presented “good” (0.75 ≤ ICC < 0.9) to “excellent” 182 (0.9 ≥ ICC) stability37 over a median span of 7 years (Fig. 2c). As a final validation, we 183 utilized the prognostic model to predict response to therapy of the samples obtained at 184 response assessment, where we obtained better performances compared to the pre-185 treatment samples (AUCVDZ = 0.97; AUCUSTE = 1.00) (Fig. 2d). Altogether, our observations 186 suggest that response-associated differences in DNA methylation detected prior to treatment 187 remain stable during treatment. 188 Patients who previously experienced treatment failure to multiple 189 biological treatments are classified correctly 190 Having established that the prognostic response prediction models for both drugs perform 191 well both pre- and into treatment, we explored the performance of our model in a further 192 independent cohort of 33 adult CD patients from whom blood had been collected after 193 sequential treatment with anti-TNF, VDZ, and/or USTE. Treatment failure was determined 194 with endoscopy in 28 (82%), MRI in 2 (6%), or the need for surgical resection in 4 (12%), 195 combined with biochemical biomarkers such as CRP in 12 (35%) and/or FCP in 16 (47%) 196 and clinical parameters in 32 (94%) patients. In this cohort, a relatively large proportion of 197 patients had extensive disease location (67.6%), perianal disease (52.9%), and/or IBD-198 related surgery (70.6%) in the past (Supplementary Table 3). 199 Both VDZ and USTE response prediction models accurately predicted NR, with 24 out of 27 200 (88.9%) VDZ-NR patients and 24 out of 26 (92.3%) USTE-NR patients correctly identified. 201 Focusing specifically on the patients with previous non-response to both VDZ and USTE the 202 models correctly classified 16 out of 19 (84%) as NR. Notably, 12 of out of these 19 patients 203 (75%) were anti-TNF experienced. 204 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Assessment of potential confounding variables 205 Through multiple linear regression analyses we observed that 22 of the 25 (88%) VDZ 206 response-associated CpGs and 38 of the 68 (55%) USTE response-associated CpGs (Fig. 207 3a and Supplementary Fig. 1-2) presented statistically significant differences. The large 208 discrepancy between the linear regression analyses and the USTE response-associated 209 CpGs indicates that a more complex non-linear relationship exists among the response-210 associated predictor CpGs and underscores that statistical p-values may not equal biologic 211 or functional relevance and/or importance at an individual level. As it has been established 212 that the peripheral blood DNA methylome is associated with certain phenotypic 213 characteristics such as sex, age, smoking status, as well as the underlying cellular 214 composition38-41, we investigated whether any of these variables confounded our results. In 215 doing so with linear regression analysis, we observed that 12 out of 22 (55%) and 30 out of 216 38 (79%) markers remained significantly associated with response for VDZ and USTE, 217 respectively (Fig. 3a). In terms of effect size, the mean percentage methylation difference 218 between R and NR of the predictor CpGs on average decreased by 15% and increased by 219 5% in VDZ and USTE, respectively (Fig. 3b). To understand whether the confounding 220 variables are capable of predicting response, we constructed a prediction model solely based 221 on the confounding variables using the discovery cohort and tested this on the validation 222 cohort. The confounder model yielded an AUC of 0.57 and 0.65 for VDZ and USTE, 223 respectively, against the validation cohort indicating worse performance than the prediction 224 model based on CpGs only (Fig. 3c). Our results indicate that while there appears to be a 225 significant association with sex, age, smoker status and blood cell distribution, the CpG-only 226 prediction model outperforms a confounder-only model. 227 We next investigated whether the predictor CpGs were significantly associated with severity 228 of systemic and intestinal inflammation at baseline measured using CRP and FCP, as was 229 previously reported by Somineni et al.42. For VDZ, 5 predictor CpGs significantly associated 230 with CRP whereas only a single CpG was associated with FCP (Supplementary Fig. 4). For 231 USTE, we observed 9 predictor CpGs associated with CRP and 2 predictor CpGs with FCP 232 (Supplementary Fig. 4). In all cases, mean differences in methylation were smaller than 233 0.5%. These observations suggest that the majority of the identified predictor CpGs are 234 independent of inflammatory status at onset of therapy. 235 RFPL2 presents concordant response-associated differential methylation 236 and expression 237 To understand the functional relevance of the predictor CpGs, we annotated the them to their 238 respective genes based on whether they were located in either a gene promoter or enhancer 239 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint resulting in the annotation of 20 VDZ response-associated predictor CpGs to 16 unique 240 genes (Supplementary Table 1) and 43 USTE response-associated predictor CpGs to 46 241 genes (Supplementary Table 2). We performed transcriptomic analyses on a subset of 242 samples from the discovery cohort (VDZ; NT1=10R10NR, NT2=10R8NR and USTE; NT1=14R11NR, 243 NT2=10R9NR). Comparing R with NR identified pretreatment differences for predictor CpG-244 associated genes TULP4 (p-valueT1 = 4.38E-02) and RFPL2 (p-valueT1 = 4.72E-02) for VDZ 245 (Fig. 4a, b and Supplementary Table 4), with RFPL2 presenting a significant inverse 246 correlation between DNA methylation and gene expression (Pearson r = -0.65; p-value = 247 1.39E-03) (Fig. 4c). For USTE, we observed significant differences in the expression of 248 predictor CpG-associated genes MRC1 (p-valueT1 = 4.25E-04) and TMEM191B (p-valueT1 = 249 7.31E-03) (Fig. 4a and d and Supplementary Table 5) with TMEM191B presenting a 250 significant positive correlation between DNA methylation and gene expression (Pearson r = 251 0.57; p-value = 3.11E-03) (Fig. 4e). We next investigated whether predictor CpG-associated 252 genes presented differential expression after the start of treatment by comparing R with NR 253 at response assessment. VDZ predictor CpG-associated genes MCM2 (p-valueT2 = 2.40E-254 03) and RFPL2 (p-valueT2 = 3.88E-03) were differentially expressed at response assessment 255 (Fig. 4a, f and Supplementary Table 4) with RFPL2, once again, presenting a significant 256 inverse correlation between DNA methylation and expression (Pearson r = -0.55; p-value = 257 0.017) (Fig. 4g). For USTE, predictor CpG-associated genes POTEF (p-valueT2 = 1.47E-02), 258 HDAC4 (p-valueT2 = 2.18E-02), PARP4 (p-valueT2 = 3.49E-02) and MARK3 (p-valueT2 = 259 2.87E-02) presented differential expression at response assessment (Fig. 4a, h and 260 Supplementary Table 5), but did not show any significant correlation with DNA methylation. 261 Taken together, our results show that some of the response predictor CpG-associated genes 262 present differential expression either pretreatment and/or during response assessment, with 263 VDZ- and USTE- response-associated genes RFPL2 and TMEM191B, respectively, 264 presenting a significant correlation between DNA methylation and gene expression. 265 Nonetheless, we acknowledge that most predictor CpG-associated genes present no 266 response-associated differential expression. 267 268 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint

269 Although the introduction of biologicals has transformed the care of patients with CD, the 270 current clinical practice of treatment selection remains suboptimal. Nonetheless, a large body 271 of real-world evidence studies suggests there could be a place for a more individualized 272 approach43,44. Therefore, in parallel with ongoing drug development, predictive biomarkers 273 that allow for selection of successful medical therapy would represent a major step forward in 274 clinical care. 275 Here, we conducted a longitudinal case-control study where we identified methylation 276 signatures composed of 25 and 68 markers associated with combined endoscopic, 277 biochemical and clinical response to VDZ and USTE, respectively, in a cohort of adult CD 278 patients in whom treatment response was assessed with stringent criteria. We were able to 279 build models with significant predictive performance at an AUC > 0.85 for both models. The 280 models demonstrated similar performance in an independent, external validation with an 281 AUC of 0.75 for both models. 282 Recent real-world data and (post-hoc) findings from both the GEMINI and UNITI trials 283 indicate superior response to both VDZ and USTE in anti-TNF naïve patients45-49. In our 284 discovery cohorts, 77% of VDZ and 98% of the USTE-treated patients were previously 285 exposed to anti-TNF medication. Stratifying the patients in the validation cohort by previous 286 anti-TNF exposure showed that both models performed noticeably better in anti-TNF naïve 287 rather than exposed patients. However, the number of patients included in both subset 288 comparisons are relatively small and further exploration using larger groups patients are 289 needed. Furthermore, we demonstrate the ability of our models to effectively identify patients 290 with previous non-response to both VDZ and USTE treatment. This holds true for both anti-291 TNF naïve and experienced patients, providing significant importance for clinical practice, as 292 the precise prediction of non-response to both drugs offers clinicians the opportunity to make 293 informed decisions. For anti-TNF experienced patients, this could mean a direct switch 294 towards newer modes-of-action (i.e. JAK-inhibitors). 295 While this study is the first to demonstrate the utility of DNA methylation profiling in whole 296 blood for objective response to VDZ and USTE in CD patients, previous analyses from two 297 separate studies explored its application to predict anti-TNF response in both CD and UC 298 patients50,51. In the first study, the authors sought to identify an anti-TNF response-associated 299 profile combining integrated methylation and gene expression data from samples taken 300 before and 2 weeks into treatment using a primary endpoint of clinical remission at week 301 1450. The observations were made using a mixed cohort of 37 IBD (18 CD, 19 UC) patients 302 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint and validated against a publicly available gene expression data of 20 CD patients. In the 303 second study, the authors did not show replication of these observations using methylation 304 data of 385 patients, as part of the previously published PANTS study51,52. In this study, bio-305 naïve patients with active CD started treatment with adalimumab or infliximab. Interestingly, 306 the authors report 323 differentially methylated positions annotated to 210 genes identified at 307 baseline that were significantly associated with serum drug concentrations at week 14. This 308 profile could potentially be of interest to identify patients in need of intensified anti-TNF drug 309 monitoring or dosing. The methodological differences and lack of endoscopic data in the 310 previous experiments preclude direct comparison with our study in which outcome 311 assessments were more stringent. 312 Through two separate stability analyses, we demonstrated both short- and long-term hyper 313 stability of the majority of our identified CpG markers indicating their independence of 314 treatment as well as the resultant difference in inflammation. The latter is further evidenced 315 by the lack of correlation between the methylation status of the predictor CpGs and both 316 baseline CRP and FCP, therapy switch or even CD-related surgery53 suggesting that the 317 CpGs are response-predictors but do not directly mediate inflammation. Nonetheless, 318 confounder analysis indicated potential confounding for some of the predictor CpGs for age, 319 sex, smoking status and the blood cell distribution. We observed however that prediction 320 modeling using the confounders only performed significantly worse than the CpG model, 321 indicating that confounders do not contribute substantially to the predictive performance. 322 We note that the predictor CpGs annotate to genes associated with MHC class I (HLA-C) 323 and cell migration (TSPEAR, NID2)54,55 for VDZ and macrophage function (RHOJ, MARK3 324 and PCGF3)56-61 and polarization (PKNOX1 and MRC1)62-64, histone remodeling and Th17-325 differentiation (HDAC4)65-69 and TGF-β signaling (SMAD1 and EBF3)70-77 for USTE. 326 Integrating DNA methylation with their transcriptomic data indicates that VDZ-response-327 associated RFPL2, which encodes an E3 ubiquitin ligase78,79, and USTE-response-328 associated TMEM191B, which encodes a transmembrane protein, present concordant 329 differential methylation and expression. However, their exact role in either IBD and/or 330 response to therapy remains unknown to date. Accordingly, beyond the utility of the predictor 331 CpGs in classifying response to therapy, identifying their role in the pathogenesis and 332 etiology of non-response remains challenging at present and hence a subject for future 333 studies. 334 Our study derives its main strength from the sampling and endpoint assessment strategy. In 335 addition, patients with anti-drug antibodies or without a measurable serum drug 336 concentration and those that stopped treatment due to adverse events were excluded prior to 337 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint the selection of this cohort. Non-responders therefore reflect a more homogenous group of 338 true biological (i.e. pharmacodynamic) non-response rather than failure due to 339 pharmacokinetics or intolerance. Second, the relatively small sample size was justified by 340 performing a sample size estimation based on prior pilot experiments in combination with 341 previous studies on statistical power in EWAS80. Lastly, besides stability of the observed 342 methylation differences during induction- and maintenance treatment, our markers 343 demonstrated stable differences between R and NR over time, which was further 344 corroborated when interrogating our 7-year longitudinal DNA methylation survey. This time-345 independent behavior of the predictor CpGs indicates that exposure over time and change in 346 inflammatory parameters does not affect the response-associated behavior of the CpGs, 347 suggesting that the CpGs are very stable, which in turn increases its utility in clinical practice. 348 There are however some limitations to this study. First, a post-hoc power analysis of the 349 mean effect size difference between responders and non-responders indicated a mean 350 absolute difference of 9.9% and 7.6% for VDZ and USTE, respectively. At the collected 351 samples sizes, this would translate to a statistical power of approximately 98.4% and 64.3% 352 for VDZ and USTE, respectively80. Despite the lower power for USTE, we note that the model 353 remained performant in predicting response in the external cohort. Second, in the external 354 cohort response assessment was less stringent in approximately 70% of the patients due to 355 the COVID-19 pandemic during which we encountered a notable reduction in the access to 356 non-essential endoscopies, particularly in the UK81. Nonetheless, the modified response 357 criteria are a reflection of clinical parameters that physicians commonly use in daily practice. 358 This pragmatic approach enhances the overall generalizability of our results to a much larger 359 IBD population. Notably, the analysis where we included the available endoscopic outcomes 360 in the subset of UK patients for whom these data were available, enhanced the performance 361 to an AUC of 0.83 for both the VDZ and USTE models, reinforcing the validity of both models 362 in identifying objective responders to these biological therapies. While we acknowledge this 363 limitation, we believe that our comprehensive approach provides valuable insights into the 364 predictive capabilities of the models under diverse clinical scenarios. Third, while we 365 purposely used PBL samples as these are minimally invasive and easily obtained during 366 daily clinical practice, PBL represents a mixed cellular population. Therefore, the specific cell 367 types responsible for the observed predictive signal remain unidentified82. It should be noted 368 however that after correcting for the blood cell distribution, several predictor CpGs remained 369 statistically significant, indicating independence of cellular composition. Third, the majority of 370 the predictor CpG loci identified are situated within gene introns, complicating the biological 371 interpretation of our findings. Although the identified predictor CpGs collectively serve as a 372 strong predictor of response, we acknowledge that the underlying biology behind this 373 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint observation is more complicated than merely the inverse correlation between DNA 374 methylation and gene expression. Lastly, while strict removal of most catalogued and 375 predicted genetic variant-binding probes, we acknowledge that there is still a possibility that 376 underlying genetic differences could have influenced our outcome83. Nonetheless both 377 models performed effectively in both the discovery- (Dutch) and the validation (UK) cohorts. 378 In summary, our findings pave the way towards personalized medicine for CD. Several US 379 and European studies report a significant reduction in pharmaco-economic burden of CD if 380 clinical and endoscopic remission can be achieved with adequate treatment compared to the 381 cost of treatment of IBD patients on suboptimal medication84-88. In the absence of a predictive 382 biomarker panel, current endoscopic response rates at week 52 have been reported around 383 45% in VDZ and 42% in USTE treated CD patients7,8. Taken together, our biomarker panels 384 could potentially increase these proportions by approximately 20% for VDZ and 24% for 385 USTE, thereby significantly impacting healthcare costs and disease burden in these patients. 386 We acknowledge that clinical validation of our findings in a randomized prospective trial, 387 comparing our method of pre-treatment selection with current clinical practice, is needed to 388 firmly demonstrate both clinical and economic benefit89. To this end, the Omicrohn trial as 389 part of the ongoing Horizon Europe funded METHYLOMIC project has been launched and is 390 currently underway90. 391 392 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Online Methods 393 Study population and design 394 We prospectively recruited adult male and female CD patients that presented with a 395 combination of clinical, biochemical and endoscopic disease activity at ileo-colonoscopy and 396 were scheduled to start VDZ or USTE treatment within 1 year after the last endoscopy at the 397 Amsterdam University Medical Centers, University of Amsterdam, Amsterdam, Netherlands 398 (discovery cohort) and the John Radcliffe Hospital, Oxford, United Kingdom (validation 399 cohort). All patients were naïve to the biological of interest. 400 Patients were treated according to standard-of-care protocols, which for VDZ meant that 401 patients were given 300 mg infusions at week 0, 2 and 6 followed by infusions at an 8 week 402 interval. For USTE, standard-of-care involved patients receiving a single intravenous infusion 403 (6 mg/kg rounded to 260 mg, 390 mg or 520 mg) at week 0 and subsequent 90 mg 404 subcutaneous injections at an 8 week interval. For both VDZ and USTE, interval 405 intensification every 6 or 4 weeks with an additional week 10 infusion for VDZ or extra 406 intravenous boost infusion for USTE at the treating physicians’ discretion. To ensure 407 assessment of mechanistic and not pharmacokinetic failures of each biological treatment, 408 only patients with measurable serum concentrations without anti-drug antibodies at response 409 assessment were used for methylation analyses. 410 The project was approved by the medical ethics committee of the Academic Medical Hospital 411 (METC NL57944.018.16 and NL53989.018.15) and written informed consent was obtained 412 from all subjects prior to sampling, as well as by the National Health Service Research Ethics 413 committee. (REC reference: 21/PR/0010 Protocol number: 14833 IRAS project ID: 266041). 414 UK patients were recruited and consented under the ethics of the Translational 415 Gastrointestinal Unit biobank IBD cohort ethics (09/H1204/30) and GI cohort ethics 416 (16/YH/0247 and 21/YH/0206). The quality of the collected data and study procedures were 417 assessed by an independent monitor. 418 Vedolizumab discovery and validation cohort characteristics 419 The VDZ discovery cohort consisted of 64 patients (NR = 36, NNR = 28) of which 49 (77%) had 420 previously been exposed to anti-TNF treatment and 9 (14%) to USTE. Response was 421 defined on endoscopic as well as clinical/biochemical criteria in this group (see section 422 “Definitions of response”). In 62 (97%) patients, a follow-up endoscopy was performed. The 423 remaining 2 patients were classified as non-responders as a result of urgent surgical 424 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint intervention due to worsening of disease without CD-associated complications, such as 425 stenosis or perforating disease. R and NR presented overall comparable clinical 426 characteristics, with no significant differences in age, sex and smoking behavior. Importantly, 427 serum VDZ concentrations at T2 were not significantly different between R and NR (median 428 15 (IQR 7.7-20.5) versus median 14 (IQR 3.6-27.5), p-value = 0.77) although more NR 429 patients received an additional VDZ infusion at week 10 (35.7% vs 13.9%, p-value = 0.04). 430 Notably, patients in the NR group had more frequently been exposed to anti-TNF treatment 431 (89.3% vs 66.7%, p-value = 0.03) and/or USTE (25% vs 5.6%, p-value = 0.03) compared to 432 R group. Thirteen of the 15 R patients received VDZ as first-line biologic, had a shorter 433 disease duration, lower rates of previous surgery and perianal disease as well as a higher 434 percentage of B1 phenotype compared to anti-TNF experienced patients. 435 The VDZ validation cohort consisted of 25 additional patients (NR = 14 and NNR = 11) from 436 Oxford with a median 9 (IQR 3-15) year disease duration. Seventeen (68%) did not undergo 437 a follow-up colonoscopy to assess endoscopic response due to restrictions during the 438 COVID-19 pandemic and therefore were assessed using the modified definition of response, 439 which was defined as a combination of clinical and biochemical parameters (see “Definitions 440 of response” below). Out of these 25 patients, 13 (52%) were biological naïve, 12 (48%) 441 were anti-TNF-experienced and 7 (28%) were previously treated with USTE, which was 442 enriched among the NR group (54.5% vs 7.1%, p-value = 0.01). All other clinical 443 characteristics were comparable in R and NR, including age, sex, and smoking behavior. As 444 in the VDZ discovery cohort, the majority (69%) of biological naïve patients were responders 445 to VDZ. 446 Ustekinumab discovery and validation cohort characteristics 447 The USTE discovery cohort consisted of 62 patients (NR=30, NNR=32) of which 60 (98%) 448 were previously exposed to anti-TNF and 26 (42%) to VDZ. A follow-up endoscopy was 449 performed in 58 (94%) patients. The remaining 4 were evaluated with serial intestinal 450 ultrasound using validated criteria assessed by an expert IBD ultrasonographist. Clinical 451 characteristics did not significantly differ between R and NR, although the R population 452 consisted of more female patients (R = 80%, NR = 56.3%, p-value = 0.05). No significant 453 differences in treatment intensification (21.9% vs 6.7%, p-value = 0.08), extra intravenous 454 boost infusions (p-value = 0.26) or serum USTE concentrations at T2 between R and NR 455 were observed (median 3.0 (IQR 1.8-5.5) versus median 4.8 (IQR 2.1-8.6), p-value = 0.31). 456 The USTE validation cohort consisted of 33 patients (NR=22 and NNR=11), with a median 457 disease duration of 11 (IQR 3-20) years. Twenty-five (76%) did not undergo follow-up 458 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint colonoscopy as a result of pandemic-related restrictions on non-essential endoscopy and 459 therefore were assessed using the modified definition of response. Of the 33 included 460 patients, 9 were biological naïve (52%), 21 (63.6%) were anti-TNF-experienced and 5 461 (15.2%) were previously treated with VDZ. Responders presented a significantly longer 462 disease duration compared to non-responders (median 15 vs 5 years, p=0.05). No significant 463 differences in age, sex, and smoking behavior were observed. 464 Sample collection and storage protocols 465 In all patients, whole peripheral blood leukocyte (PBL) samples were collected for 466 measurement of epigenome-wide DNA methylation prior to the start of VDZ/USTE, before 467 the baseline endoscopy or the first infusion (time point 1) and after an interval of 6-9 months 468 into treatment (time point 2) using 4.0-6.0mL BD ethylenediaminetetraacetic acid (EDTA) 469 vacutainer tubes. For the Amsterdam discovery cohort, samples were subsequently 470 aliquoted into 1.10mL micronic tubes before storing at -80 ºC until further handling. The 471 Oxford validation samples were directly frozen at -80 ºC in preparation for later extraction. At 472 both time points, additional PBL samples were stored for the purpose of targeted gene 473 expression analyses from a subset of the Amsterdam discovery patients using 2.0-mL 474 PAXgene Blood RNA tubes and were frozen at –20 ºC for 24 h before storing at –80 ºC until 475 further handling. 476 Definitions of response 477 At response assessment, patients were classified as responders (R) or non-responders (NR) 478 based on a strict combination of endoscopic, biochemical and clinical criteria: ≥ 50% 479 reduction in the endoscopic SES-CD score, corticosteroid-free clinical remission (≥ 3 point 480 drop91 in HBI or HBI ≤ 4 and no systemic steroids) and/or biochemical response (CRP 481 reduction ≥ 50% or CRP≤ 5 mg/L and FCP reduction ≥ 50% or FCP ≤ 250 µg/g). Modified 482 response was defined as a combination of corticosteroid-free clinical- (HBI ≤ 4) and 483 biochemical (CRP ≤ 5 mg/L and/or FCP ≤ 250 µg/g) remission between week 26-52 without 484 treatment change through week 52. 485 DNA isolation and in vitro DNA methylation analysis 486 For the Amsterdam discovery cohorts, genomic DNA was extracted using the QIAsymphony. 487 We next assessed the quantity of DNA using the FLUOstar OMEGA and quality of the high-488 molecular weight DNA on a 0.8% agarose gel. Following these steps, 750ng of DNA per 489 sample was randomized per plate, to limit batch effects, after which genomic DNA was 490 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint bisulfite converted using the Zymo EZ DNA Methylation kit and analyzed on the Illumina 491 HumanMethylation EPIC BeadChip array. Aforementioned work was performed at the Core 492 Facility Genomics, Amsterdam UMC, Amsterdam, the Netherlands. 493 For the Oxford validation cohorts, genomic DNA was extracted using the Qiagen Puregene 494 Blood Core Kit C at the Oxford Translational Gastroenterology Unit, University of Oxford. 495 DNA samples were assessed for quality using the NanoDrop spectrometer (NanoDrop 1000, 496 Thermo Scientific). A total of 750ng of DNA per sample was randomized per plate, to limit 497 batch effects, after which genomic DNA was bisulfite converted using the Zymo EZ DNA 498 Methylation kit and analyzed on the Illumina HumanMethylation EPIC BeadChip array at 499 UCL Genomics, University College London, London, United Kingdom. 500 Raw methylation data pre-processing 501 DNA methylation data processing was orchestrated in Snakemake (v7.14.1)92. Raw 502 methylation data was imported into the R statistical environment (v4.3.1) using the 503 Bioconductor minfi93,94 package (v1.44). Quality control was performed using shinyMethyl 504 (v1.38.0)95 for probe level quality control, ewastools (v1.7.2)96 to ensure paired samples were 505 correctly labeled, and the Horvath clock97 to ensure a proper match with the metadata. One 506 patient sample from the VDZ and one patient sample from the USTE discovery cohort were 507 removed due to a discrepancy in the predicted sex with the annotated sex. Raw signals were 508 normalized using functional normalization98. Probes were annotated to their gene of interest 509 using the provided Illumina annotations, which were further enhanced with enhancer data as 510 obtained using the publicly accessible promoter-capture Hi-C data99. We subsequently 511 calculated the methylation signal in the form of percentage methylation. Technical artefacts 512 as a result of batch, plate and plate position were removed using ComBat (v0.9.7)100 as 513 implemented in the sva (v3.50.0)101 package, a tool for removing known batch effects in 514 microarray data using the parametric empirical Bayes framework, using the default 515 parameters. Probes hybridizing to allosomes were removed to identify sex-independent 516 differences. Moreover, probes hybridizing to known and tentative genetic variants were 517 removed to identify true methylation signals. Known genetic variants were identified based 518 on their presence in dbSNP (v137), whereas tentative genetic variants were based on the 519 methylation signal displaying a tri- or bimodal distribution, a hallmark of genetic variants 520 underlying DNA methylation102, as determined using gaphunter103 set to a threshold of 0.25. 521 The eventual number of CpGs used for training the prediction models were 806,308 and 522 808,815 for VDZ and USTE, respectively. For the validation of the prediction models, raw 523 methylation data from the validation cohort was pre-processed together with the discovery 524 cohort using functional normalization and ComBat to mitigate batch effects introduced by the 525 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint different experimental setup. From the combined dataset, the predictor CpGs were extracted 526 and the prediction model was recalibrated against the discovery dataset where after 527 predictions were made on the validation dataset. 528 Machine learning models: stability selected gradient boosting analyses 529 The machine learning modeling was divided into two steps, namely feature selection and 530 validation (Fig. 1a). Feature selection was performed on the discovery cohort, collected at 531 the AmsterdamUMC, whereas validation was performed on the validation cohort, collected at 532 the John Radcliffe Hospital. During the feature selection procedure, the model was at no 533 point exposed to the validation cohort. To identify baseline epigenetic markers associated 534 with response/non-response to treatment, we implemented a rigorous supervised machine 535 learning approach using stability selected gradient boosting34,35,104 combined with covered 536 information disentanglement (CID)105. Gradient boosting is an algorithm for supervised 537 learning, which operates through stepwise improvement of weak learners thereby minimizing 538 the overall prediction error against the observed data. We opted for gradient boosting as it 539 captures linear, non-linear and interaction effects better than traditional linear regression 540 approaches104. CID, in conjunction with stability selected gradient boosting, represents an 541 approach for feature selection that assigns permutation-based feature importance that, unlike 542 other methods for assigning feature importance, is unbiased by multicollinearity105. 543 We first removed CpGs with low variance and then applied univariate feature selection using 544 a F-test. Stability selection was subsequently employed to identify reliable biomarkers by 545 splitting randomly the discovery data into an 80% training data and a 20% test data using 546 stratified shuffle split and repeating this process 100 times to mitigate overfitting106. During 547 each split, we computed the CID for each CpG by randomly permuting it 100 times and 548 calculating the mean feature importance and assessing the average effect of permutation on 549 the model's performance (i.e., predicted versus true outcome)105. After all 100 iterations, the 550 mean feature importance per iteration was averaged and compared against a randomly 551 generated noise variable that was included throughout the entire modeling process where 552 CpGs with an aggregated feature importance ranked above a random variable were termed 553 predictor CpGs and retained for future analyses. 554 Having determined the predictor CpGs, we subsequently validated the predictive 555 performance internally and externally. Internal validation was performed by extracting the 556 predictor CpGs, training an ensemble of 100 gradient boost models on the 80% discovery 557 training data and utilizing all models to predict against the withheld 20% discovery test data. 558 External validation was performed using a similar approach where the discovery and 559 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint validation cohort were merged, the predictor CpGs extracted, and an ensemble of 100 560 models trained on the discovery cohort. The resultant models were then requested to predict 561 the response in the external validation cohort. In both cases, the output of each model 562 yielded a prediction score on a scale of 0 to 1 per sample. The prediction scores were 563 aggregated by calculating the mean prediction score per sample, representing the final, 564 ensembled, output of the model. This final prediction score was used to calculate the 565 receiver operator characteristic (ROC) to assess performance. The resultant prediction 566 scores were subsequently converted into classes by freezing the model at the Youden index, 567 thereby balancing the true positive rate relative to the false positive rate. Analyses were not 568 conducted separately for males and females as the cohorts would become too small. 569 Aforementioned analyses were performed in Python (v3.10) (https://www.python.org), with 570 packages scikit-learn and xgboost for the model development. 571 In silico DNA methylation analysis 572 Differential methylation analyses were performed in R using limma107 (v3.46) and eBayes108. 573 Separate analyses were run either regressing against the slide and slide position only, or in 574 combination with common confounding variables: age, sex, smoking behavior and the 575 estimated blood cell distribution109,110. The blood cell distribution was estimated using the 576

described by Houseman et al. 109 and implemented by Salas et al.110 for the Illumina 577 HumanMethylation EPIC BeadChip array, where methylation profiles of the current data are 578 compared against a reference methylation profiles of cell-sorted neutrophils, B cells, 579 monocytes, NK cells, CD4+ T cells and CD8+ T cells, enabling the inference of the cellular 580 composition. Intra-class correlation analyses were performed by conducting a two-way, 581 mixed, consistency analysis comparing samples obtained during response assessment with 582 samples obtained pretreatment using irr (v0.84.1). Statistical significance was defined as a 583 false discovery rate-adjusted p-value < 0.05. Visualizations were generated using ggplot2111 584 (v3.3.5). 585 RNA expression and data processing 586 Transcriptomic analyses was conducted through RNA sequencing, wherein mRNA was 587 extracted utilizing the QIAsymphony system, converted into cDNA and sequenced in a 588 paired-end format on the Illumina NovaSeq6000 at the Amsterdam UMC Core Facility 589 Genomics, generating a dataset comprising 40 million 150 bp-reads. Following that, in silico 590 data processing was orchestrated in Snakemake where read quality was assessed using 591 FastQC (v0.11.8) and summarized using MultiQC (v1.0)92,112. Raw reads were aligned to the 592 human genome (GRCh38) using the STAR aligner (v2.7.0), with annotations provided by the 593 Ensembl v95 annotation. Post-alignment processing was done in SAMtools (v1.9), followed 594 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint by read counting using the featureCounts function from the Subread package (v1.6.3)113-115. 595 Differential expression analysis (DE) analysis, was carried out within the R statistical 596 environment (v4.3) using the Bioconductor package DESeq2 (v1.38.3). We specifically 597 focused on genes associated with the predictor CpG loci based on the latter’s location in 598 either promoter or enhancer regions. Differentially expressed genes (DEGs) were identified 599 based on their significant differences, defined as those with a Benjamini–Hochberg-adjusted 600 p-value <0.05. Visualization of the results was accomplished using ggplot2 (v3.4.0)116. 601 Sample size estimation 602 We based our sample size on an initial pilot experiment and the calculations performed by 603 Tsai and Bell80 where we maintained a nominal p-value threshold of 0.05. A pilot experiment 604 using a subset of VDZ treated patients (7 R and 5 NR) indicated on average that the most 605 differentially methylated CpGs presented a mean difference in percentage methylation of 606 10% when comparing R with NR. We subsequently consulted the power calculations 607 reported by Tsai and Bell80, who conducted a case-control epigenome wide simulation study. 608 Assuming an approximately equal number of cases and controls and a mean difference in 609 percentage methylation of at least 10% at a nominal p-value threshold of 0.05, a statistical 610 power of at least 80% would be achieved if we included 40 patients (20 R and 20 NR) per 611 drug. To further eliminate the possibility of being underpowered, we aimed to collect at least 612 60 patients for the discovery cohort per drug, which would be supplemented by at least 20 613 patients for the validation cohort. 614 Statistical analysis of clinical variables 615 Baseline characteristics of all included patients were summarized using descriptive statistics. 616 Categorical variables are presented as percentages and continuous variables as median 617 annotated with the interquartile range (IQR). Differences in distribution between responders, 618 non-responders and the different cohorts were assessed using a chi-square test (categorical 619 variables) or Mann-Whitney U (continuous variables). Two-tailed probabilities were used with 620 a p-value ≤ 0.05 were considered statistically significant. Analyses of clinical data were 621 performed in IBM SPSS statistics (v26). To estimate the probability of a patient responding to 622 either VDZ or USTE after being predicted to be a responder, we calculated the post-test 623 probability using the following formulas: 624 (1) /g1868/g1870/g1857 /g1872/g1857/g1871/g1872 /g1867/g1856/g1856/g1871 /g3404 /g3043/g3045/g3032 /g3047/g3032/g3046/g3047 /g3043/g3045/g3042/g3029/g3028/g3029/g3036/g3039/g3036/g3047/g3052 /g2869/g2879/g3043/g3045/g3032 /g3047/g3032/g3046/g3047 /g3043/g3045/g3042/g3029/g3028/g3029/g3036/g3039/g3036/g3047/g3052 625 (2) /g1864/g1861/g1863/g1857/g1864/g1861/g1860/g1867/g1867/g1856 /g1870/g1853/g1872/g1861/g1867 /g3404 /g3046/g3032/g3041/g3046/g3036/g3047/g3036/g3049/g3036/g3047/g3052 /g2869/g2879/g3046/g3043/g3032/g3030/g3036/g3033/g3036/g3030/g3036/g3047/g3052 626 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint (3) /g1868/g1867/g1871/g1872 /g1872/g1857/g1871/g1872 /g1867/g1856/g1856/g1871 /g3404 /g1864/g1861/g1863/g1857/g1864/g1861/g1860/g1867/g1867/g1856 /g1870/g1853/g1872/g1861/g1867 /g1499 /g1868/g1870/g1857 /g1872/g1857/g1871/g1872 /g1867/g1856/g1856/g1871 627 (4) /g1868/g1867/g1871/g1872 /g1872/g1857/g1871/g1872 /g1868/g1870/g1867/g1854/g1853/g1854/g1861/g1864/g1861/g1872/g1877 /g3404 /g3043/g3042/g3046/g3047 /g3047/g3032/g3046/g3047 /g3042/g3031/g3031/g3046 /g3043/g3042/g3046/g3047 /g3047/g3032/g3046/g3047 /g3042/g3031/g3031/g3046/g2878/g2869 628 The pre-test probabilities were obtained from the largest VDZ7, and USTE8 treatment-629 response studies to date. The sensitivity and specificity were calculated by determining the 630 number of true positives, true negatives, false positives and false negatives when predicting 631 response in the validation cohort. 632 Data availability 633 The raw DNA methylation- (.idat) and gene expression (.fastq.gz) data alongside the de-634 identified patient metadata as reported on in this study have been published under controlled 635 access for research purposes at the European Genome-phenome Archive (EGA). The DNA 636 methylation data for the VDZ discovery cohort can be found under accession ID: 637 EGAD00010002651. The DNA methylation data for the VDZ validation cohort can be found 638 under accession ID: EGAD00010002652. The DNA methylation data for the USTE discovery 639 cohort can be found under accession ID: EGAD00010002649. The DNA methylation data for 640 the USTE validation cohort can be found under accession ID: EGAD00010002650. The 641 RNA-sequencing data for the VDZ discovery cohort can be found under accession ID: 642 EGAD50000000385. The RNA-sequencing data for the USTE discovery cohort can be found 643 under accession ID: EGAD50000000386. 644 Code availability 645 All Snakemake, bash calls and R scripts have been made available on GitHub and can be 646 found at https://github.com/ND91/HGPRJ0000008_EPICCD_multi_drug.git. The machine 647 learning modeling was performed using proprietary algorithms founded on the same 648 statistical principles as those of gradient boosting, permutation importance, and covered 649 information disentanglement. The code for these techniques is openly available at 650 https://xgboost.ai and https://github.com/JBPereira/CID or https://scikit-learn.org/stable. 651 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Tables 2 Table 1: Baseline characteristics discovery- and validation cohorts. Values in bold are significant, percentages shown are valid percentages. ADA: adalimumab. VDZ: 3 vedolizumab. USTE: ustekinumab. R: responder. NR: non-responder, SD: standard deviation. IQR: interquartile range, CRP: C-reactive protein. FCP: Fecal calprotectin. HBI: 4 Harvey Bradshaw Index. SES-CD: simple endoscopic disease activity score, Immunomodulator: azathioprine, mercaptopurine, thioguanine, methotrexate. Anti-TNF: 5 infliximab, adalimumab or golimumab. 6 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Amsterdam (discovery) Oxford (validation) Amsterdam (discovery) Oxford (validation) VDZ R (n=36) VDZ NR (n=28) P- val VDZ R (n=14) VDZ NR (n=11) P- val USTE R (n=30) USTE NR (n=32) P- val USTE R (n=22) USTE NR (n=11) P-val Sex, n (%), female 17 (47.2) 18 (64.3) 0.17 6 (42.9) 3 (27.3) 0.42 24 (80) 18 (56.3) 0.05 13 (59.1) 5 (45.5) 0.63 Age, years, median (IQR) 36 (25-52) 28 (32-56) 0.98 41 (21-65) 42 (24-70) 0.56 38 (29-56) 35 (23-45) 0.80 43 (35-55) 30 (24- 52) 0.25 Disease duration, years, median (IQR) 12 (3-19) 9 (4-20) 0.80 12 (2-15) 8 (3-17) 0.41 12 (7-23) 10 (5-21) 0.81 15 (6-24) 5 (2-12) 0.05 Ethnic background, n (%), White European 29 (80.6) 19 (67.9) 0.24 13 (92.9) 9 (81.8) 0.80 25 (83.3) 22 (68.8) 0.18 21 (95.5) 7 (87.5) 0.10 CRP, mg/L, median (IQR)a 3.9 (2.0- 11.7) 6.2 (3.2- 19.2) 0.26 5.2 (1.4- 11.6) 5.2 (1.3- 23.4) 0.97 5.7 (1.8- 19.2) 7.8 (3.4- 19.1) 0.45 6.1 (1.6- 18.0) 4.8 (2.3- 10.9) 0.30 FCP, ug/g, median (IQR)b 1003 (412- 1802) 1322 (389- 2654) 0.79 - - - 809 (180- 1599) 729 (221- 2347) 0.69 - - - Total baseline HBI, mean (±SD) or median (IQR) c 7.4 (±4.2) 9.3 (±4.9) 0.11 3 (1-8) 4 (1-8) 0.53 8.4 (±5.3) 8.2 (±5.9) 0.90 6 (4-9) 4 (2-8) 0.69 Total baseline SES-CD, median (IQR) 8 (6-12) 10 (5-13) 0.27 7 (6-12) 7 (6-22) 0.92 10 (6-17) 9 (6-14) 0.62 8 (5-16) 8 (3-19) 0.73 Endoscopic evaluation at follow-up 36 (100) 26 (93) 0.19 4 (29) 4 (36) 1.00 29 (97) 29 (91) 0.61 6 (27) 2 (18) 0.69 Disease location, n (%) - Ileal disease (L1) - Colonic disease (L2) - Ileocolonic disease (L3) - Upper GI involvement (L4) 13 (36.1) 9 (25.0) 14 (38.9) 1 (2.8) 8 (28.6) 5 (17.9) 15 (53.6) 2 (7.1) 0.50 0.41 3 (21) 2 (14) 9 (64) - 2 (18) 5 (45) 4 (36) - 0.21 - 6 (20.0) 7 (23.3) 17 (56.7) - 6 (18.8) 10 (31.3) 16 (50.0) - 0.78 - 7 (32) 8 (36) 6 (27) - 5 (46) 2 (18) 3 (27) 1 (9) 0.43 0.13 Disease behavior, n (%) 0.60 0.53 0.39 0.06 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint 7 - Non stricturing non-penetrating (B1) - Stricturing (B2) - Penetrating (B3) - Perianal disease (p) 16 (44.4) 10 (27.8) 10 (27.8) 8 (22.2) 9 (32.1) 10 (35.7) 9 (32.1) 11 (39.3) 0.14 9 (82) 2 (18) - 1 (7) 10 (91) 1 (9) - 2 (18) 0.48 10 (33.3) 11 (36.7) 9 (30.0) 10 (33.3) 12 (37.5) 15 (46.9) 5 (15.6) 12 (37.5) 0.73 15 (68) 5 (23) - 6 (27) 5 (46) 2 (18) 3 (27) 3 (27) 0.88 Previous IBD related surgery, n (%) 14 (38.9) 16 (57.1) 0.15 4 (28.6) 4 (36.4) 0.88 23 (76.7) 21 (65.6) 0.34 11 (50) 7 (63) 0.29 Concomitant medication, n (%) - Immunomodulators - Prednisone taper scheme 2 (5.6) - 2 (7.1) - 0.80 - - 2 - - - 1 (3.3) 7 (23.3) - 4 (12.5) 0.23 0.26 1 - - - - - Previous treatment exposure, n (%) - Immunomodulators - Anti-TNFs - Vedolizumab - Ustekinumab 28 (77.8) 24 (66.7) - 2 (5.6) 24 (85.7) 25 (89.3) - 7 (25.0) 0.42 0.03 - 0.03 13 (92.9) 5 (35.7) - 1 (7.1) 8 (72.7) 7 (63.9) - 6 (54.5) 0.17 0.25 - 0.01 29 (96.7) 30 (100) 12 (40.0) - 32 (100) 31 (96.9) 14 (43.8) - 0.23 0.25 0.77 - 17 (77.3) 14 (63.6) 2 (9.1) - 5 (45.5) 7 (63.3) 3 (27.3) - 0.07 1.00 0.15 - Smoking, n (%)d - Active 6 (17.1) 2 (7.1) 0.22 3 (21.4) 1 (9.1) 0.30 6 (20) 5 (15.6) 0.65 5 (22.7) 1 (9.1) 0.24 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Table 2. Predictive performance metrics of the prognostic biomarkers on the discovery cohort acquired at the AmsterdamUMC (n=126) and the validation cohort acquired at 8 the John Radcliffe Hospital (n=58). TP = True positives. TN = True negatives. FP = False positives. FN = False negatives. AUC = Area under the receiver operator 9 characteristic curve. VDZ: vedolizumab. USTE: ustekinumab. 6 0 VDZ discovery VDZ validation USTE discovery USTE validation TP 33 10 27 16 TN 22 8 29 8 FP 4 4 2 3 FN 3 3 2 6 AUC 0.87 0.75 0.89 0.75 Recall 0.92 0.77 0.93 0.73 Precision 0.89 0.71 0.93 0.84 F1-score 0.90 0.74 0.93 0.78 6 1 6 2 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Legends 663 Figures 664 Fig. 1: Predictive model using stability selected gradient boosting for response to therapy. a) An overview 665 of the feature selection and supervised machine learning approach for predicting response to VDZ and USTE as 666 used in the current study. Int. train: Internal training set. Int. test: Internal test set. ROC: Receiver operator 667 characteristic. b) Receiver operating characteristics plots showing the mean area under the curve (AUC) 668 performance of the discovery (n=126) and validation (n=58) cohorts. c) Left: radar plots presenting the 669 standardized difference in methylation between R (purple) and NR (green) for the top 15 predictor CpGs. Right: 670 Aggregated feature importance of the top 15 predictor CpGs. 671 Fig. 2: Longitudinal stability analyses. a) Volcano plot representing the differential methylation analyses when 672 comparing into treatment (T2) with pretreatment (T1) where grey dots represent CpG loci located on the Illumina 673 HumanMethylation EPIC BeadChip array and black dots represent response-associated predictor CpGs. X-axis 674 represents mean difference in percentage methylation, Y-axis represents the statistical significance as depicted 675 in –log10(p-value) with a dashed line at p-value = 0.05. b) Scatterplot showing the correlation of differential DNA 676 methylation between R and NR pretreatment (T1) and into treatment (T2). Grey dots represent CpG loci located 677 on the Illumina HumanMethylation EPIC BeadChip array and black dots represent response-associated predictor 678 CpGs. c) Boxplot of the two-way, consistency, intra-class correlation (ICC) coefficients of the predictor CpGs 679 calculated when comparing pretreatment and into treatment as well as the ICC coefficients of the predictor CpGs 680 obtained from a previous study on long-term stability of DNA methylation in IBD patients53. The vertical dashed 681 grey lines represent classification boundaries introduced by Koo and Li37, with blocks representing poor (ICC < 682 0.5), moderate (0.5 ≤ ICC < 0.75), good (0.75 ≤ ICC < 0.9), and excellent (0.9 ≥ ICC). d) Receiver operating 683 characteristic plots representing the predictive performance into treatment (T2; black) and pre-treatment (T1; 684 grey) as reference. 685 Fig. 3: Analyses of potential confounding variables. a) Volcano plot representing the change in response-686 associated differential methylation when correcting for the potential confounding variables age, sex, and 687 estimated cellular composition (black) or not (grey). X-axis represents mean difference in percentage 688 methylation, Y-axis represents the statistical significance as depicted in –log10(p-value). The dotted line 689 represents a threshold set at p = 0.05. b) Boxplot of the change in effect size after correcting for the confounding 690 variables as calculated by (β corrected- β uncorrected)/ β uncorrected. c) Receiver operator characteristic curve comparing 691 the response-prediction model for the CpG model (grey) with the confounder model (blue). 692 Fig. 4: Integrative analyses of predictor CpG-associated genes. a) Scatterplot showing the effect size (Wald 693 statistic) of the response-associated difference in pretreatment (T1; X-axis) and into treatment (T2; Y-axis). Grey 694 dots represent all genes measured, black dots represent non-differentially expressed predictor CpG-associated 695 genes, red, blue and purple dots represent predictor associated genes that are differentially expressed at T1, T2, 696 or T1 and T2, respectively. b) Boxplot of the T1 log2(expression) for VDZ response predictor CpG-associated 697 genes TULP4 and RFPL2 stratified by response. c) Scatterplot of the T1 log2(expression) of RFPL2 (X-axis) 698 relative to the percentage DNA methylation (Y-axis) for the predictor CpG cg12906381 annotated with the 699 Pearson correlation coefficient and the associated p-value. d) Boxplot of the T1 log2(expression) for USTE 700 response predictor CpG-associated genes MRC1 and TMEM191B stratified by response. e) Scatterplot of the T1 701 log2(expression) of TMEM191B (X-axis) relative to the percentage DNA methylation (Y-axis) for the predictor 702 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint CpG cg13982436 annotated with the Pearson correlation coefficient and the associated p-value. f) Boxplot of the 703 T2 log2(expression) for VDZ response predictor CpG-associated genes MCM2 and RFPL2 stratified by response. 704 g) Scatterplot of the T2 log2(expression) of RFPL2 (X-axis) relative to the percentage DNA methylation (Y-axis) 705 for the predictor CpG cg12906381 annotated with the Pearson correlation coefficient and the associated p-value. 706 h) Boxplot of the T2 log2(expression) for USTE response predictor CpG-associated genes POTEF, HDAC4, 707 PARP4 and MARK3 stratified by response. 708 Supplementary Tables 709 Supplementary Table 1: Differential methylation summary statistics of the VDZ-response predictor CpGs. 710 Columns represent the Illumina CpG identifier, the associated HGNC gene symbol, the chromosome, the location 711 on the chromosome (build: hg19), the mean difference in percentage methylation between responders and non-712 responders, the nominal p-value associated with the mean difference, the intra-class correlation coefficient 713 between pretreatment and during response assessment. 714 Supplementary Table 2: Differential methylation summary statistics of the USTE-response predictor CpGs. 715 Columns represent the Illumina CpG identifier, the mean difference in percentage methylation between 716 responders and non-responders, the nominal p-value associated with the mean difference, the genomic 717 coordinates on the human genome (build: hg19), the annotated gene name. 718 Supplementary Table 3: Clinical characteristics multi-biological failure cohort. 719 Supplementary Table 4: Differential expression summary statistics of the VDZ-response predictor CpGs-720 associated genes. Columns represent the Ensembl gene ID, the HGNC gene name, the mean difference in 721 log2(fold change) between responders and non-responders and the nominal p-value associated with the mean 722 difference. 723 Supplementary Table 5: Differential expression summary statistics of the USTE-response predictor CpGs-724 associated genes. Columns represent the Ensembl gene ID, the HGNC gene name, the mean difference in 725 log2(fold change) between responders and non-responders and the nominal p-value associated with the mean 726 difference. 727 Supplementary Figures 728 Supplementary Fig. 1: Longitudinal differential methylation of vedolizumab response-associated predictor CpGs. 729 Individual boxplots and temporal stability of the 25 VDZ-associated predictor CpGs over time. Top: jitterplot 730 visualization with samples obtained from the same donor connected by a dotted line annotated with the two-way 731 consistency intra-class correlation coefficient and the 95% confidence intervals. Bottom: boxplot visualization 732 grouped by response and stratified by timepoint. Annotations include nominal p-value as determined through 733 response-associated differential methylation analyses for samples obtained pretreatment (T1) or during response 734 assessment (T2). Green and red represent responders (R) and non-responders (NR), respectively. 735 Supplementary Fig. 2: Longitudinal differential methylation of ustekinumab response-associated predictor 736 CpGs. Individual boxplots and temporal stability of the 68 USTE-associated predictor CpGs over time. Top: 737 jitterplot visualization with samples obtained from the same donor connected by a dotted line annotated with the 738 two-way consistency intra-class correlation coefficient and the 95% confidence intervals. Bottom: boxplot 739 visualization grouped by response and stratified by timepoint. Annotations include nominal p-value as determined 740 through response-associated differential methylation analyses for samples obtained pretreatment (T1) or during 741 response assessment (T2). Green and red represent responders (R) and non-responders (NR), respectively. 742 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Supplementary Fig. 3: Stratified prediction analyses. Receiver operating characteristics plots showing the mean 743 area under the curve (AUC) performance of the validation cohort patients that were a) assessed using either the 744 strict (NVDZ = 8, NUSTE = 8; blue) or modified (NVDZ =17, NUSTE = 25; pink) response criteria, and b) either 745 anti-TNF exposed (NVDZ = 12, NUSTE = 21; black) or non-exposed (NVDZ = 13, NUSTE = 12; yellow) for VDZ 746 (left) and USTE (right). 747 Supplementary Fig. 4: Analyses of markers of IBD-associated inflammation. Volcano plot representing the 748 differential methylation analyses regressing against a) C-reactive protein (CRP) and b) fecal calprotectin (FCP) 749 pretreatment (T1) for VDZ (left) and USTE (right). Grey dots represent all CpG loci located on the Illumina 750 HumanMethylation EPIC BeadChip array and black dots represent response-associated predictor CpGs. X-axis 751 represents mean difference in percentage methylation, Y-axis represents the statistical significance as depicted 752 in –log10(p-value). CpGs with a nominal p-value < 0.05 are annotated. 753 754 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint

755 1. de Souza, H.S. & Fiocchi, C. Immunopathogenesis of IBD: current state of the art. 756 Nat Rev Gastroenterol Hepatol 13, 13-27 (2016). 757 2. Liu, J.Z. , et al. Association analyses identify 38 susceptibility loci for inflammatory 758 bowel disease and highlight shared genetic risk across populations. Nat Genet 47, 759 979-986 (2015). 760 3. Alsoud, D., Vermeire, S. & Verstockt, B. Biomarker discovery for personalized 761 therapy selection in inflammatory bowel diseases: Challenges and promises. Curr 762 Res Pharmacol Drug Discov 3, 100089 (2022). 763 4. Juillerat, P., Grueber, M.M., Ruetsch, R., Santi, G., Vuillemoz, M. & Michetti, P. 764 Positioning biologics in the treatment of IBD: A practical guide - Which mechanism of 765 action for whom? Curr Res Pharmacol Drug Discov 3, 100104 (2022). 766 5. Loftus, E.V., Jr., et al. Upadacitinib Induction and Maintenance Therapy for Crohn's 767 Disease. N Engl J Med 388, 1966-1980 (2023). 768 6. D'Haens, G., et al. Risankizumab as induction therapy for Crohn's disease: results 769 from the phase 3 ADVANCE and MOTIVATE induction trials. Lancet 399, 2015-2030 770 (2022). 771 7. Lowenberg, M. , et al. Vedolizumab Induces Endoscopic and Histologic Remission in 772 Patients With Crohn's Disease. Gastroenterology 157, 997-1006 e1006 (2019). 773 8. Sands, B.E., et al. Ustekinumab versus adalimumab for induction and maintenance 774 therapy in biologic-naive patients with moderately to severely active Crohn's disease: 775 a multicentre, randomised, double-blind, parallel-group, phase 3b trial. Lancet 399, 776 2200-2211 (2022). 777 9. Colombel, J.F., et al. Infliximab, azathioprine, or combination therapy for Crohn's 778 disease. N Engl J Med 362, 1383-1395 (2010). 779 10. Dulai, P.S. , et al. Development and Validation of Clinical Scoring Tool to Predict 780 Outcomes of Treatment With Vedolizumab in Patients With Ulcerative Colitis. Clin 781 Gastroenterol Hepatol 18, 2952-2961 e2958 (2020). 782 11. Verstockt, B., et al. Expression Levels of 4 Genes in Colon Tissue Might Be Used to 783 Predict Which Patients Will Enter Endoscopic Remission After Vedolizumab Therapy 784 for Inflammatory Bowel Diseases. Clin Gastroenterol Hepatol 18, 1142-1151 e1110 785 (2020). 786 12. Verstockt, B., et al. Low TREM1 expression in whole blood predicts anti-TNF 787 response in inflammatory bowel disease. EBioMedicine 40, 733-742 (2019). 788 13. Verstockt, B., et al. DOP70 An integrated multi-omics biomarker predicting 789 endoscopic response in ustekinumab treated patients with Crohn's disease. Journal 790 of Crohn's and Colitis 13, S072-S073 (2019). 791 14. Soendergaard, C., Seidelin, J.B., Steenholdt, C. & Nielsen, O.H. Putative biomarkers 792 of vedolizumab resistance and underlying inflammatory pathways involved in IBD. 793 BMJ Open Gastroenterol 5, e000208 (2018). 794 15. Ananthakrishnan, A.N. , et al. Gut Microbiome Function Predicts Response to Anti-795 integrin Biologic Therapy in Inflammatory Bowel Diseases. Cell Host Microbe 21, 796 603-610 e603 (2017). 797 16. Lee, J.W.J. , et al. Multi-omics reveal microbial determinants impacting responses to 798 biologic therapies in inflammatory bowel disease. Cell Host Microbe 29, 1294-1304 799 e1294 (2021). 800 17. Stevens, T.W. , et al. Systematic review: predictive biomarkers of therapeutic 801 response in inflammatory bowel disease-personalised medicine in its infancy. Aliment 802 Pharmacol Ther 48, 1213-1231 (2018). 803 18. Gisbert, J.P. & Chaparro, M. Predictors of Primary Response to Biologic Treatment 804 [Anti-TNF, Vedolizumab, and Ustekinumab] in Patients With Inflammatory Bowel 805 Disease: From Basic Science to Clinical Practice. J Crohns Colitis 14, 694-709 806 (2020). 807 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint 19. Zilbauer, M. & Heuschkel, R. Disease Prognostic Biomarkers in Inflammatory Bowel 808 Diseases-A Reality Check. J Crohns Colitis 16, 162-165 (2022). 809 20. Bird, A. DNA methylation patterns and epigenetic memory. Genes Dev 16, 6-21 810 (2002). 811 21. Kaluscha, S., et al. Evidence that direct inhibition of transcription factor binding is the 812 prevailing mode of gene and repeat repression by DNA methylation. Nat Genet 54, 813 1895-1906 (2022). 814 22. de Mendoza, A. , et al. Large-scale manipulation of promoter DNA methylation 815 reveals context-specific transcriptional responses and stability. Genome Biol 23, 163 816 (2022). 817 23. Hornschuh, M., Wirthgen, E., Wolfien, M., Singh, K.P., Wolkenhauer, O. & Dabritz, J. 818 The role of epigenetic modifications for the pathogenesis of Crohn's disease. Clin 819 Epigenetics 13, 108 (2021). 820 24. Peery, R.C., Pammi, M., Claud, E. & Shen, L. Epigenome - A mediator for host-821 microbiome crosstalk. Semin Perinatol 45, 151455 (2021). 822 25. Qin, Y. & Wade, P.A. Crosstalk between the microbiome and epigenome: messages 823 from bugs. J Biochem 163, 105-112 (2018). 824 26. Li, Y., et al. Intestinal mucosa-derived DNA methylation signatures in the penetrating 825 intestinal mucosal lesions of Crohn's disease. Sci Rep 11, 9771 (2021). 826 27. Sadler, T. , et al. Genome-wide analysis of DNA methylation and gene expression 827 defines molecular characteristics of Crohn's disease-associated fibrosis. Clin 828 Epigenetics 8, 30 (2016). 829 28. Adams, A.T. , et al. Two-stage Genome-wide Methylation Profiling in Childhood-onset 830 Crohn/i1s Disease Implicates Epigenetic Alterations at the VMP1/MIR21 and HLA 831 Loci. Inflammatory Bowel Diseases 20, 1784--1793 (2014). 832 29. Ventham, N.T., et al. Integrative epigenome-wide analysis demonstrates that DNA 833 methylation may mediate genetic risk in inflammatory bowel disease. Nature 834 Communications 7, 13507 (2016). 835 30. Li Yim, A.Y.F. , et al. Peripheral blood methylation profiling of female Crohn's disease 836 patients. Clin Epigenetics 8, 65 (2016). 837 31. Howell, K.J. , et al. DNA Methylation and Transcription Patterns in Intestinal Epithelial 838 Cells From Pediatric Patients With Inflammatory Bowel Diseases Differentiate 839 Disease Subtypes and Associate With Outcome. Gastroenterology 154, 585--598 840 (2018). 841 32. Gasparetto, M., et al. Transcription and DNA Methylation Patterns of Blood-Derived 842 CD8(+) T Cells Are Associated With Age and Inflammatory Bowel Disease But Do 843 Not Predict Prognosis. Gastroenterology (2020). 844 33. Li Yim, A.Y.F. , et al. Whole-Genome DNA Methylation Profiling of CD14+ Monocytes 845 Reveals Disease Status and Activity Differences in Crohn's Disease Patients. J Clin 846 Med 9(2020). 847 34. Meinshausen, N. & Bühlmann, P. Stability Selection Journal of the Royal Statistical 848 Society Series B: Statistical Methodology 72, 417/473 (2010). 849 35. Chen, T. & Guestrin, C. XGBoost: A Scalable Tree Boosting System. . in 850 Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge 851 Discovery and Data Mining 785–794 (New York, 2016). 852 36. Joustra V., L.Y.A. Long-term temporal stability of peripheral blood DNA methylation 853 alterations in patients with inflammatory bowel disease. (BioRxiv, 2022). 854 37. Koo, T.K. & Li, M.Y. A Guideline of Selecting and Reporting Intraclass Correlation 855 Coefficients for Reliability Research. J Chiropr Med 15, 155-163 (2016). 856 38. Tsai, P.C., et al. Smoking induces coordinated DNA methylation and gene 857 expression changes in adipose tissue with consequences for metabolic health. Clin 858 Epigenetics 10, 126 (2018). 859 39. Houseman, E.A., Kelsey, K.T., Wiencke, J.K. & Marsit, C.J. Cell-composition effects 860 in the analysis of DNA methylation array data: a mathematical perspective. BMC 861 Bioinformatics 16, 95 (2015). 862 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint 40. Solomon, O. , et al. Meta-analysis of epigenome-wide association studies in 863 newborns and children show widespread sex differences in blood DNA methylation. 864 Mutat Res-Rev Mutat 789(2022). 865 41. Dobbs, K.R., et al. Age-related differences in monocyte DNA methylation and 866 immune function in healthy Kenyan adults and children. Immun Ageing 18, 11 (2021). 867 42. Somineni, H.K. , et al. Blood-Derived DNA Methylation Signatures of Crohn Disease 868 and Severity of Intestinal Inflammation. Gastroenterology (2019). 869 43. Kappelman, M.D., et al. Real-World Evidence Comparing Vedolizumab and 870 Ustekinumab in Antitumor Necrosis Factor-Experienced Patients With Crohn's 871 Disease. Am J Gastroenterol 118, 674-684 (2023). 872 44. Torres, J. , et al. ECCO Guidelines on Therapeutics in Crohn's Disease: Medical 873 Treatment. J Crohns Colitis 14, 4-22 (2020). 874 45. Siegel, C., et al. P714 Mucosal healing in CD with vedolizumab versus other 875 biologics: endoscopic outcomes during long-term routine care in a multinational 876 observational study. Journal of Crohn's and Colitis 17, i844-i844 (2023). 877 46. Danese, S., et al. Treat to target versus standard of care for patients with Crohn's 878 disease treated with ustekinumab (STARDUST): an open-label, multicentre, 879 randomised phase 3b trial. Lancet Gastroenterol Hepatol 7, 294-306 (2022). 880 47. Feagan, B.G. , et al. Rapid Response to Vedolizumab Therapy in Biologic-Naive 881 Patients With Inflammatory Bowel Diseases. Clin Gastroenterol Hepatol 17, 130-138 882 e137 (2019). 883 48. Sandborn, W.J. , et al. Long-term efficacy and safety of ustekinumab for Crohn's 884 disease through the second year of therapy. Aliment Pharmacol Ther 48, 65-77 885 (2018). 886 49. Feagan, B.G. , et al. Ustekinumab as Induction and Maintenance Therapy for Crohn's 887 Disease. N Engl J Med 375, 1946-1960 (2016). 888 50. Mishra, N., et al. Longitudinal multi-omics analysis identifies early blood-based 889 predictors of anti-TNF therapy response in inflammatory bowel disease. Genome 890 Med 14, 110 (2022). 891 51. Lin, S., et al. Whole blood DNA methylation changes are associated with anti-TNF 892 drug concentration in patients with Crohn's disease. J Crohns Colitis (2023). 893 52. Kennedy, N.A., et al. Predictors of anti-TNF treatment failure in anti-TNF-naive 894 patients with active luminal Crohn's disease: a prospective, multicentre, cohort study. 895 Lancet Gastroenterol Hepatol 4, 341-353 (2019). 896 53. Joustra, V., et al. Long-term Temporal Stability of Peripheral Blood DNA Methylation 897 Profiles in Patients With Inflammatory Bowel Disease. Cell Mol Gastroenterol Hepatol 898 15, 869-885 (2023). 899 54. Ido, H., et al. Molecular dissection of the alpha-dystroglycan- and integrin-binding 900 sites within the globular domain of human laminin-10. J Biol Chem 279, 10946-10954 901 (2004). 902 55. Ulazzi, L., et al. Nidogen 1 and 2 gene promoters are aberrantly methylated in human 903 gastrointestinal cancer. Mol Cancer 6, 17 (2007). 904 56. Yuan, L. , et al. RhoJ is an endothelial cell-restricted Rho GTPase that mediates 905 vascular morphogenesis and is regulated by the transcription factor ERG. Blood 118, 906 1145-1153 (2011). 907 57. Park, S.Y. , et al. RhoA/ROCK-dependent pathway is required for TLR2-mediated IL-908 23 production in human synovial macrophages: Suppression by cilostazol. Biochem 909 Pharmacol 86, 1320-1327 (2013). 910 58. Zeng, R.J., Zhuo, Z.W., Luo, Y.J., Sha, W.H. & Chen, H. Rho GTPase signaling in 911 rheumatic diseases. Iscience 25(2022). 912 59. Sandi, M.J. , et al. MARK3-mediated phosphorylation of ARHGEF2 couples 913 microtubules to the actin cytoskeleton to establish cell polarity. Sci Signal 10(2017). 914 60. Hu, Y.J. , et al. PCGF3 promotes the proliferation and migration of non-small cell lung 915 cancer cells via the PI3K/AKT signaling pathway. Exp Cell Res 400(2021). 916 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint 61. Zhuang, G.Q. , et al. A Novel Regulator of Macrophage Activation miR-223 in 917 Obesity-Associated Adipose Tissue Inflammation. Circulation 125, 2892-+ (2012). 918 62. Jiao, P., et al. miR-223: An Effective Regulator of Immune Cell Differentiation and 919 Inflammation. Int J Biol Sci 17, 2308-2322 (2021). 920 63. Wright, P.B. , et al. The mannose receptor (CD206) identifies a population of colonic 921 macrophages in health and inflammatory bowel disease. Sci Rep 11, 19616 (2021). 922 64. Koelink, P.J. , et al. Anti-TNF therapy in IBD exerts its therapeutic effect through 923 macrophage IL-10 signalling. Gut 69, 1053-1063 (2020). 924 65. Lu, W. , et al. The microRNA miR-22 inhibits the histone deacetylase HDAC4 to 925 promote T(H)17 cell-dependent emphysema. Nat Immunol 16, 1185-1194 (2015). 926 66. Glauben, R., Sonnenberg, E., Wetzel, M., Mascagni, P. & Siegmund, B. Histone 927 Deacetylase Inhibitors Modulate Interleukin 6-dependent CD4(+) T Cell Polarization 928 in Vitro and in Vivo. J Biol Chem 289, 6142-6151 (2014). 929 67. Brusselle, G.G. & Bracke, K.R. MicroRNA miR-22 drives T(H)17 responses in 930 emphysema. Nat Immunol 16, 1109-1110 (2015). 931 68. Dou, B., Ma, F.Z., Jiang, Z.Y. & Zhao, L. Blood HDAC4 Variation Links With Disease 932 Activity and Response to Tumor Necrosis Factor Inhibitor and Regulates CD4+T Cell 933 Differentiation in Ankylosing Spondylitis. Front Med-Lausanne 9(2022). 934 69. Brusselle, G.G. & Bracke, K.R. MicroRNA miR-22 drives T(H)17 responses in 935 emphysema. Nat Immunol 16, 1109-1110 (2015). 936 70. Browning, L.M., et al. TGF-beta-mediated enhancement of TH17 cell generation is 937 inhibited by bone morphogenetic protein receptor 1alpha signaling. Sci Signal 938 11(2018). 939 71. Kuczma, M. & Kraj, P. Bone Morphogenetic Protein Signaling Regulates 940 Development and Activation of CD4(+) T Cells. Vitam Horm 99, 171-193 (2015). 941 72. Lu, L., et al. Synergistic effect of TGF-beta superfamily members on the induction of 942 Foxp3(+) Treg. Eur J Immunol 40, 142-152 (2010). 943 73. Park, S.H. , et al. Hypermethylation of EBF3 and IRX1 genes in synovial fibroblasts of 944 patients with rheumatoid arthritis. Mol Cells 35, 298-304 (2013). 945 74. Stockinger, B. & Veldhoen, M. Differentiation and function of Th17 T cells. Curr Opin 946 Immunol 19, 281-286 (2007). 947 75. Lyakh, L., Trinchieri, G., Provezza, L., Carra, G. & Gerosa, F. Regulation of 948 interleukin-12/interleukin-23 production and the T-helper 17 response in humans. 949 Immunol Rev 226, 112-131 (2008). 950 76. Zhao, L., Ghetie, D., Jiang, Z. & Chu, C.-Q. How Can We Manipulate the IL-23/IL-17 951 Axis? Current Treatment Options in Rheumatology 1, 182-196 (2015). 952 77. Morishima, N., Mizoguchi, I., Takeda, K., Mizuguchi, J. & Yoshimoto, T. TGF-beta is 953 necessary for induction of IL-23R and Th17 differentiation by IL-6 and IL-23. Biochem 954 Bioph Res Co 386, 105-110 (2009). 955 78. Rajkovic, A., Lee, J.H., Yan, C. & Matzuk, M.M. The ret finger protein-like 4 gene, 956 Rfpl4, encodes a putative E3 ubiquitin-protein ligase expressed in adult germ cells. 957 Mech Dev 112, 173-177 (2002). 958 79. Huang, C. Roles of E3 ubiquitin ligases in cell adhesion and migration. Cell Adh Migr 959 4, 10-18 (2010). 960 80. Tsai, P.C. & Bell, J.T. Power and sample size estimation for epigenome-wide 961 association scans to detect differential DNA methylation. Int J Epidemiol 44, 1429-962 1441 (2015). 963 81. Rees, C. & Penman, I. Has the COVID-19 pandemic changed endoscopy in the UK 964 forever? Lancet Gastroenterol Hepatol 8, 6-8 (2023). 965 82. Lappalainen, T. & Greally, J.M. Associating cellular epigenetic models with human 966 phenotypes. Nat Rev Genet 18, 441-451 (2017). 967 83. Villicana, S., et al. Genetic impacts on DNA methylation help elucidate regulatory 968 genomic processes. Genome Biol 24, 176 (2023). 969 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint 84. Yu, A.P., Cabanilla, L.A., Wu, E.Q., Mulani, P.M. & Chao, J. The costs of Crohn's 970 disease in the United States and other Western countries: a systematic review. Curr 971 Med Res Opin 24, 319-328 (2008). 972 85. Mehta, F. Report: economic implications of inflammatory bowel disease and its 973 management. Am J Manag Care 22, s51-60 (2016). 974 86. Jackson, B., Con, D., Ma, R., Gorelik, A., Liew, D. & De Cruz, P. Health care costs 975 associated with Australian tertiary inflammatory bowel disease care. Scand J 976 Gastroenterol 52, 851-856 (2017). 977 87. Ghosh, N. & Premchand, P. A UK cost of care model for inflammatory bowel disease. 978 Frontline Gastroenterol 6, 169-174 (2015). 979 88. Rubin, D.T., Mody, R., Davis, K.L. & Wang, C.C. Real-world assessment of therapy 980 changes, suboptimal treatment and associated costs in patients with ulcerative colitis 981 or Crohn's disease. Aliment Pharmacol Ther 39, 1143-1155 (2014). 982 89. Noor, N.M., et al. A biomarker-stratified comparison of top-down versus accelerated 983 step-up treatment strategies for patients with newly diagnosed Crohn's disease 984 (PROFILE): a multicentre, open-label randomised controlled trial. Lancet 985 Gastroenterol Hepatol 9, 415-427 (2024). 986 90. DNA methylation markers to predict tr eatment success of biologicals in Crohn’s 987 disease. 988 91. Vermeire, S., Schreiber, S., Sandborn, W.J., Dubois, C. & Rutgeerts, P. Correlation 989 between the Crohn's disease activity and Harvey-Bradshaw indices in assessing 990 Crohn's disease severity. Clin Gastroenterol Hepatol 8, 357-363 (2010). 991 92. Molder, F. , et al. Sustainable data analysis with Snakemake. F1000Res 10, 33 992 (2021). 993 93. Aryee, M.J., et al. Minfi: a flexible and comprehensive Bioconductor package for the 994 analysis of Infinium DNA methylation microarrays. Bioinformatics 30, 1363-1369 995 (2014). 996 94. Fortin, J.P., Triche, T.J., Jr. & Hansen, K.D. Preprocessing, normalization and 997 integration of the Illumina HumanMethylationEPIC array with minfi. Bioinformatics 33, 998 558-560 (2017). 999 95. Fortin, J.P., Fertig, E. & Hansen, K. shinyMethyl: interactive quality control of Illumina 1000 450k DNA methylation arrays in R. F1000Res 3, 175 (2014). 1001 96. Heiss, J.A. & Just, A.C. Identifying mislabeled and contaminated DNA methylation 1002 microarray data: an extended quality control toolset with examples from GEO. Clin 1003 Epigenetics 10, 73 (2018). 1004 97. Horvath, S. DNA methylation age of human tissues and cell types. Genome Biol 14, 1005 R115 (2013). 1006 98. Fortin, J.P., et al. Functional normalization of 450k methylation array data improves 1007 replication in large cancer studies. Genome Biol 15, 503 (2014). 1008 99. Javierre, B.M. , et al. Lineage-Specific Genome Architecture Links Enhancers and 1009 Non-coding Disease Variants to Target Gene Promoters. Cell 167, 1369-1384 e1319 1010 (2016). 1011 100. Johnson, W.E., Li, C. & Rabinovic, A. Adjusting batch effects in microarray 1012 expression data using empirical Bayes methods. Biostatistics 8, 118-127 (2007). 1013 101. Leek, J.T., Johnson, W.E., Parker, H.S., Jaffe, A.E. & Storey, J.D. The sva package 1014 for removing batch effects and other unwanted variation in high-throughput 1015 experiments. Bioinformatics 28, 882-883 (2012). 1016 102. Daca-Roszak, P. , et al. Impact of SNPs on methylation readouts by Illumina Infinium 1017 HumanMethylation450 BeadChip Array: implications for comparative population 1018 studies. BMC Genomics 16, 1003 (2015). 1019 103. Andrews, S.V., Ladd-Acosta, C., Feinberg, A.P., Hansen, K.D. & Fallin, M.D. "Gap 1020 hunting" to characterize clustered probe signals in Illumina methylation array data. 1021 Epigenetics Chromatin 9, 56 (2016). 1022 104. Friedman, J.H. Greedy Function Approximation: A Gradient Boosting Machine. The 1023 Annals of Statistics 29, 1189-1232 (2001). 1024 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint 105. J., P., E., S., A., Z. & E., L. Covered Information Disentanglement: Model 1025 Transparency via Unbiased Permutation Importance. in The Thirty-Sixth AAAI 1026 Conference on Artificial Intelligence (AAAI-22) (2022). 1027 106. Demsar, J. & Zupan, B. Hands-on training about overfitting. PLoS Comput Biol 17, 1028 e1008671 (2021). 1029 107. Ritchie, M.E., et al. limma powers differential expression analyses for RNA-1030 sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015). 1031 108. Smyth, G.K. Linear models and empirical bayes methods for assessing differential 1032 expression in microarray experiments. Stat Appl Genet Mol Biol 3, Article3 (2004). 1033 109. Houseman, E.A., et al. DNA methylation arrays as surrogate measures of cell 1034 mixture distribution. BMC Bioinformatics 13, 86 (2012). 1035 110. Salas, L.A., et al. An optimized library for reference-based deconvolution of whole-1036 blood biospecimens assayed using the Illumina HumanMethylationEPIC BeadArray. 1037 Genome Biol 19, 64 (2018). 1038 111. Wickham, H. ggplot2: Elegant Graphics for Data Analysis, (Springer-Verlag New 1039 York, 2016). 1040 112. Ewels, P., Magnusson, M., Lundin, S. & Kaller, M. MultiQC: summarize analysis 1041

for multiple tools and samples in a single report. Bioinformatics 32, 3047-3048 1042 (2016). 1043 113. Dobin, A., et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21 1044 (2013). 1045 114. Li, H., et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 1046 2078-2079 (2009). 1047 115. Liao, Y., Smyth, G.K. & Shi, W. featureCounts: an efficient general purpose program 1048 for assigning sequence reads to genomic features. Bioinformatics 30, 923-930 1049 (2014). 1050 116. H, W. ggplot2: Elegant Graphics for Data Analysis [Internet]. New York, NY: 1051 Springer-Verlag New York (2009). 1052 1053 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint Int. train (80%) Int. test (20%)#1 #2 #100 Feature importanceCpGs Predictor CpGs Random Prediction New data Training #1 #2 #100 Discovery (AmsterdamUMC) Validation (John Radcliffe Hospital) Train (80%) Test (20%) Feature selection Mean ROC ROCXGboostXGboost Youden index #1 #2 #100 a 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 Vedolizumab b 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 AUCDiscovery = 0.87 AUCValidation = 0.75 AUCDiscovery = 0.89 AUCValidation = 0.75 c Ustekinumab ROC curve Predictor CpGs NR R FPR TPR TPR FPR Feature importance Mean feature importance Mean feature importance . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint −5 0 5 10 0 2 4 6 −10 −5 0 5 10 15 0 2 4 6 -log10(p-value) Ustekinumab T2vT1Difference % methylation a −20 0 20 −20 0 20 Vedolizumab Spearman ρ = 0.971; p = 7.832E-07 T2: RvNRDifference % methylation −30 −20 −10 0 10 20 −30 −20 −10 0 10 20 Ustekinumab Spearman ρ = 0.941; p < 2.2E-16 T1: RvNRDifference % methylation b 0.00 0.25 0.50 0.75 1.00 Current study Joustra et al. 2022 Vedolizumab 0.00 0.25 0.50 0.75 1.00 Ustekinumabc ICC 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00TPR FPR Vedolizumab 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 AUCDiscovery T1 = 0.97 AUCDiscovery T2 = 1.00 Ustekinumabd AUCDiscovery T1 = 0.89 AUCDiscovery T2 = 0.97 Vedolizumab 8 15-10 . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint T1: RvNRDifference % methylation -log10(p-value) a Uncorrected Corrected −10 0 10 0 2 4 6 Vedolizumab −10 0 10 20 0 2 4 6 Ustekinumab Vedolizumab 0.0-0.5 (βcorrected-βuncorrected)/βuncorrected Ustekinumab 0 2-2 b 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 Vedolizumab Ustekinumab AUCCpG model = 0.75 AUCConfounder model = 0.57 AUCCpG model = 0.75 AUCConfounder model = 0.65 c . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint T1: RvNRWald statistic T1: RvNRWald statistic T2: RvNRWald statistic T2: RvNRWald statistic b c d e f g h TULP4 p = 0.0439 MCM2 p = 2.397E-03 RFPL2 p = 3.884E-03 POTEF p = 0.0147 HDAC4 p = 0.0219 PARP4 p = 0.0349 MARK3 p = 0.0287 MRC1 p = 4.255E-04 TMEM191B p = 7.318E-03 RFPL2 p = 0.0439 RFPL2 - cg12906381 Pearson r = -0.65; p = 1.394E-03 RFPL2 - cg12906381 Pearson r = -0.55; p = 0.017 TMEM191B - cg13982436 Pearson r = 0.57; p = 3.105E-03 NS T1 T2 T1&T2 log 2 (expression) log 2 (expression) % methylation log 2 (expression) log 2 (expression) % methylation log 2 (expression) log 2 (expression) log 2 (expression) % methylation 3.6 3.4 3.2 3.0 2.8 2.6 70 60 50 40 30 20 60 40 20 6 5 4 7.2 7.0 6.8 6.6 6.4 NR R NR R NR R NR R 2.6 2.8 3.0 3.2 3.4 3.6 1.5 2.0 2.5 NR R NR R NR R NR R NR R NR R 2.5 2.0 1.5 120 100 80 60 40 5.5 5.0 4.5 4.0 3.6 3.2 2.8 2.4 2.4 2.8 3.2 3.6 6.0 5.5 5.0 4.5 4.0 11.0 10.5 10.0 9.5 12.0 11.5 11.0 11.6 11.2 10.8 MCM2 TULP4 RFPL2 −6 −3 0 3 6 −6 −3 0 3 6 Vedolizumab POTEF HDAC4 MRC1 PARP4 MARK3 TMEM191B −6 −3 0 3 6 −6 −3 0 3 6 Ustekinumaba . CC-BY-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted July 25, 2024. ; https://doi.org/10.1101/2024.07.25.24310949doi: medRxiv preprint