Genetic and Epigenetic Characterization of a Discordant KMT2A/AFF1-Rearranged Infant Monozygotic Twin Pair

The KMT2A/AFF1 rearrangement is associated with an unfavorable prognosis in infant acute lymphocytic leukemia (ALL). Discordant ALL in monozygotic twins is uncommon and represents an attractive resource to evaluate intrauterine environment–genetic interplay in ALL. Mutational and epigenetic profiles were characterized for a discordant KMT2A/AFF1-rearranged infant monozygotic twin pair and their parents, and they were compared to three independent KMT2A/AFF1-positive ALL infants, in which the DNA methylation and gene expression profiles were investigated. A de novo Q61H NRAS mutation was detected in the affected twin at diagnosis and backtracked in both twins at birth. The KMT2A/AFF1 rearrangement was absent at birth in both twins. Genetic analyses conducted at birth gave more insights into the timing of the mutation hit. We identified correlations between DNA methylation and gene expression changes for 32 genes in the three independent affected versus remitted patients. The strongest correlations were observed for the RAB32, PDK4, CXCL3, RANBP17, and MACROD2 genes. This epigenetic signature could be a putative target for the development of novel epigenetic-based therapies and could help in explaining the molecular mechanisms characterizing ALL infants with KMT2A/AFF1 fusions.


Introduction
Acute lymphoblastic leukemia (ALL) represents 25% of all diagnosed tumors among 0-14-year-old patients [1]. Leukemia in children less than one year old (infants) represents about 3% of all pediatric ALL and most frequently develops a B-cell ALL form, which is characterized by an uncontrolled proliferation of immature B-cell precursors in the peripheral blood (PB) and bone marrow [2].
The triggering genetic alteration underlying leukemogenesis is a chromosomal translocation occurring at the KMT2A (formerly MLL) gene locus (11q23), with the KMT2A/AFF1 fusion being the most frequent cytogenetic abnormality in infant B-cell ALL. For children in

Exome Sequencing on Monozygotic Twins Discordant for ALL
To identify ALL-linked de novo somatic mutations, the siblings and parents were subjected to WES. No high-confidence de novo somatic indel changes were detected in the proband. Among the de novo somatic variants of the proband with global allele frequencies lower than 1%, we detected the missense heterozygous mutation NRAS c.183A > T p.Q61H (variant allele frequency = 26.7%) as being absent in all the healthy samples and confirmed this using Sanger sequencing (Supplementary Table S1) ( Figure 1B).
In order to determine whether the NRAS mutation was already present at birth or acquired afterward, neonatal blood spots of the twins were scrutinized to search for NRAS p.Q61H mutation by ddPCR, and 0.13 copies/µL and 0.17 copies/µL of a mutant allele were detected in both the proband and the healthy twin sister, respectively. In order to determine whether the NRAS mutation was already present at birth or acquired afterward, neonatal blood spots of the twins were scrutinized to search for NRAS p.Q61H mutation by ddPCR, and 0.13 copies/µL and 0.17 copies/µL of a mutant allele were detected in both the proband and the healthy twin sister, respectively.

Differential DNA Methylation and Gene Expression
In all five family samples (S1-S5), a panel of 393,797 CpG was obtained. In order to specifically identify ALL-related DNA methylation changes in the affected child (S1), we focused our analysis on 8120 CpGs, satisfying the following criteria: ∆β-values of at least 30% between the proband (S1) and both her healthy sister (S2) and the remitted sample (S3) and ∆β-values ≤10% among the healthy samples (S2 and S3). Most of the CpGs (6599) were consistently hypermethylated in the proband at diagnosis compared with the healthy samples (∆β-values between 30% and 82%). Among the 8120 CpG sites identified within the twin samples, 5139 were also differentially methylated with similar effects (i.e., t-test mean differences of at least 30%) and the same directions among the three diagnosis-remission-matched samples (couples A1d-A1r, A2d-A2r, and A3d-A3r) which were derived from three independent (and unrelated) B-ALL infants carrying the KMT2A/AFF1 chromosomal rearrangement (Supplementary Table S2). The majority of these CpGs (4396) were hypermethylated at diagnosis (30-85% difference) in accordance with what was observed within the twin samples.
In total, the 5139 differentially methylated CpG sites were annotated to 1904 gene regions. The top 20 GSEA results are listed in Supplementary Table S3. All samples collected at diagnosis (S1, A1d, A2d, and A3d) consisted predominantly of immature B-cells, whereas the remission samples (S3, A1r, A2r, and A3r) included a mixture of different mature blood cell populations. To take this difference into account,

