**DNA Methylation of** *FKBP5* **as Predictor of Overall Survival in Malignant Pleural Mesothelioma**

**Giovanni Cugliari 1,\* , Chiara Catalano <sup>1</sup> , Simonetta Guarrera 2,3 , Alessandra Allione <sup>1</sup> , Elisabetta Casalone <sup>1</sup> , Alessia Russo <sup>1</sup> , Federica Grosso <sup>4</sup> , Daniela Ferrante 5,6, Clara Viberti <sup>1</sup> , Anna Aspesi <sup>7</sup> , Marika Sculco <sup>7</sup> , Chiara Pirazzini <sup>8</sup> , Roberta Libener <sup>9</sup> , Dario Mirabelli 10,11,**† **, Corrado Magnani 5,6,11,**† **, Irma Dianzani 7,11 and Giuseppe Matullo 1,11,12,\***


Received: 9 October 2020; Accepted: 18 November 2020; Published: 21 November 2020

**Simple Summary:** Our study is the first one to investigate DNA methylation changes in white blood cells (WBCs) from easily accessible peripheral blood as malignant pleural mesothelioma (MPM) survival biomarker. The Cox proportional hazards regression model highlighted that the methylation status of the CpG dinucleotide cg03546163 is an independent marker of prognosis in MPM patients with a better performance than traditional inflammation-based scores such as lymphocyte-to-monocyte ratio (LMR). Biological validation and replication showed that epigenetic changes at the *FKBP5* gene were robustly associated with overall survival (OS) in MPM cases. The identification of simple and valuable prognostic markers for MPM will enable clinicians to select patients who are most likely to benefit from aggressive therapies and avoid subjecting non-responder patients to ineffective treatment.

**Abstract:** Malignant pleural mesothelioma (MPM) is an aggressive tumor with median survival of 12 months and limited effective treatments. The scope of this study was to study the relationship between blood DNA methylation (DNAm) and overall survival (OS) aiming at a noninvasive prognostic test. We investigated a cohort of 159 incident asbestos exposed MPM cases enrolled in an Italian area with high incidence of mesothelioma. Considering 12 months as a cut-off for OS,

epigenome-wide association study (EWAS) revealed statistically significant (*p* value = 7.7 × 10−<sup>9</sup> ) OS-related differential methylation of a single-CpG (cg03546163), located in the 5′UTR region of the *FKBP5* gene. This is an independent marker of prognosis in MPM patients with a better performance than traditional inflammation-based scores such as lymphocyte-to-monocyte ratio (LMR). Cases with DNAm < 0.45 at the cg03546163 had significantly poor survival compared with those showing DNAm ≥ 0.45 (mean: 243 versus 534 days; *p* value< 0.001). Epigenetic changes at the *FKBP5* gene were robustly associated with OS in MPM cases. Our results showed that blood DNA methylation levels could be promising and dynamic prognostic biomarkers in MPM.

**Keywords:** malignant pleural mesothelioma; asbestos exposure; DNA methylation; lymphocyte-tomonocyte ratio; epigenome-wide analysis; survival analysis

#### **1. Introduction**

Malignant pleural mesothelioma (MPM) is an aggressive tumor. The disease usually develops after a long latency (20–40 years) following asbestos exposure [1]. Although MPM is considered a rare malignancy (prevalence 1–9/100,000), about 40,000 deaths have been estimated to occur each year globally [2,3]. The World Health Organization estimates that 125 million people annually around the world are exposed to asbestos. The International Agency for Research on Cancer confirmed that all fibrous forms of asbestos are carcinogenic to humans, causing mainly mesothelioma, respiratory-tract tumors, mesothelioma, and cancer at other tissue sites [4].

The prognosis of MPM is poor with a median survival of about 12 months from the diagnosis [5]. Generally, the first-line treatment is a combination of a multitargeted anti folate (pemetrexed or raltitrexed) drug and a platinum compound (cisplatin or carboplatin) [6] Currently, only a single randomized trial demonstrated an increase in survival time when comparing cisplatin and pemetrexed versus cisplatin alone [7]; unfortunately, most patients became resistant to this treatment and relapsed rapidly. No oncogenic driver has been identified and molecular pathways leading to MPM have not yet been clearly determined. Other therapeutic strategies such as immunotherapy are promising but require further investigation and improvement [8].

Recent research on the pathogenesis of MPM indicated that (i) both genetic and epigenetic alterations contribute to asbestos-induced tumorigenesis [9,10], (ii) inflammation-based prognostic scores that include lymphocyte counts are associated with survival [11].

MPM has a low frequency of protein-altering mutations (~25 mutations per tumor), compared to many other tumors [12]. Moreover, germline mutations in different genes mainly involved in DNA damage repair confer moderate-to-high genetic risk of MPM development [13]. The BAP1-tumor predisposition syndrome is the most studied genetic condition associated with MPM development and is caused by mutations in the BRCA1-associated protein 1 (*BAP1*) gene [13].

In the last 10 years, epigenetic markers, such as DNA methylation (DNAm) and microRNAs (miRNAs), have gained popularity as possible early diagnostic and prognostic biomarkers in cancer research, including MPM. While genetic markers may differ from case to case in most cancer patients (i.e., each patient may carry a different mutation within the same gene), different subjects show variable levels of epigenetic biomarkers in specific target regions and different tissues depending on disease status [14].

DNA methylation is one of the epigenetic factors [15] that can be altered in cancer tissues. However, regarding mechanisms and clinical outcome of epigenetic derangements in MPM, less information is available [16,17] Although DNAm is stable, it can be modified throughout life by several factors such as ageing, lifestyle, environmental exposures, and diseases. It thus represents an adaptive phenomenon linking environmental factors and the development of pathologic phenotypes such as

cancers. DNAm changes are considered to possibly play a role in MPM progression, and have therefore been suggested as a potential tool for prognosis [18].

The fact that epigenetic modifications, unlike genetic changes, are potentially reversible, may open new perspectives for patient clustering and novel therapeutic options. A reliable prognostic biomarker that offers high sensitivity and specificity would be a major advancement for MPM. Blood-based biomarkers that have been explored in MPM include megakaryocyte potentiating factor (an alternative cleavage product of the mesothelin precursor protein) [19], and Fibulin 3 which is also found in pleural fluid, and whose high levels appear to correlate with advanced disease [20].

Considering clinical end-point, low pleural fluid glucose and high C-reactive protein and pleural thickening represent the main prognosis factors [21]. Recent studies confirm that using also a combination of epigenetic alterations as biomarkers is more informative with respect to an only genetic approach on overall survival (OS) [17].

This study was undertaken with the goal of better characterizing the MPM OS evaluating the potential predictive value of peripheral blood DNAm profiles. The second goal was the comparison of the DNAm prognostic performance with the broadly used lymphocyte-to-monocyte ratio (LMR) method.

#### **2. Results**

#### *2.1. Epigenome-Wide Association Study (EWAS)*

EWAS revealed a statistically significant hypo-methylated single-CpG (cg03546163) in the *FKBP5* gene in the low survival group after Bonferroni post-hoc correction (Figure 1).

**Figure 1.** Manhattan plot for epigenome-wide association study (EWAS) test on 450k single CpGs. Overall survival was used as dependent variable considering 12 months as cut-off adjusting for age, gender, histological subtype, asbestos exposure, WBCs estimation, population stratification, and technical variability. Bonferroni post hoc line highlights statistically significant differences on OS at single CpG level.

Bootstrap was computed to estimate the measures of accuracy, using random sampling methods. The other five CpGs in the *FKBP5* gene showed hypomethylation in poor MPM survivors, with unadjusted *p* value < 0.05 (Table 1); instead, no CpGs in the *FKBP5* gene showed statistically significant hypermethylation in poor MPM survivors.

