New DNA Methylation Signals for Malignant Pleural Mesothelioma Risk Assessment

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 asbestos-exposed 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−7), and (ii) cg06633438 in MLLT1, statistically hypermethylated in cases (MD = 0.07, 95% CI = 0.04|0.10, p = 1.0 × 10−6). 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 (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−8) 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−7; BM + cg06633438 at MLLT1. AUC = 0.89, 6.3 × 10−8. 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.


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][11][12][13][14][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][18][19]. Several studies have evaluated DNAm alterations in MM samples [20][21][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.
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.

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.
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. 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.
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. 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.

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).

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.  . 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. 0.049 * § 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. ₼: statistically significant at beta = 0.01.

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).

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. considered as predictors of MPM risk in the uantile information) were used considering tos exposure and DNAm levels at cg03546163 : 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 −3 ; granulocytes, p = 2.2 × 10 −16 ; B cells, 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.

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).

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. Reference for cg03546163 in FKBP5: hypermethylation and low asbestos exposure levels; Reference for cg06633438 in MLLT1: hypomethylation and low asbestos exposure levels.

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.

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].
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].  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 Cancers 2021, 13, 2636 7 of 14 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.

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 Cancers 2021, 13, 2636 8 of 14 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.   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).

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.

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.

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].

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.

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 (P FDR ) ≤ 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).

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.

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. 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.