Differential DNA Methylation and Gene Expression
In all five family samples (S1-S5), a panel of 393,797 CpG was obtained. In order to specifically identify ALL-related DNA methylation changes in the affected child (S1), we focused our analysis on 8120 CpGs, satisfying the following criteria: ∆βvalues of at least 30% between the proband (S1) and both her healthy sister (S2) and the remitted sample (S3) and ∆β-values ≤10% among the healthy samples (S2 and S3). Most of the CpGs (6599) were consistently hypermethylated in the proband at diagnosis compared with the healthy samples (∆β-values between 30% and 82%). Among the 8120 CpG sites identified within the twin samples, 5139 were also differentially methylated with similar effects (i.e., t-test mean differences of at least 30%) and the same directions among the three diagnosis-remission-matched samples (couples A1 d -A1 r , A2 d -A2 r , and A3 d -A3 r ) which were derived from three independent (and unrelated) B-ALL infants carrying the KMT2A/AFF1 chromosomal rearrangement (Supplementary Table S2). The majority of these CpGs (4396) were hypermethylated at diagnosis (30-85% difference) in accordance with what was observed within the twin samples.
In total, the 5139 differentially methylated CpG sites were annotated to 1904 gene regions. The top 20 GSEA results are listed in Supplementary Table S3. All samples collected at diagnosis (S1, A1 d , A2 d , and A3 d ) consisted predominantly of immature B-cells, whereas the remission samples (S3, A1 r , A2 r , and A3 r ) included a mixture of different mature blood cell populations. To take this difference into account, we subjected the entire methylation dataset shared by the twins (S1, S2, and S3), the three diagnosis-remission-paired samples (couples A1 d -A1 r , A2 d -A2 r , and A3 d -A3 r ) and the CD34 + hematopoietic stem cells (398,224 CpG sites) to unsupervised clustering analysis. As shown in Figure 2, the first separation clearly divided the affected samples from the healthy samples.
we subjected the entire methylation dataset shared by the twins (S1, S2, and S3), the three diagnosis-remission-paired samples (couples A1d-A1r, A2d-A2r, and A3d-A3r) and the CD34 + hematopoietic stem cells (398,224 CpG sites) to unsupervised clustering analysis. As shown in Figure 2, the first separation clearly divided the affected samples from the healthy samples. To verify that the "top" 5139 differentially methylated CpG sites corresponded to a signature characteristic of B-ALL, independent from the different proportion of blood cells in the samples, we compared the CpG site methylation levels among the leukemic samples (S1, A1d, A2d, and A3d) and non-leukemic samples (S2, S3, A1r, A2r, and A3r) and the non-leukemic CD19 + B cells and CD34 + cells isolated from the adult donors and cord blood samples, respectively. This set of non-leukemic reference cells represented mature B-cells (CD19 + ) and multipotent progenitor cells (CD34 + ). The "targeted" PCA ( Figure 3) revealed a cluster of B-ALL-affected individuals distinct from the non-leukemic individuals and reference samples. PC2 separated the CD19 + cells from all the other samples. We imputed this result to the different ages between the CD19 + cell donors (adults) versus the infant samples. To verify that the "top" 5139 differentially methylated CpG sites corresponded to a signature characteristic of B-ALL, independent from the different proportion of blood cells in the samples, we compared the CpG site methylation levels among the leukemic samples (S1, A1 d , A2 d , and A3 d ) and non-leukemic samples (S2, S3, A1 r , A2 r , and A3 r ) and the non-leukemic CD19 + B cells and CD34 + cells isolated from the adult donors and cord blood samples, respectively. This set of non-leukemic reference cells represented mature B-cells (CD19 + ) and multipotent progenitor cells (CD34 + ). The "targeted" PCA ( Figure 3) revealed a cluster of B-ALL-affected individuals distinct from the non-leukemic individuals and reference samples. PC2 separated the CD19 + cells from all the other samples. We imputed this result to the different ages between the CD19 + cell donors (adults) versus the infant samples.  Whole transcriptome profiles were determined for the four diagnosis-remission samples through mRNA sequencing. Unsupervised hierarchical clustering revealed expression profiles distinguishing the diagnosis samples from the remission samples (Supplementary Figure S2). Specifically, we identified 4376 differentially expressed genes in Whole transcriptome profiles were determined for the four diagnosis-remission samples through mRNA sequencing. Unsupervised hierarchical clustering revealed expression profiles distinguishing the diagnosis samples from the remission samples (Supplementary Figure S2). Specifically, we identified 4376 differentially expressed genes in the diagnosis versus remission samples (1.44 < |Log2FoldChange (FC)| < 12.34; 6.03 × 10 −42 < PFDR < 0.05; Supplementary Table S4). The downregulated (2536) and upregulated (1840) genes at diagnosis included genes involved in the development or functioning of the immune system, response to stimuli, and characterization of the immune cell subtypes. The top 20 significant overlaps computed for both the down-and upregulated genes are provided in Supplementary Table S5A,B.
Among the 4376 genes significantly up-and downregulated at diagnosis, 363 (8.3%) overlapped those harboring the 5139 differentially methylated CpG sites in both of the twin samples and in the three paired diagnosis-remission samples. We correlated the CpG methylation changes and expression levels. Thirty-two out of 363 genes harbored at least 3 CpGs (134 CpG sites in total; Supplementary Table S6), with rho ≥ |0.7|. Notably, in the healthy parents (S4 and S5), the methylation profiles of these 134 CpG sites were comparable to those of an internal control group composed of healthy adult subjects (data not shown), which excluded possible deregulated patterns potentially transmitted from parents to their daughters.
GSEA analysis performed on the 32 genes whose expression levels correlated with the DNA methylation profiles showed significant overrepresentation of genes upregulated in embryonic stem cells from TCEB3 knock-out mice, cell-cell junctions, and cell movements (FDR q-value < 9.36 × 10 −3 ; Supplementary Table S7). The differentially methylated CpG sites corresponding to these 32 genes were all hypermethylated at diagnosis, and 65% of them were within the gene promoter. Negative correlations were observed between methylation and gene expression for all CpGs except for those within the MACROD2 gene and one CpG in the PLEKHG5 3 UTR. The expression levels of five genes showed the strongest significance between the affected and healthy subjects (PFDR < 10 −6 ) as well as a correlation higher than 80% with corresponding DNA methylation changes. In particular, the methylation levels of 11 CpG sites located within 1500bp or 200bp upstream of the transcription starting site (TSS) of the RAB32 gene were inversely correlated with their expression levels (log2FC = −4.42; PFDR = 4.08 × 10 −8 ; −0.97 ≤ rho ≤ −0.84; Pcorr ≤ 0.04; Figure 4).