**Table 1.** Differential DNAm analyses of the *FKBP5* gene ordered by effect size (low survival group was used as reference). Information about single CpGs including location-related values and model outputs (effect size, standard error,*p*values).


Low survival group was set as reference. Adjustment covariates: age, gender, asbestos exposure, histological subtype, smoke, population stratification, WBCs estimation, and technical variability. \*: statistically significant at*p*value<0.05; §: statistically significant at Bonferroni and FDR post hoc adjustments.

#### *2.2. Survival Analysis* ≥

CpG sites and LMR were considered as predictors in the regression model. Categorical variables (quantile information) were used. ≥

Cox model was computed considering the same list of covariates included in the EWAS. Patients with DNAm < 0.45 at the cg03546163 had significantly poorer survival compared with subjects with DNAm ≥ 0.45 (mean, 243 versus 534 days; *p* value < 0.001). Survival at the 1st and the 3rd Quartiles was 135 versus 209 days and 401 versus 842 days, respectively, comparing patients with single CpG DNAm < 0.45 with those with single CpG DNAm ≥ 0.45. The multivariate analysis showed that cg03546163 DNAm at *FKBP5* was independently associated with OS. Kaplan–Meier curves revealed that a decrease of methylation at cg03546163 (<0.45) was significantly associated with worse OS (HR = 2.14 *p* value < 0.0001) (Figure 2a). ≥

≥ ≥ **Figure 2.** K-M survival curves show (**a**) cg03546163: patients with a DNAm < 0.45 had significantly poor survival compared with a DNAm ≥ 0.45 (mean, 243 versus 534 days; *p* value < 0.001); (**b**) LMR: patients with values < 2.86 had significantly poor survival compared with patients with values ≥ 2.86 (mean, 310 versus 528 days; *p* value < 0.001). cg03546163 is an independent marker of prognosis in patients with MPM and performs better than LMR (HRcg03546163 = 2.14 vs HRLMR = 1.66).

Patients with LMR<2.86 had significantly poorer survival compared with patients with LMR≥2.86 (mean, 310 versus 528 days; *p* value < 0.001). Survival at 1st Quartile was 175 versus 262 days whereas at 3rd Quartile was 484 versus 969 days comparing patients with LMR < 2.86 with those with LMR > 2.86. LMR was independently associated with OS: Kaplan–Meier curves showed that decreased LMR (<2.86) was significantly associated with decreased OS (HR = 1.66; *p* value < 0.01) (Figure 2b).

Histological subtype (epithelioid versus non-epithelioid), smoking status (current, never, and former), and asbestos exposure showed no statistically significant results on survival.

#### *2.3. Validation and Replication*

The statistically significant association between cg03546163 DNAm and OS was confirmed in an independent sample of patients (replication) and using a different targeted DNAm analysis technique (validation). A sample of 133 MPM cases (58 low survivors and 75 high survivors) was recruited and stratified in low and high OS considering the same cut-off (365 days).

The same model used for the discovery phase was performed. Patients with below median OS had significantly lower DNAm at the cg03546163 compared with those with above median OS (mean, 188 versus 786 days; *p* value < 0.001). The 1st Quartile was 113 versus 482 days and the 3rd Quartile was 262 versus 862 days comparing patients with DNAm difference (reference above median OS, MD: −0.04, 95%CI: −0.07|−0.01, *p* value: 0.04) at the cg03546163. The multivariate analysis confirmed that cg03546163 DNAm at *FKBP5* was independently associated with OS.

#### **3. Discussion**

A growing number of studies reported on the identification of epigenetic prognostic biomarkers in several cancers [17,22].

This study focused on the exploration of epigenetic factors related to MPM survival in MPM incident cases from Piedmont (Italy), a region with a well-documented history of asbestos exposure [23].

More than 450k methylation sites were evaluated in DNA from whole blood looking for new insights related to overall survival in MPM. The main result was the hypomethylation of a single CpG (cg03546163) in the 5′ UTR region of the FKBP5 gene in patients with poorer survival compared to patients with longer survival; it also showed to be an independent marker of prognosis in MPM patients. This result was replicated in a different series of patients belonging to the same cohort using the Sequenom Quantitative DNAm analysis.

In general, a combination of epigenetic and clinical factors is under investigation in clinical prognosis and survival, including tumor histology, gender, hemoglobin level, platelet and white blood cell count, and lactate dehydrogenase level [24].

Recently, due to the important role of inflammation in the development of MPM, several studies investigated the effect of inflammation-based biomarkers on the prognosis [11,22]. We selected the LMR for the comparison because its performance was previously reported to be higher than other inflammation-based markers in MPM [25].

To validate the prognostic value of the observed CpG methylation site, we compared our result with the LMR score.

Kaplan–Meier survival curves for MPM patients highlighted cg03546163 methylation at FKBP5 gene as a prognostic factor superior to the LMR score.

The FKBP Prolyl Isomerase 5 (FKBP5), also known as FK506 binding protein 51 (FKBP51), is a member of the immunophilin protein family, which contributes to the immunoregulation and to the basic cellular processes involving protein folding and trafficking. Together with other members of the FKBPs family, this protein participates in transcriptional complexes and acts as a co-transcription factor.

Although no studies have investigated the methylation of *FKBP5* as prognostic factor in MPM, a growing number of whole-blood studies investigated its DNA methylation levels in order to explain the impact of environmental stress in the etiology and treatment of several diseases [26]. Interestingly, in a recent study on the Behcet's disease (BD) hypomethylation in the 5′UTR region (including cg03546163) of FKBP5 characterized cases was demonstrated and it was strongly associated with high gene expression, suggesting a possible role of DNA methylation in the pathogenesis [27].

Other five single CpGs at FKBP5 showed hypomethylation in poor survivors: this evidence supports the potential overall contribution of *FKBP5* methylation on the patient classification by OS.

In several human cancer tissues, a relevant role for FKBP5 in sustaining cancer cell growth and aggressiveness has been documented. In particular, for glioma [28], prostate cancer and melanoma [29] a strict correlation between protein abundance and aggressiveness has been demonstrated.

Probably, the relationship between FKBP5 and tumor progression and aggressiveness, is represented by its implication in NF-kB and AKT signaling pathways, with key roles in tumorigenesis and response to antineoplastic chemotherapy [30].

Moreover, a well characterized antiapoptotic effect is mediated by NF-κB transcription factors and FKBP5 has documented antiapoptotic effects: recent studies hypothesized that FKBP5 could promote inflammation, by activating the master immune regulator NF-kB, after an epigenetic upregulation due to aging and stress [31,32].

Previous studies conducted on various cancer types, showed that upregulation of FKBP5 gene expression is associated with drug resistance [33]. In a study on an ovarian cancer cell line, the upregulation of FKBP5 increased the resistance to chemotherapeutic agents, whereas the gene silencing sensitized ovarian cancer cells to taxol [34]. In the present study we could not evaluate FKBP5 gene expression due to the lack of available RNA, which was not collected in the study. However, this should be further addressed and verified in future studies.

One study demonstrated that overexpression of FKBP5 increased the chemosensitivity through the AKT pathway [31]. A similar study supported this observation making FKBP5 an effective biomarker for sensitivity to chemotherapy; patient responses to chemotherapy may be determined by the variation in FKBP5 levels [35].

#### *Limitation of the Study*

Being able to identify the direction of causality will greatly aid in determining the usefulness of epigenetic variation.

Leukocyte DNA methylation could mainly represent a nonspecific marker related to a general inflammatory status due to the presence of a tumor rather than a specific MPM biomarker and further studies should be carried out to support our findings.

As additional limitation, we had therapy information only for a small subset of patients and we could not test treatment-specific OS differences in relation to FKBP5 methylation levels.

#### **4. Material and Methods**

#### *4.1. Study Population*

Study subjects belong to a wider ongoing collaborative study on MPM, which is actively enrolling MPM cases in the municipalities of Casale Monferrato (Piedmont region, Italy), an area with an exceptionally high incidence of mesothelioma caused by widespread asbestos exposure for locals, both occupational and environmental, due to the asbestos-cement Eternit plant that was operational until 1986 [36]. Additional MPM cases were recruited in the main hospitals of the municipalities of Turin, Novara, and Alessandria (Piedmont region, Italy). The study included incident MPM cases diagnosed between 2000 and 2010 after histological and/or cytological confirmation of MPM diagnosis [37,38].

No peritoneal cases were considered with the aim to better identify epigenetics characteristics of MPM.

In the present study, 159 MPM cases belonging to a larger case–control study with genetic [10,39] and blood DNAm data [9] were selected according to the following criteria: (i) availability of good quality DNA at the time of the analyses and (ii) asbestos exposure above the background level, as defined in [40]. An additional 133 independent samples from the same cohort were included for the validation/replication analyses.

Descriptive information of MPM patients are shown in Table 2. Median survival (365 days) was used as cut-off value to stratify patients in high and low survivors.

No differences in categorical (center, gender, smoke, histotype) and continuous (asbestos exposure, WBCs composition) variables among low and high survivors were found.

Our study complies with the Declaration of Helsinki principles and conforms to ethical requirements. All volunteers signed an informed consent form at enrollment. The study protocol was approved by the Ethics Committee of the Italian Institute for Genomic Medicine (prot.n.CE-2015-GM-2, 30/10/2015, HUGEF, Turin, Italy).

#### *4.2. Exposure Assessment*

For all subjects, occupational history and lifestyle habits information were collected through interviewer-administered questionnaires filled out at enrollment during a face-to-face interview. Job titles were coded according to the International Standard Classification of Occupations [40] and according to the Statistical Classification of Economic Activities in the European Community.

Frequency, duration, and intensity of exposure were estimated, then a cumulative exposure index was computed. The evaluation of asbestos exposure (occupational, environmental, and domestic) was conducted by an experienced occupational epidemiologist. For the selection criteria and descriptive evaluation, asbestos exposure doses (fibers/mL years) were rank transformed to remove skewness.


**Table 2.** Descriptive information of MPM patients. Median survival (365 days) was used as cut-off value to stratify patients in high and low survivors.

Asbestos exposure (occupational, environmental, and domestic) was normalized considering frequency, duration, and intensity.

#### *4.3. Blood DNAm Analysis*

Genomic DNA was extracted from whole blood collected in EDTA by an on-column DNA purification method (QIAamp DNA Blood Mini Kit, QIAGEN GmbH, Germany), according to manufacturer's instructions. DNA integrity was checked by an electrophoretic run in standard TBE 0.5× buffer on a 1% low melting agarose gel (Sigma-Aldrich GmbH, Schnelldorf, Germany); DNA purity and concentration were assessed by a NanoDrop 8000 Spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). Five hundred nanograms of genomic DNA for each sample were bisulfite treated (EZ-96 DNA Methylation-Gold Kit, Zymo Research Corporation, Irvine, CA, USA) to convert un-methylated cytosine to uracil. Cases were randomly and blindly distributed across conversion plates.

The Infinium HumanMethylation450 BeadChip (Illumina Inc., San Diego, CA, USA) was used to measure the methylation level of more than 485,000 individual CpG loci at a genome-wide resolution [41].

Twelve samples were analyzed on each BeadChip. As a "position effect" was reported for Illumina Methylation BeadChips, each sample position on the BeadChip was completely random as well. We further verified the randomization of the position on each BeadChip was effective by checking for a

position effect, and we found no occurrence of it. BeadChips were processed according to manufacturer protocols. Data were inspected with the dedicated GenomeStudio software v2011.1 with Methylation module 1.9.0 (Illumina Inc., San Diego, CA), and quality checked as previously described [42].

#### *4.4. Beta-Value Extraction*