Discussion
In the current study, to identify new potential molecular mechanisms characterizing KMT2A/AFF1 + B-ALL, we investigated, for the first time, the case of monozygotic twins with the discordant diagnosis of KMT2A/AFF1 leukemia by exploiting different "-omics" approaches. Exome sequencing showed only the de novo NRAS Q61H mutation. Whole-

Discussion
In the current study, to identify new potential molecular mechanisms characterizing KMT2A/AFF1 + B-ALL, we investigated, for the first time, the case of monozygotic twins with the discordant diagnosis of KMT2A/AFF1 leukemia by exploiting different "-omics" approaches. Exome sequencing showed only the de novo NRAS Q61H mutation. Wholegenome sequencing has previously revealed a mutational landscape characterized by few non-silent alterations in ALL infants with rearrangements of the KMT2A gene [7,8]. This is consistent with the very short latency period between KMT2A rearrangement onset (often occurring in utero) and the clinically overt leukemia [9,10]. The NRAS Q61H mutant allele frequency well below 50% in the proband confirms the subclonal nature of the mutation [7]. A high frequency of RAS mutations has been reported in infant patients, supporting the hypothesis that RAS activation may be responsible for the shortened latency [11]. Moreover, NRAS mutations seem to contribute to the adverse prognosis typical of KMT2A-rearranged infant patients [12][13][14][15].
In our study, the NRAS mutation was also backtracked in the siblings' neonatal blood spots, suggesting that the mutation event arose during the prenatal period. In contrast, no pre-leukemic cells with the KMT2A/AFF1 rearrangement were detected in either twins at birth. These results offer two hypothetical scenarios: (1) a post-natal chromosomal rearrangement event, which would indicate that the NRAS subclonal mutation preceded the chromosomal aberration, a phenomenon recently described in an infant with acute myelomonocytic leukemia [16], or (2) undetectable rearranged blasts could have been present in both the neonatal twins. Considering the fact that most ALL infants could have less than one pre-leukemic blast of over 104 nucleated cells at birth [17], and we did not even observe the patient-specific clonal IGH VH6-JH3 rearrangement, we could not exclude the presence of a silent pre-leukemic clone. We could speculate the existence of a clone with the KMT2A/AFF1 fusion, which competes with the NRAS-mutated clone and initiates the leukemic process, and a minor, which may accelerate the development of the overt leukemia.
The discordancy for leukemia in monozygotic twins is very rare [6,[18][19][20]. Placental anastomoses together with the early onset of the disease support a prenatal origin of leukemia, with a concordance rate among identical twins reaching 100% [6]. Three hypotheses could explain the discordancy for ALL: (1) The pre-leukemic clones remain in a dormant state in the healthy child, as it was shown that, occasionally, pre-leukemic clones can last up to 14 years [21]. (2) The clearance of the pre-leukemic clones happened within the healthy infant. This speculation is based on the adrenal hypothesis that proposes an anti-leukemic effect as a consequence of cortisol secretion occurring during an early infectious process [22]. (3) The MLL rearrangement was post-natal in the ALL twin only [16].
Aside from genetic alterations, we reported B-ALL-associated DNA hypermethylation in the KMT2A/AFF1-positive twin infant at diagnosis, which reverted after chemotherapy treatment. DNA methylation changes were replicated in three independent diagnosisremission-matched samples carrying the KMT2A/AFF1 rearrangement.
GSEA performed on the diagnosis vs. remission differentially methylated genes highlighted a significant enrichment for target genes of polycomb repressive complex 2, a histone methyltransferase associated with gene repression which plays a fundamental role in the fine regulation of pluripotent stem cells [23]. Among significantly different expressed genes, crucial leukemogenic genes known to be targets of the KMT2A/AFF1 fusion protein (for example, HOXAs and MEIS1 genes) were upregulated at diagnosis, as expected [24]. We also reported DNA methylation changes that mostly translated into significantly downregulated genes related to TCEB3 expression (Supplementary Table S7). The protein product of TCEB3, Elongin A, forms a stable complex with Elongins B and C, which are required to maintain low expression levels of PCR2-repressed genomic targets [25,26]. We observed the strongest correlations between methylation changes and corresponding gene expression for the RAB32, PDK4, RANBP17, CXCL3, and MACROD2 genes. Rab small GTPases are a member of the RAS oncogene superfamily [27]. Aberrant DNA methylation of RAB32 was already described in a heterogeneous group of B-ALL pediatric patients harboring different chromosomal rearrangements, including KMT2A/AFF1. However, in comparing leukemic samples with non-leukemic samples, RAB32 gene expression and methylation were not significantly correlated, probably due to the presence of mixed B-ALL subtypes [28]. Higher RAB32 mRNA expression was also described in good-prognosis pediatric patients affected by AML compared with those with poor prognosis [29] and in adult AML patients versus healthy volunteers [30]. Therefore, lower RAB32 gene expression seems to be associated with a malignant phenotype.
PDK4 hypermethylation has been reported for pediatric B-ALL patients with the ETV6-RUNX1 translocation [31]. Bohl and colleagues found an association between PDK4 expression and the lack of a response to the hypomethylating agent decitabine in AML patients. This gene is a poor prognostic marker and has been associated with EVI1and FLT3-ITD-mediated signaling [32]. Disruption of the RanBP17/Hox11L2 region was described in ALL children [33].
The chemokine CXCL3 was found to be epigenetically repressed in B-cell ALL [28]. Consistently, this gene was downregulated in our patients at diagnosis. Chemokine signaling in ALL may have a role in the localization of leukemic blasts in specific niches, and it may confer resistance to chemotherapy [34].
In the samples at diagnosis, MACROD2 expression was considerably higher than in the remission samples. Translocations affecting the MACROD2 gene were reported in a recurrent hyperdiploid ALL pediatric patient at relapse, characterized by several chromothriptic events [35,36]. A number of studies suggest a potential role of MACROD2 in cancer; however, it is not clear whether high or low MACROD2 expression levels contribute to tumorigenesis [37]. We reported the lack of an inverse correlation between MACROD2 expression levels and the promoter methylation, which may indicate a different mechanism of gene expression regulation, such as histone modifications in the promoter or the involvement of different isoforms. In summary, these five top-ranked genes were reported here for the first time as deregulated genes in KMT2A/AFF1 infants.
Among the other genes with strong correlations between DNA methylation and gene expression (Supplementary Table S6), aberrant DNA methylation and downregulation of FKBP9, CA8, MSC, TRPC6, ZNF492, ZNF229, and ZNF471 genes in leukemia patients are reported here for the first time.
New potential therapeutic targets are necessary, particularly for infants who respond poorly to conventional chemotherapy, and allogeneic hematopoietic stem cell transplantation with an HLA-matched donor is recommended [38]. Given the reversible nature of methylation changes, they represent attractive therapeutic targets. DNA methylation profiling could be important to monitor the efficacy of these novel therapies. Moreover, the accurate assignment of patients to specific risk groups is currently a difficult and expensive process, requiring immunophenotyping, cytogenetics, and molecular diagnostics [39,40]. To conclude, the signature consisting of 32 genes whose methylation levels were reported here to be strongly and significantly correlated with expression changes could be a putative target for the development of novel epigenetic-based therapies and are worthy of further investigations to better understand the molecular mechanisms underlining the dismal prognosis of ALL infants with KMT2A/AFF1 fusions.

Cases and Biological Samples
Monochorionic twin sisters were born by cesarean section at 35 weeks of gestation, conceived after monochorionic in vitro fertilization. At 5 months of age, one child was diagnosed as pro-B ALL positive for the translocation KMT2A/AFF1, with 75% of blasts at the cytofluorimetric analysis (CD10 − , CD19 + ) and 80% of blasts at the morphological analysis. The patient was treated as a high-risk infant ALL according to the INTERFANT-06 protocol (EudraCT Number: 2005-004599-19; ClinicalTrials.gov Identifier: NCT00550992) in AIEOP (Associazione Italiana Ematologia Oncologia Pediatrica) centers in Italy. After induction chemotherapy, she achieved complete remission and then received cord blood allogenic hematopoietic stem cell transplantation from an unrelated donor.
Mononuclear cells were collected from the affected child at diagnosis (S1) and at complete remission prior to allogeneic hematopoietic cell transplantation (S3), as well as from her healthy monozygotic twin sister (S2) ( Figure 1A). PB was also collected from their healthy parents (S4, S5) ( Figure 1A). DNA samples were extracted using a QIAamp DNA Blood Kit (Qiagen, Hilden, Germany).
PB from the twin sisters was drawn during the first hours post-partum and spotted onto Guthrie cards. Neonatal DNA was extracted using a QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) from a quarter of the Guthrie cards. Details on breakpoint sequencing and searching for rearrangements in the twin sisters at birth and at diagnosis are provided in Appendix A.1.
To perform NRAS mutation analysis in the twin sisters' neonatal spots, DNA was extracted from half of the Guthrie cards using the QIAamp DNA Investigator FTA and Guthrie cards kit (Qiagen, Hilden, Germany). DNA was quantified using the Quant-iT High-Sensitivity dsDNA Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA).
In addition, we recruited four unrelated samples from infant patients with clinical characteristics in line with the twin specimens. The samples (i.e., A1 d -A1 r , A2 d -A2 r , A3 d -A3 r , and A4 d -A4 r ) were collected at the time of diagnosis ( d ) and remission ( r ); DNA (A1 d -A1 r , A2 d -A2 r , and A3 d -A3 r, ) or RNA (A1 d -A1 r , A2 d -A2 r , A3 d -A3 r , and A4 d -A4 r ) were isolated from the total bone marrow or PB (Supplementary Table S8).
Informed consent was obtained from the parents. The present study was carried out in accordance with the ethical principles of the Declaration of Helsinki.