Raw DNAm data were analyzed with the R package (methylumi'). The average methylation value at each locus was computed as the ratio of the intensity of the methylated signal over the total signal (un-methylated + methylated) [43]. Beta-values represent the percentage of methylation at each individual CpG locus, ranging from 0 (no methylation) to 1 (full methylation).

We excluded from the analyses (i) single Beta-values with detection *p* value ≥ 0.01; (ii) CpG loci with missing Beta-values in more than 20% of the assayed samples; (iii) CpG loci detected by probes containing SNPs with MAF ≥ 0.05 in the CEPH (Utah residents with ancestry from northern and western Europe, CEU) population; (iv) samples with a global call rate ≤ 95%. Lastly, CpGs on chromosomes X and Y were excluded from the analysis.

#### *4.5. Batch E*ff*ect, Population Stratification, and White Blood Cell Estimations*

To account for methylation assay variability and batch effects, we corrected all differential methylation analyses for "control probes" principal components (PCs). Using PCs assessed by principal component analysis of the BeadChip's built-in control probes as a correction factor for statistical analyses of microarray data is a method that allows to account for the technical variability of several steps in the DNAm analysis, from the bisulfite conversion to BeadChip processing [44].

Geographic origins of subjects may influence DNAm profiles. To consider this source of potential bias, we took advantage of the whole genome genotyping dataset from the same subjects from our previous study [10]. The first PCs calculated based on genome-wide genotyping were shown to correlate with different geographic origins of people [45,46].

WBC subtype percentages calculated based on genome-wide methylation data [47] for each subject were extracted. This method quantifies the normally mixed composition of leukocytes beyond what is possible by simple histological or flow cytometric assessments. In a diverse array of diseases and following numerous immune-toxic exposures, leukocyte composition will critically inform the underlying immune-biology to most chronic medical conditions. Then, it is necessary to extract and control for the percentage of involved WBCs with the aim to infer about a functional biological pathway.

LMR score was calculated from the DNAm-estimated WBCs by dividing the total lymphocyte count by the monocyte count.

#### *4.6. Statistical Analyses*

#### Epigenome-Wide Association Study

Association test was used to analyze the mean differences (MD) at single-CpG methylation between low and high survival. Multiple regression analysis adjusted for age, gender, histological subtype, asbestos exposure, smoke, estimated WBCs, population stratification (first 2 PCs) and technical variability (first 10 PCs) was implemented. For multiple comparisons tests, Bonferroni *p* value ≤ 0.05 was considered statistically significant.

Using random sampling methods, bootstrap was implemented to estimate the measures of accuracy defined in terms of bias, variance, confidence intervals, and prediction error. Bootstrap is also an appropriate way to control and check the stability of the results. The bias-corrected and accelerated (BCa) bootstrap interval was calculated with regard to single CpGs.

#### *4.7. Survival Analysis*

The survival time was determined as the time between the date of diagnosis and the date of death. If patients were still alive at the last follow-up (2016), survival was defined as the time from the date of

diagnosis until June 2016. The time and the median event times with 95% confidence intervals were estimated according to the Kaplan–Meier method. The proportional hazards regression model was used for both the univariate and multivariate analyses (Cox model).

Comparison of OS curves was performed using two-tailed log-rank tests with a 0.05 level of significance. Only variables with *p* value < 0.1 in the univariate analysis were included in the final model for the multivariate analysis. In the Cox regression analysis, the backward conditional method (stepwise-AIC) was used. LMR and CpG sites were considered as predictors in regression model.

#### *4.8. Statistical Power*

To ensure a power of the study greater than 80% (two-tailed test at 0.05 alpha error), only CpGs with mean difference (MD) of Beta-value between low and high survival of ≥ |0.035| were selected. Covariates were included step-by-step in sensitivity analysis to validate the association output considering effect size, standard error, 95% confidence interval and *p* value variations.

CpGs with Bonferroni *p* value ≤ 0.05 underwent gene set enrichment analysis to identify pathways potentially affected by MPM related methylation changes.

All statistical analyses were conducted using the open source software R (4.0.2).

#### *4.9. Validation and Replication*

Sequenom MassARRAY for the DNAm signal validation and replication was used. In detail, the EpiTYPER assay (Sequenom) uses a MALDI-TOF mass spectrometry-based method to quantitatively assess the DNA methylation state of CpG sites of interest [48]. DNA (500 ng) was bisulfite-converted using the EZ-96 DNA Methylation Kit (Zymo Research) with the following modifications: incubation in CT buffer for 21 cycles of 15 min at 55 ◦C and 30 s at 95 ◦C, elution of bisulfite-treated DNA in 100 µL of water. The treatment converts unmethylated Cytosine into Uracil, leaving methylated Cytosine unchanged. In this way, variations in the sequence are produced depending on DNA methylation status of the original DNA molecule.

PCR amplification, treatment with SAP solution, and Transcription/RNase A cocktails were performed according to the protocol provided by Sequenom and the mass spectra were analyzed by EpiTYPER analyzer (Sequenom, San Diego, CA, USA). As the MassARRAY assay is unable to discriminate between CpGs located at close vicinity to each other in the sequence, the close neighboring CpGs were analyzed as "Units", i.e., the measured methylation level is the average of the methylation levels of the CpGs cumulatively analyzed within the Unit. In the case of cg03546163 the measured methylation level is the average between two CpG sites located very close (Figure S1).

The amplicon for cg03546163 (chr6:35,654,364) encompasses 196bp (chr6:35,654,222-chr6:35,654,418 (GRCh37/hg19)) and PCR was performed on 10 ng of converted DNA using the following primers:


#### **5. Conclusions**

Our results suggest the potential use of DNAm analysis in blood to develop noninvasive tests for prognostic evaluation in MPM; our study is the first to demonstrate that a single CpG in *FKBP5* gene is an independent marker of prognosis in patients with MPM and is superior to the LMR inflammation-based prognostic score. The identification of simple and valuable prognostic markers for MPM will enable clinicians to select patients who are most likely to benefit from aggressive therapies and avoid subjecting nonresponder patients to ineffective treatment. Moreover, epigenetic modifications such as DNAm are potentially reversible and can open new perspectives for epigenetic therapies in MPM. Knowledge of epigenetic changes has provided new therapeutic opportunities against cancer. To allow better approach of cancer cell inhibitory strategies, the understanding of molecular mechanisms that underlie cellular DNA epigenetic alterations may be useful. In this context, we reported epigenetic deregulations in blood samples from MPM patients in relation to OS, paving the road to both patients' stratification and the possible discovery of new combined therapeutic options in MPM. Studies of a large population are needed to investigate the relationship between prognostic markers and treatment regimens. The usage of methylation alterations in clinical specimens as biomarkers could be recognized. Noninvasively obtained, methylation-based biomarkers detected in blood cells from cancer patients offer significant practical advantages, being promising and dynamic prognostic markers.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-6694/12/11/3470/s1, Figure S1. Locations of cg03546163 (CpG2, in bold) and a second CpG site very close (CpG1, in red), investigated by Sequenom MassARRAY.

**Author Contributions:** Conceptualization, G.C., G.M., C.M., D.M., I.D. and F.G.; methodology, G.C., G.M. and S.G.; software, G.C.; validation, G.C., G.M., A.A. (Alessandra Allione), A.R., E.C. and C.P.; formal analysis, G.C.; investigation, G.C., G.M., S.G., A.A. (Alessandra Allione), A.R., E.C. and C.C.; resources, G.C., G.M., C.M., D.M., I.D., F.G., C.V., D.F., A.A. (Anna Aspesi), M.S. and R.L.; data curation, G.C., C.M., D.M., I.D., F.G., C.P., C.V., D.F., A.A. (Anna Aspesi), M.S. and R.L.; writing—original draft preparation, G.C.; writing—review and editing, G.C., G.M., C.M., D.M., I.D., F.G., S.G., A.A. (Alessandra Allione), A.R., E.C., C.C., C.P., C.V., D.F., A.A. (Anna Aspesi), M.S. and R.L.; visualization, G.C., G.M., S.G., A.A. (Alessandra Allione), A.R., E.C., C.C., C.P., C.V., D.F., A.A. (Anna Aspesi), M.S. and R.L; supervision, G.C. and G.M.; project administration, G.C. and G.M.; funding acquisition, G.M., C.M., D.M., I.D. and F.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research leading to these results has received funding from AIRC under IG 2018—[ID. 21390 project—P.I. GM] and partly by the HERMES (Hereditary Risk in MESothelioma) Project, funded by the offer of compensation to the inhabitants of Casale Monferrato deceased or affected by mesothelioma (to I.D. and C.M.). and the Ministero dell'Istruzione, dell'Università e della Ricerca—MIUR project "Dipartimenti di Eccellenza 2018–2022" (n◦ D15D18000410001, to GM) to the Department of Medical Sciences, University of Turin.

**Acknowledgments:** The authors wish to thank all the patients contributing to this research, Caterina Casadio and Ezio Piccolini for their past clinical contribution and Paolo Garagnani for the Sequenom analysis supervision.

**Conflicts of Interest:** The authors declare that they have no competing interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **New DNA Methylation Signals for Malignant Pleural Mesothelioma Risk Assessment**

**Giovanni Cugliari 1,\* , Alessandra Allione <sup>1</sup> , Alessia Russo <sup>1</sup> , Chiara Catalano <sup>1</sup> , Elisabetta Casalone <sup>1</sup> , Simonetta Guarrera 2,3, Federica Grosso <sup>4</sup> , Daniela Ferrante 5,6, Marika Sculco <sup>7</sup> , Marta La Vecchia <sup>7</sup> , Chiara Pirazzini <sup>8</sup> , Roberta Libener <sup>9</sup> , Dario Mirabelli 10,11, Corrado Magnani 5,6,11 , Irma Dianzani 7,11 and Giuseppe Matullo 1,11,12,\***


**Simple Summary:** Our study investigated DNA methylation differences in easily accessible white blood cells (WBCs) between malignant pleural mesothelioma (MPM) cases and asbestos-exposed cancer-free controls. A multiple regression model highlighted that the methylation level of two single CpGs (cg03546163 in *FKBP5* and cg06633438 in *MLLT1*) are independent MPM markers. The epigenetic changes at the *FKBP5* and *MLLT1* genes were robustly associated with MPM in asbestosexposed subjects. Interaction analyses showed that MPM cases and cancer-free controls showed DNAm differences which may be linked to asbestos exposure.

**Abstract:** Malignant pleural mesothelioma (MPM) is a rare and aggressive neoplasm. Patients are usually diagnosed when current treatments have limited benefits, highlighting the need for noninvasive tests aimed at an MPM risk assessment tool that might improve life expectancy. Three hundred asbestos-exposed subjects (163 MPM cases and 137 cancer-free controls), from the same geographical region in Italy, were recruited. The evaluation of asbestos exposure was conducted considering the frequency, the duration and the intensity of occupational, environmental and domestic exposure. A genome-wide methylation array was performed to identify novel blood DNA methylation (DNAm) markers of MPM. Multiple regression analyses adjusting for potential confounding factors and interaction between asbestos exposure and DNAm on the MPM odds ratio were applied. Epigenome-wide analysis (EWAS) revealed 12 single-CpGs associated with the disease. Two of these showed high statistical power (99%) and effect size (>0.05) after false discovery rate (FDR) multiple comparison corrections: (i) cg03546163 in *FKBP5,* significantly hypomethylated in cases (Mean Difference in beta values (MD) = −0.09, 95% CI = −0.12|−0.06, *p* = 1.2 × 10−<sup>7</sup> ), and (ii) cg06633438 in *MLLT1*, statistically hypermethylated in cases (MD = 0.07, 95% CI = 0.04|0.10, *p* = 1.0 × 10−<sup>6</sup> ). Based on the interaction analysis, asbestos exposure and epigenetic profile together may improve MPM risk assessment. Above-median asbestos exposure and hypomethylation of cg03546163 in *FKBP5*

**Citation:** Cugliari, G.; Allione, A.; Russo, A.; Catalano, C.; Casalone, E.; Guarrera, S.; Grosso, F.; Ferrante, D.; Sculco, M.; La Vecchia, M.; et al. New DNA Methylation Signals for Malignant Pleural Mesothelioma Risk Assessment. *Cancers* **2021**, *13*, 2636. https://doi.org/10.3390/ cancers13112636

Academic Editors: Daniel L. Pouliquen and Joanna Kopecka

Received: 30 March 2021 Accepted: 24 May 2021 Published: 27 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

(OR = 20.84, 95% CI = 8.71|53.96, *p* = 5.5 × 10−11) and hypermethylation of cg06633438 in *MLLT1* (OR = 11.71, 95% CI = 4.97|29.64, *p* = 5.9 × 10−<sup>8</sup> ) genes compared to below-median asbestos exposure and hyper/hypomethylation of single-CpG DNAm, respectively. Receiver Operation Characteristics (ROC) for Case-Control Discrimination showed a significant increase in MPM discrimination when DNAm information was added in the model (baseline model, BM: asbestos exposure, age, gender and white blood cells); area under the curve, AUC = 0.75; BM + cg03546163 at *FKBP5*. AUC = 0.89, 2.1 × 10−<sup>7</sup> ; BM + cg06633438 at *MLLT1*. AUC = 0.89, 6.3 × 10−<sup>8</sup> . Validation and replication procedures, considering independent sample size and a different DNAm analysis technique, confirmed the observed associations. Our results suggest the potential application of DNAm profiles in blood to develop noninvasive tests for MPM risk assessment in asbestos-exposed subjects.

**Keywords:** malignant pleural mesothelioma; asbestos exposure; DNA methylation; epigenome-wide analysis; interaction analysis

#### **1. Introduction**

Mesothelioma has a long latency period, usually emerging 20–40 years after asbestos exposure [1]. Malignant pleural mesothelioma (MPM) is rare (prevalence 1–9/100,000), but the corresponding annual death toll worldwide is still estimated at about 40,000 [2,3]. Each year, 125 million people are exposed to asbestos, according to a World Health Organization report [4]. The International Agency for Research on Cancer confirmed that all fibrous forms of asbestos are carcinogenic to humans. The main outcome of exposure is mesothelioma, but cancer at other sites, such as respiratory-tract tumors, are moderately frequent [4]. Previous in vitro studies have demonstrated the cytotoxic effects of asbestos fibers [5,6].

A significant association between MPM and asbestos exposure has been reported, showing a clear, increasing trend in the odds ratio (OR) with increasing cumulative exposure among subjects exposed to over 10 fiber/mL-years [7]. Another study reported that the incidence of malignant mesothelioma (MM) was strongly associated with the proximity of one's residence to an asbestos exposure source [8].

DNA methylation (DNAm) is an epigenetic mechanism involved in gene expression regulation. In particular, dysregulation of promoter DNAm and histone modification are epigenetic mechanisms involved in human malignancies [9].

According to recent papers, both DNAm and genetic alterations may contribute to MPM tumorigenesis [10–15]. Whereas the genome remains consistent throughout one's lifetime, factors like ageing, lifestyle, environmental exposures and diseases can modify DNAm. The adaptive nature of DNAm means that it can be used to link environmental factors to the development of pathologic phenotypes such as cancers. Fasanelli et al. observed an association between exposure to tobacco and site-specific CpG methylation. They also used peripheral blood DNA to evaluate the importance of these epigenetic alterations in the aetiology of lung cancer [16].

There is less information on the mechanisms and clinical outcomes of epigenetic derangements in MPM [17–19]. Several studies have evaluated DNAm alterations in MM samples [20–22], but few of them focused on DNAm alteration in blood as a circulating marker. Fischer et al. examined serum DNAm of nine gene-specific promoters from MM cases [23]. A more recent paper identified hypomethylation of a single CpG in *FKBP5* in whole blood cells as a predictor of overall survival in MPM cases [13]. Guarrera et al. evaluated methylation levels in DNA from whole blood leukocytes as potential diagnostic markers for MPM and found a differential methylation between asbestos-exposed MPM cases and controls, mainly in genes related to the immune system [11]. The identification of reliable DNAm biomarkers with high sensitivity and specificity for MPM risk assessment would be a major advancement.

This study was undertaken with the goal to identify new biomarkers for MPM risk assessment and to determine if peripheral blood DNAm profiles have any predictive value. The second goal was to evaluate the interaction effect of asbestos exposure with DNAm on MPM risk. Currently, there are no sensitive testing methods available for the screening of asbestos-exposed individuals who are at high risk of developing MPM. Thus, the identification of reliable MPM diagnostic biomarkers in peripheral blood might provide a tool for detecting the disease at an early stage.

#### **2. Results**

#### *2.1. Epigenome-Wide Association Study (EWAS)*

CpGs (445,254) passed quality control procedures and were considered for statistical analyses. EWAS revealed two statistically significant differentially methylated single-CpGs between case and control groups: cg03546163 in the *FKBP5* gene (Mean Difference in beta values (MD) = 0.09, 95% CI = −0.12|−0.06, *p* = 1.2 × 10−<sup>7</sup> , *p* = 0.028) and cg06633438 in the *MLLT1* gene (MD = 0.07, 95% CI = 0.04|0.10, *p* = 1.0 × 10−<sup>6</sup> , *p* = 0.049) after False Discovery Rate (FDR) post hoc correction (Figure 1; Table 1). − − − − ≥

**Figure 1.** Manhattan plot for EWAS test on 450 k single CpGs. Single-CpG DNAm was used as dependent variable adjusting for age, gender, White blood cells (WBCs: monocytes, granulocytes, natural killer, B cells, CD4+ T and CD8+ T) estimation, population stratification and technical variability. FDR post hoc line highlights statistically significant differences between cases and controls at single CpG level.

Another 10 CpGs showed hypo/hypermethylation in MPM considering FDR < 0.05 but not effect size (MD) cut off ≥ |0.05| (Table 1).

− − − − − − Bootstrap was computed to estimate measures of accuracy using random sampling methods. The bias-corrected and accelerated (BCa) bootstrap interval was calculated for cg03546163 in *FKBP5* (95% CIBCa = −0.16|−0.10, z0 = −0.008, a = 0.002) and cg06633438 in *MLLT1* (95% CIBCa = −0.06|−0.1, z0 = −0.011, a = 0.0004) genes, confirming the robustness of the results considering the sample under study.

− − − − − − −

− <sup>−</sup>

₼

− −

− − − −

−


**Table 1.** Differential DNAm analyses ordered by effect size. Information about single-CpGs, including location-related values and model outputs (effect size, standard error, *p* values). −

− − − −

− − − −

−

−

− − <sup>−</sup> − − Control group was set as reference. Adjustment covariates: age, gender, population stratification, WBCs (monocytes, granulocytes, natural killer, B cells, CD4+ T and CD8+ T) estimation and technical variability. \*: statistically significant at *p* value < 0.05; §: statistically significant at FDR post hoc adjustments. − <sup>−</sup> \*§₼ : statistically significant at beta = 0.01.

> ₼ ₼ Statistically significant differences in MD between cases and controls were found in the WBCs estimated (monocytes, *p* = 6.0 × 10−<sup>3</sup> ; granulocytes, *p* = 2.2 × 10−16; B cells, *p* = 1.1 × 10−12; NK cells, *p* = 3.6 × 10−<sup>4</sup> ; CD4+ T, *p* = 2.2 × 10−16; CD8+ T, *p* = 6.8 × 10−11; Naïve CD4T, *p* = 0.012; Naïve CD8T, *p* = 7.0 × 10−<sup>3</sup> ).

− <sup>−</sup>

− −

− −

In order to assess if smoking status, classified as current, former and never-smokers, could modify DNAm profiles, we performed a multivariate regression analysis with the same model used for the discovery phase. No evidence of methylation differences linked to different smoking levels was found for any of the twelve statistically significant CpGs.

#### *2.2. Receiver Operation Characteristics (ROC) for Case-Control Discrimination*

− The baseline model (BM) including age, gender, asbestos exposure and WBCs was compared with BM adding the DNAm levels of cg03546163 or cg06633438. Receiver Operation Characteristics 8ROC9 curves showed a significant increase in MPM discrimination when DNAm information was added in the model (Table 2).

− − − **Table 2.** Disease discrimination test considering (AUC) comparison between baseline model and models additionally including single-CpG.


− Models are shown as baseline model (BM) or BM + Single CpG DNAm. AUC Differences between considered model and BM were estimated with the DeLong's test.

#### *2.3. Interaction Analysis*

−

CpG sites and asbestos exposure were considered as predictors of MPM risk in the interaction model. Categorical variables (quantile information) were used considering median values.

We tested the interaction between asbestos exposure and DNAm levels at cg03546163 in *FKBP5* and cg06633438 in *MLLT1.*

Considering cg03546163 in *FKBP5,* DNA hypermethylation and low asbestos exposure levels were used as references, while for cg06633438 in *MLLT1,* DNA hypomethylation and low asbestos exposure levels were set as references (Table 3).

The OR was estimated as the relationship between the combination of single-CpGs DNAm levels and asbestos exposure quantile, and the reference (low median asbestos exposure and hypermethylation status for cg03546163, or hypomethylation status for cg06633438). Age, gender, population stratification, and WBCs were included in the GLM (family = binomial) to adjust the interaction effect.


**Table 3.** Interaction between asbestos exposure and single CpG DNAm on the MPM Odds ratios.

Reference for cg03546163 in *FKBP5*: hypermethylation and low asbestos exposure levels; Reference for cg06633438 in *MLLT1*: hypomethylation and low asbestos exposure levels.

> The relationship between asbestos exposures and single-CpG DNAm levels was evaluated. An increase of one unit of asbestos exposure (rank transformed fibers/mL years) was related to the *FKBP5* gene (β = −0.016, 95% CI = −0.031|−0.001, *p* = 0.044) and *MLLT1* gene (β = −0.014, 95% CI = 0.001|0.026, *p* = 0.035) methylation level variations.

> Strong association between asbestos exposure and MPM risk, considering dichotomous distribution of asbestos exposure, was found (OR = 6.11, 95% CI = 3.73|10.20, *p* = 1.8 × 10−12). Quartile distribution of asbestos exposure was evaluated to estimate the potential incremental association with MPM risk (1st quartile: used as reference; 2nd quartile: OR = 1.83, 95% CI = 0.93|3.69, *p* = 0.09; 3rd quartile: OR = 6.63, 95% CI = 3.30|13.81, *p* = 2.1 × 10−<sup>7</sup> ; 4th quartile: OR = 11.00, 95% CI = 5.26|24.30, *p* = 7.3 × 10−10).

#### *2.4. Validation and Replication*

For the replication and validation approaches, an independent sample of 140 MPM cases and 104 cancer-free asbestos-exposed controls from the same areas were considered, using a targeted DNAm analysis technique.

The direction and magnitude of the association was consistent for cg03546163 and cg06633438 DNAm. Patients showed significantly lower DNAm at cg03546163 (MD = −0.061, 95% CI = −0.087|−0.036, *p* = 4.5 × 10−<sup>6</sup> ) and higher DNAm at cg06633438 (MD = 0.024, 95% CI = 0.061|0.013, *p* = 4.0 × 10−<sup>2</sup> ) compared with controls. A multivariate analysis confirmed that DNAm at cg03546163 in *FKBP5* and cg06633438 in *MLLT1* were independently associated with MPM detection.

#### **3. Discussion**

In the present study, we used a whole genome microarray approach to investigate DNAm in WBCs from MPM cases and asbestos-exposed cancer-free controls from a region with a history of asbestos exposure (Piedmont, Italy) [10] in order to identify new noninvasive epigenetic markers related to MPM. The identification of reliable MPM diagnostic biomarkers in peripheral blood might improve risk assessment.

We observed hypomethylation of CpG cg03546163, located in the 5′ UTR region of *FKBP5* gene, in MPM cases compared to controls.

Epigenetic activation of the FKBP Prolyl Isomerase 5 (*FKBP5*) gene has been shown to be associated with increased stress sensitivity and the risk of psychiatric disorders [24]. *FKBP5* is an immunophilin and has an important role in immunoregulation and in protein folding and trafficking. It plays a role in transcriptional complexes and acts as a cotranscription factor, along with other proteins in the *FKBP* family [25]. The suggestion of a possible role of *FKBP5* in the development and progression of different types of cancer has stemmed from several studies. In particular, high protein expression has been linked to either suppression or promotion of tumour growth, depending on tumour type and microenvironment [26,27].

*FKBP5* is involved in the NF-kB and AKT signaling pathways, both of which are implicated in tumorigenesis [28]. Notably, NF-kB appears to be frequently constitutively activated in malignant tumours and involved in the modulation of genes linked to cell motility, neoangiogenesis, proliferation and programmed cell death [29]. An epigenetic upregulation of *FKBP5* could promote NF-kB activation [30]. STAT3-NFkB activity is involved in chemoresistance in MM cells [31], and NFkB was shown to be constitutively active as a result of asbestos-induced chronic inflammation [32].

CpG cg06633438 located in the body region of the *MLLT1* gene was hypermethylated in cases compared to controls.

The *MLLT1* gene encodes the ENL protein, a histone acetylation reader component of the super elongation complex (SEC), which promotes transcription at the elongation stage by suppressing transient pausing by the polymerase at multiple sites along the DNA. In acute myeloid leukemia, *MLLT1* regulates chromatin remodeling and gene expression of many important proto-oncogenes [31]. Yoshikawa and colleagues suggested that mesothelioma may be the consequence of the somatic inactivation of chromatin-remodeling complexes and/or histone modifiers, including *MLLT1* [30].

In mesothelioma patients with short-term recurrence after surgery, frequent 19p13.2 loss was reported. This region encompasses several putative tumor suppressors or oncogenes, including *MLLT1* [32].

Interestingly, *MLLT1* and *FKBP5* showed opposite behavior, increasing and decreasing DNAm levels, respectively, in relation to MPM. This finding could reflect the opposite expression profiles of the two genes among all the different subtypes of white blood cells in normal human hematopoiesis, as reported in the Blood Spot database (http://servers. binf.ku.dk/bloodspot/, accessed on 26 May 2021) (Figure 2) [33].

**Figure 2.** Expression profiles in normal human haematopoiesis. *MLLT1* (**A**) and *FKBP5* (**B**) expression profiles in normal human haematopoiesis as reported in the Blood Spot database (http://servers.binf.ku.dk/bloodspot/, accessed on 26 May 2021).

Our interaction analysis showed that considering DNAm levels at *FKBP5* and *MLLT1* genes together with asbestos exposure levels may help to better define MPM risk for asbestos-exposed subjects.

Six single-CpGs showed differential methylation in patients, including those located in *C13orf30, CDC42BPB, ZNF629* and *ALG9* genes; the other six were not annotated to named genes. *ALG9* is a glycogene whose reduced expression has been described during the epithelial-to-mesenchymal transition, an essential process also involved in cancer progression [34]. The *CDC42BPB* gene is ubiquitously expressed in mammals and encodes a

serine/threonine protein kinase, a member of the MRCK family [35]. The role of MRCKs in cytoskeletal reorganization during cell migration and invasion has been characterized [36]. The biological function of *C13orf30* and *ZNF629*, a DNA-binding transcription factor, is still to be established.

MPM cases and asbestos-exposed controls showed different proportions of estimated WBCs, which may denote the crucial implication of the immune system. It is known that in cancer, including mesothelioma, the immune system is affected [37], and there is evidence that asbestos directs antigen overstimulation, and that reactive oxygen species production induces functional changes in WBCs [38]. Indeed, in MPM cases, we showed a reduction of estimated CD4+ and CD8+ T lymphocytes, suggesting a weaker adaptive immune system [39]. This may reflect the possible occurrence of functional changes in WBC subtypes in MPM [40,41].

The need for reliable biomarkers is of extreme relevance for a disease such as MPM, which is characterized by the accumulation and persistence of asbestos fibers in the lungs, leading to a long latency period before clear clinical signs of the tumor are detectable. Several biomarkers for early MPM detection (e.g., mesothelin, osteopontin and fibulin-3) have been proposed so far; however, some of them are still under investigation [42]. In this context, DNAm changes in easily-accessible WBCs may provide a useful tool to better assess MPM risk in asbestos-exposed subjects.

Our findings that DNAm levels in single-CpGs in *FKBP5* and *MLLT1* genes are independent markers of MPM in asbestos-exposed subjects suggest the potential use of blood DNAm analysis as a noninvasive test for MPM detection.

Some somatic gene alterations in lung cancer have been linked to tobacco smoke, but few data are available on the role of asbestos fibers: Andujar and colleagues investigate the mechanism of P16/CDKN2A alterations in lung cancer including asbestos-exposed patients. P16/CDKN2A gene inactivation in asbestos-exposed non-small-cell lung carcinoma (NSCLC) cases, a tumor independent of tobacco smoking but associated with asbestos exposure, mainly occurs via promoter hypermethylation, loss of heterozygosity and homozygous deletion, suggesting a possible relationship with an effect of asbestos fibers [43].We observed epigenetic deregulations in the blood of MPM patients compared to that of cancer-free controls, suggesting the potential use of DNAm for risk stratification among asbestos-exposed individuals.

If this observation can be verified in prospectively collected samples, it may be possible to use CpGs methylation to further improve MPM risk estimation for subjects with occupational and/or environmental asbestos exposure.

#### *Limitation of the Study*

Leukocyte DNAm may be a nonspecific marker related to a general, tumour-induced inflammatory status rather than a specific MPM biomarker. Further studies are therefore needed to support our findings.

One main limitation of the functional interpretation of our results is that all our cases had already developed MPM at recruitment: thus, our findings likely reflect disease status rather than being markers of the dynamic processes leading to MPM onset. The lack of MPM tissue from the same subjects also poses major constraints to the functional interpretation of our findings.

Notwithstanding the above limitations, the discrimination between MPM cases and asbestos-exposed cancer-free controls improved when DNAm levels were considered together with asbestos exposure levels.

#### **4. Material and Methods**

#### *4.1. Study Population*

Study subjects were part of a wider, ongoing collaborative study, which is actively enrolling MPM cases and cancer-free controls in the municipality of Casale Monferrato (Piedmont Region, Italy). This area was chosen due to its exceptionally high incidence of

mesothelioma, caused by widespread occupational and environmental asbestos exposure originating from the Eternit asbestos-cement plant, which was operational until 1986 [44]. Additional MPM cases and cancer-free controls were recruited from other main hospitals of the Piedmont Region (in the municipalities of Turin, Novara and Alessandria). The ongoing collaborative study includes MPM cases diagnosed between incident MPM cases diagnosed between 2000 and 2010 after histological and/or cytological confirmation, and matched controls [45].

The present study included 159 MPM cases and 137 cancer-free controls from a larger case-control study, all of whom had genetic and blood DNAm data [46], good quality DNA at the time of the analyses, and information on asbestos exposure above the background level, as defined in Ferrante et al. [47]. MPM cases and asbestos-exposed cancer-free controls were matched by date of birth (±18 months) and gender. An additional 244 (140 MPM cases and 104 cancer-free controls) independent samples from the same case-control study were included for validation/replication analyses.

Tables 4 and 5 shows the descriptive characteristics of controls and cases (Min, 1st Q, Median, Mean, 3rd Q and Max) that were considered in the statistical analysis (gender, age, asbestos exposure and WBC estimates: monocytes, granulocytes, natural killer, B cells, CD4+ T and CD8+ T). Asbestos exposure (occupational, environmental and domestic) was normalized considering frequency, duration and intensity. Smoking status (current, former and never smokers) is also explained in Table 6.

**Table 4.** Descriptive characteristics of cancer-free control group.


Minimum (Min), First Quartile (1st Q), Median, Mean, Third Quartile (3rt Q) and Maximum (Max) of variables related to cancer-free controls.



Minimum (Min), First Quartile (1st Q), Median, Mean, Third Quartile (3rt Q) and Maximum (Max) of variables related to MPM cases.


**Table 6.** Descriptive characteristics of smoking status stratified by disease.

*n* and % of the three levels of smoking status stratified by disease.

Our study complied with the Declaration of Helsinki principles and conformed to ethical requirements. All volunteers signed an informed consent form at enrollment. The study protocol was approved by the Ethics Committee of the Italian Institute for Genomic Medicine (IIGM, Candiolo, Italy).

#### *4.2. Exposure Assessment*

Information on occupational history and lifestyle habits were collected from all subjects through interviewer-administered questionnaires, which were completed during face-to-face interviews at enrollment. Job titles were coded in two ways according to the International Standard Classification of Occupations [47] and the Statistical Classification of Economic Activities in the European Community.

A cumulative exposure index was computed considering frequency, duration and intensity of asbestos exposure. Occupational, environmental and domestic asbestos exposure were evaluated by an experienced occupational epidemiologist [47], and exposure doses (fibers/mL years) were rank-transformed to remove skewness.

#### *4.3. Blood DNAm Analysis and Beta-Value Extraction*

DNAm levels were measured in DNA from whole blood collected at enrollment using the Infinium HumanMethylation450 BeadChip (Illumina, San Diego, CA, USA). For blood DNAm analysis (including quality control) please refer to the previous work of the same group [11].

We used the R package 'methylumi' to analyze DNAm data. The average methylation value at each locus was computed as the ratio of the intensity of the methylated signal over the total signal (unmethylated + methylated) [48]. Beta-values ranging from 0 (no methylation) to 1 (full methylation) represent the percentage of methylation at each individual CpG locus.

We excluded the following from the analyses: (i) single beta-values with a *p*-value for detection ≥ 0.01; (ii) CpG loci that had missing beta-values in more than 20% of the assayed samples; (iii) CpG loci detected by probes containing single nucleotide polymorphisms (SNPs) with MAF ≥ 0.05 in the CEPH (Utah residents with ancestry from northern and western Europe, CEU) population; and iv) samples with a global call rate ≤ 95%. We also excluded CpGs on chromosomes X and Y.

#### *4.4. Batch Effect, Population Stratification and White Blood Cells Estimations*

All differential methylation analyses were corrected for "control probes" Principal Components (PCs) to account for variability and batch effects in methylation assays. We used PCs assessed by principal component analysis of the BeadChip's built-in control probes as a correction factor for statistical analyses of microarray data. This method allows researchers to account for the technical variability in the different steps in DNAm analysis, from bisulfite conversion to BeadChip processing [49].

An individual's geographic origins may influence DNAm profiles, which could potentially introduce bias. To take this into consideration, we took advantage of the available data from our previous study, which includes a genome-wide genotyping dataset from the same study subjects [50]. When genome-wide genotyping was used to calculate the first PCs, they were shown to correlate with different geographic origins [51].

For each subject, we extracted WBC subtype percentages, estimated based on genomewide methylation data. This method provides quantification of the composition of leukocytes than can be achieved by simple histological or flow cytometric assessments, with an admissible range of variability [52].

#### *4.5. Statistical Analyses*

#### 4.5.1. Epigenome-Wide Association Study

An association test was used to analyze the mean differences (MD) in single-CpG methylation between MPM cases and asbestos-exposed cancer-free controls. We performed multiple regression analysis adjusted for age, gender, estimated WBCs (monocytes, granulocytes, natural killer, B cells, CD4+ T and CD8+ T), population stratification (first 2 PCs) and technical variability (first 10 PCs). For multiple comparison tests, a FDR *p* value ≤ 0.05 was considered statistically significant.

Bootstrapping was performed using random sampling methods to estimate the measures of accuracy defined in terms of bias, variance, confidence intervals and prediction error. Bootstrapping can also be applied to control and check the results for stability. The bias-corrected and accelerated (BCa) bootstrap interval was calculated with regard to single CpGs.

ROC for Case-Control Discrimination was implemented, and the AUC metric was applied to estimate the predictive performance of a binary classification (cases/controls). The baseline model (BM) included age, gender, asbestos exposure and WBCs, and was compared with the BM after adding the DNAm levels of statistically significant, single-CpGs at EWAS. AUC differences between BMs before and after the addition of DNAm levels were estimated with DeLong's test.

#### 4.5.2. Statistical Power

To ensure a study power greater than 99% (two-tailed test at α = 0.05 and β = 0.01), only CpGs with a MD between cases and controls ≥ |0.05| were selected.

Covariates were included step-by-step in a sensitivity analysis to validate the association output considering effect size, standard error, 95% confidence interval and *p* value variations.

Gene set enrichment analyses were carried out on CpGs with a False Discovery Rate *p* value (PFDR) ≤ 0.05 to identify pathways that may be affected by MPM-related changes in methylation.

All statistical analyses were conducted using the open source software R (4.0.2).

#### 4.5.3. Interaction Analysis

Logistic regression was used to analyze the relationship between CpGs and asbestos exposure in MPM risk (odds ratio), adjusting for age, gender, SNP PCs and WBCs estimates. Asbestos exposure was classified as above-median or below-median, and CpG methylation was categorized as above-median or below-median.

MPM risk for a given CpG level and asbestos exposure was expressed by ORij, where i indicates the asbestos exposure (below-median or above-median) and j indicates the CpG (above-median or below-median). Considering the direction of the effect, the same approach was used: for hypomethylated CpGs, above-median was used as the reference level, whereas below-median was used for hypermethylated CpGs.

Subjects with below-median asbestos exposure and reference-level CpG DNAm were considered the baseline group, and their MPM risk was coded as OR00 = 1. Interaction was analyzed with respect to both additive and multiplicative models based on the ORs obtained by logistic regression.

Synergistic interaction (positive interaction) implies that the combined action of two factors in an additive model is greater than the sum of their individual effects. Antagonistic interaction, on the other hand, means that when two factors are present in an additive model, the action of one reduces the effect of the other.

Multivariable logistic regression models were used to explore any deviations from a multiplicative model, including asbestos exposure, CpG and the corresponding interaction term (CpG × exposure). All models were adjusted for age, gender, SNP PCs, technical covariates and WBCs estimates. *p*-values < 0.05 were considered statistically significant.

#### *4.6. Validation and Replication*

DNAm signal validation and replication was done by the EpiTYPER MassARRAY assay (Agena Bioscience). This assay uses a MALDI-TOF mass spectrometry-based method to quantitatively assess the DNA methylation state of the CpG sites of interest [53]. DNA (500 ng) was bisulfite-converted as indicated in Section 4.3.

PCR amplification, treatment with SAP solution and Transcription/RNase A cocktails were performed according to the manufacturer's instructions, and the mass spectra were analyzed by an EpiTYPER analyzer. The MassARRAY assay cannot discriminate between CpGs located in close proximity in the sequence, so instead, the close neighboring CpGs are analyzed as "Units", i.e., the measured methylation level is the average of the methylation levels of the CpGs cumulatively analyzed within the Unit. In the case of cg03546163, the measured methylation level is the average between two CpG sites located in very close proximity (Figure S1). For cg06633438, the two adjacent signals were considered, since the results for the model did not differ for effect size, standard error, 95% CI or *p* value (Figure S2).

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/cancers13112636/s1, Figure S1: Location of cg03546163 investigated by EpiTYPER MassARRAY, Figure S2: Location of cg06633438 investigated by EpiTYPER MassARRAY.

**Author Contributions:** Conceptualization, G.C., G.M., C.M., D.M., I.D. and F.G.; methodology, G.C.; software, G.C.; validation, G.C., A.A. and C.P.; formal analysis, G.C.; investigation, G.C., G.M., S.G., A.A., A.R., E.C. and C.C.; resources, G.C., G.M., C.M., D.M., I.D., F.G., D.F., M.S., M.L.V. and R.L.; data curation, G.C., C.M., D.M., I.D., F.G., C.P., D.F., M.S., M.L.V. and R.L.; writing—original draft preparation, G.C., S.G., A.A., A.R., E.C. and C.C.; writing—review and editing, G.C., G.M., C.M., D.M., I.D., F.G., S.G., A.A., A.R., E.C., C.C., C.P., D.F., M.S., M.L.V. and R.L.; visualization, G.C., G.M., S.G., A.A., A.R., E.C., C.C., C.P., D.F., M.S., M.L.V. and R.L.; supervision, G.C. and G.M.; project administration, G.C. and G.M.; funding acquisition, G.M., C.M., D.M., I.D. and F.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research leading to these results received funding from AIRC under IG 2018-ID. 21390 project (to G.M.), from the HERMES (Hereditary Risk in MESothelioma) Project funded by the offer of compensation to the inhabitants of Casale Monferrato deceased or affected by mesothelioma (to I.D. and C.M.), and from the Ministero dell'Istruzione, dell'Università e della Ricerca–MIUR project "Dipartimenti di Eccellenza 2018–2022" (n◦ D15D18000410001, to G.M.) to the Department of Medical Sciences, University of Torino.

**Institutional Review Board Statement:** Ethics approval and consent to participate: The study protocol was approved by the Ethics Committee of the Italian Institute for Genomic Medicine (IIGM, Candiolo, Italy).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The methylation of single individuals cannot be published due to informed consent limitations.

**Acknowledgments:** The authors wish to thank all the patients contributing to this research, Caterina Casadio and Ezio Piccolini for their past clinical contribution and Paolo Garagnani for the Sequenom analysis supervision.

**Conflicts of Interest:** The authors declare that they have no competing interest.

#### **References**


*Article*