Sanger Sequencing
NRAS exon 3 was subjected to Sanger sequencing at the Eurofins Genomics facilities (Ebersberg, Germany) to confirm the presence of the p.Q61H mutation identified by WES in the proband sample at diagnosis as well as its absence in her healthy twin sister and in their parents.

NRAS Mutation Analysis in the Twin Sisters' Neonatal Spots
The presence of the NRAS c.183A > T p.Q61H mutation was assessed using the QuantStudio 3D Digital PCR instrument and chips (Thermo Fisher Scientific, Waltham, MA, USA). Absolute quantification of the target sequence was obtained through 20,000 simultaneous PCR reactions. Details are provided in Appendix A.1. The results are reported as copies per microliters (copies/µL) calculated from the fluorescence signal generated by the amplification of the mutant allele targeted by the FAM dye-labeled probe and the wild-type allele targeted with the VIC dye-labeled probe.

DNA Methylation and mRNA Expression Analysis
Samples S1-S5 and A1 d -A1 r , A2 d -A2 r , and A3 d -A3 r underwent DNA methylation analysis using the HumanMethylation450 BeadChip (Illumina, San Diego, CA, USA). The total RNA of four B-ALL patients at diagnosis (A1 d , A2 d , A3 d , and A4 d ) and the same four patients at remission (A1 r , A2 r , A3 r , and A4 r ) were also processed using the TruSeq RNA Sample Prep V2 (Illumina, San Diego, CA, USA) to perform whole-transcriptome analysis. Additional information and statistical analysis are presented in Appendix A.1.

Statistical Analysis
In order to identify the most differentially methylated CpG sites among the proband, the remitted proband, and the healthy monozygotic sister, a threshold of at least 30% methylation difference was set (i.e., | β-value proband − β-value remission | ≥ 0.3 and | β-value proband β-value healthy sister | ≥ 0.3). Similarly, to ensure that the remaining sites were not differentially methylated between the healthy twin and the remitted proband, a threshold of at most 10% methylation difference was fixed. We named the ∆β-value the difference in the methylation level of the same CpG site between two different individuals. The resulting CpG sites were then intersected with the differentially methylated ones in the three B-ALL diagnosis-remission-matched samples, previously obtained through a paired t-test (|effect size| > 30%).
Principal component analysis (PCA), together with cluster heatmaps, was conducted considering the S1, S2, and S3 samples together with independent B-ALL diagnosisremission-matched samples (A1 d -A1 r , A2 d -A2 r , and A3 d -A3 r ) and methylation profiles from CD34 + hematopoietic stem cells isolated from fresh umbilical cord blood derived from three healthy blood donors (GEO Accession Number: GSE40799). CD19 + cells from 20 healthy adult blood donors were also used as non-leukemic reference cells (GEO Accession Number: GSE49031).
The "top" CpG sites (i.e., those identified within the twin samples and the three diagnosis-remission-matched samples) were mapped to genes also significantly differentially expressed and underwent subsequent analyses of the Pearson correlation between their methylation level and the expression variation of the corresponding gene. This possible relationship was assessed in the three B-ALL diagnosis-remission couples, because both DNA and RNA were available for them.
All the analyses were performed with the open source R (R 3.1.0) package. Gene enrichment analysis was performed using the gene set enrichment analysis computational method (GSEA; Broad Institute; http://software.broadinstitute.org/gsea (accessed on 31 December 2019)) to search for a statistically significant set of genes. The annotated gene sets are collected in the Molecular Signatures Database (MSigDB) v6.2 [41]. Enrichments with p-values < 0.05 after correction for multiple comparisons by FDR were considered statistically significant.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ijms22189740/s1. Figure S1: Breakpoint sequence for KMT2A/AFF1 rearrangement. Figure S2: Unsupervised clustering analysis of the Euclidean distances of gene expression data (30,864 genes) in the four diagnosis-remission-paired samples (couples Ai d -Ai r , i = 1, 2, 3, 4). Table S1: Sample S1. De novo single-nucleotide variants within gene regions and with a population frequency <1%. Table S2: CpG sites with at least 30% methylation differences between the twin samples (S1 and S2) and between three diagnosis-remission-matched samples (A1 d -A1 r , A2 d -A2 r , and A3 d -A3 r ). Table  S3: Gene set enrichment analysis. The top 20 overlap with the false discovery rate (FDR) (q-value < 0.05) computed for 1904 genes.  Table S6: Correlation between DNA methylation and gene expression levels. Table S7: Gene set enrichment analysis performed on 32 genes with significant correlations between DNA methylation and gene expression. Table S8: Clinical characteristics of the four unrelated leukemia patients. The presence of the NRAS c.183A > T p.Q61H mutation was assessed using the QuantStudio 3D Digital PCR instrument and chips (Thermo Fisher Scientific, Waltham, MA, USA). Absolute quantification of the target sequence was obtained through 20,000 simultaneous PCR reactions.
The purchased primers and probes for the NRAS Q61H mutation analysis were the TaqMan ® SNP Genotyping Assay (Thermo Fisher Scientific, assay ID C_166032675_10). The digital PCR mixture (15 µL) included 7.5 µL 2 × QuantStudio ® 3D Digital PCR Master Mix v2 (Thermo Fisher Scientific, Waltham, MA, USA), 0.75 µL TaqMan ® SNP Genotyping Assay, 6.75 µL DNA (10 ng), and distilled water. Using the QuantStudio 3D Digital PCR Chip Loader, the mixture was loaded on a chip of the QuantStudio 3D Digital PCR 20K Chip Kit v2. The chip was amplified using the ProFlex 2 × flat PCR system (Thermo Fisher Scientific, Waltham, MA, USA) under the following conditions: 96 • C for 10 min, followed by 40 amplification cycles (2 min at 60 • C, 30 s at 98 • C) and 2 min at 60 • C, followed by a cooling step at 10 • C. After amplification, the chip was read using QuantStudio 3D and analyzed with QuantStudio 3D Analysis Suite cloud software. The results were reported as copies per microliters (copies/µL), calculated from the fluorescence signal generated by the amplification of the mutant allele targeted by the FAM dye-labeled probe and the wild-type allele targeted with the VIC dye-labeled probe.
The quality-and quantity-checked DNA (500 ng) were treated with bisulfite (EZ-96 DNA Methylation-Gold Kit, Zymo Research Corporation, Irvine, CA, USA), and the bisulfite-converted DNA was then hybridized onto the HumanMethylation450 BeadChip (Illumina, San Diego, CA, USA) [45]. The BeadChips were scanned with a HiScan SQ instrument (Illumina, San Diego, CA, USA), and the fluorescence intensities corresponding to the methylation level of each CpG were recorded.
The raw data quality check and methylation Beta (β)-value computation were previously described in [46].
Methylation β-values with a detection p-value ≥0.01 were excluded from the analysis, as well as CpG loci with a detection p-value ≥0.01 in more than 20% of the assayed samples. All the samples had a global call rate higher than 95%.
The CpG sites were annotated to gene regions according to the Infinium Human-Methylation450K Manifest.