Epigenetic Deregulation of Telomere-Related Genes in Newly Diagnosed Multiple Myeloma Patients

Simple Summary Aberrant expression of telomere-related genes (TRGs) and telomere-shortening facilitates development of different types of cancer. Given the fact that high-risk multiple myeloma (MM) patients harbor critically shortened telomeres, we investigated the epigenetic basis of TRG-expression in newly diagnosed MM patients. We demonstrated that DNA methylation alone or in cooperation with overlapping chromatin marks and underlying genomic sequence build up an epigenetic network of aberrant gene expression in TRGs. Furthermore, we identified five TRGs, which are epigenetically controlled and consistently deregulated across the major molecular subgroups of MM, identifiable at the early stage of the disease. Abstract High-risk Multiple Myeloma (MM) patients were found to maintain telomere length (TL), below the margin of short critical length, consistent with proactive overexpression of telomerase. Previously, DNA methylation has been shown as a determinant of telomere-related gene (TRG) expression and TL to assess risk in different types of cancer. We mapped genome-wide DNA methylation in a cohort of newly diagnosed MM (NDMM; n = 53) patients of major molecular subgroups, compared to age-matched healthy donors (n = 4). Differential methylation and expression at TRG-loci were analyzed in combination with overlapping chromatin marks and underlying DNA-sequences. We observed a strong correlation (R2 ≥ 0.5) between DNA methylation and expression amongst selective TRGs, such that demethylation at the promoters of DDX1 and TERF1 were associated to their oncogenic upregulation, while demethylation at the bodies of two key tumor suppressors ZNF208 and RAP1A led to downregulation of the genes. We demonstrated that TRG expression may be controlled by DNA methylation alone or in cooperation with chromatin modifications or CCCTC-binding factor at the regulatory regions. Additionally, we showed that hypomethylated DMRs of TRGs in NDMM are stabilized with G-quadruplex forming sequences, suggesting a crucial role of these epigenetically vulnerable loci in MM pathogenesis. We have identified a panel of five TRGs, which are epigenetically deregulated in NDMM patients and may serve as early detection biomarkers or therapeutic targets in the disease.


Background
Multiple Myeloma (MM) is a B-cell neoplasm, manifested with abnormal infiltration of plasma cells (PCs) in the bone marrow matrix, which leads to the development of lytic bone lesions and myelosuppression [1]. The underlying polymorphic genetic and epigenetic variants in tumor cells strongly influence MM development and survival outcome [2][3][4][5].
In addition to etiological genetic and epigenetic features, telomere length (TL) has been identified as a crucial determinant of survival in patients with MM.
Telomeres are TTAGGG-repeat enriched nucleoprotein complexes positioned at the termini of eukaryotic chromosomes and contribute to the maintenance of genomic stability and controlled cell-proliferation [6]. In contrast, cancer cells typically maintain shortened telomeres, favored by proactive catalytic activities of defective telomerase or through an alternative lengthening of telomere mechanism [7][8][9]. Previous studies have demonstrated a strong correlation between the shortening in leukocyte telomere length and the predisposition of developing different types of cancer, including hematological malignancies [10][11][12][13]. Particularly, patients with MM harboring high telomerase activity and short TL were marked as a 'high-risk' group, having significantly (p < 0.0001; hazard ratio = 3.4) shorter overall survival [14,15]. A progressive shortening in TL was also observed in patients with MM when compared to the patients at the stage defined as premalignant monoclonal gammopathy of undetermined significance (MGUS) [16]. Thus, an increased number of telomere aggregates were ascribed as a prognostic biomarker for high-risk in the smoldering multiple myeloma (SMM) patients, consistent with their poor drug susceptibility and short-term survival [17]. While TL maintenance is found to be crucially linked to the prognosis and survival of the disease, both long-and short-term inhibition of telomerase resulted in reduced growth of MM cells [18,19].
The regulation of telomerase activity and the maintenance of TL in MM could be influenced by a set of extrinsic and intrinsic factors. For instance, cytokines such as interleukin 6 and insulin-like growth factor 1 were reported to increase telomerase activity in MM cells without alterations to the human telomerase encoding gene, telomerase reverse transcriptase (hTERT) [20]. In contrast, bortezomib, the anticancer medication used to treat MM, consistently reduced both hTERT expression and telomerase activity in MM cell lines [21]. Intrinsic factors such as genetic mutations or DNA methylation at subtelomeres and telomere-related genes (TRGs) were also reported to have a strong association with TL maintenance in different types of cancer [22][23][24][25][26]. A few previous studies have attempted to establish the relationship between differential expression of TRGs and TL maintenance in MM [27,28]. However, to our knowledge, no previous study had been reported to investigate the epigenetic basis of differential TRG expression in the disease.
To address this deficiency, we have combined epigenetic modifications, such as DNA methylation and chromatin marks, to determine their impact on the expression of a complete catalog of TRGs [29]. For the present study, we selected a cohort of patients with NDMM belonging to six major molecular subgroups, categorized based on immunoglobulin heavy chain locus (IgH) translocation (TC) with a number of oncogenes such as FGFR3/MMSET, t(4;14); CCND1, t (11;14); MAF, t(14;16); or MAFB, t (14;20); and non-TC hyperdiploid (HY) subgroups with oncogenic overexpression of multiple copies of CCND1 (D1) or CCND2 (D2). We also investigated overlaps of G-quadruplex forming sequences (G4FS) to the epigenetic risk-loci spanning the TRGs in patients with NDMM and in representative MM cell lines.

Sample Information and Genomic DNA/RNA Preparation
After approval by the institutional review board, CD138-positive bone marrow PCs were collected from a cohort of consented patients with NDMM (n = 53) and compared to pools of age-matched healthy donors (controls; n = 4). Bone marrow aspirates were enriched for >90% tumor cells using RoboSep automated cell separation techniques (Stemcell Technologies, Vancouver, British Columbia, Canada). We included patients with NDMM from the t(4;14), n = 8; t(11;14), n = 10; t(14;16), n = 9; and t(14;20), n = 7, TC subgroups and the D1, n = 12, and D2, n = 7, HY subgroups. Demographic information of the enrolled subjects is summarized in Table S1. We have also included MM cell lines, representative of the major IgH translocations found in MM patients, such as H929 for t (4;14), U266 for t (11;14), MM.1S for t (14;16), and SACHI for t (14;20) subgroups. MM Cancers 2021, 13, 6348 3 of 17 cell lines were maintained in RPMI 1640 with 20% fetal bovine serum, 1% L-glutamine, and 1% penicillin-streptomycin. Authenticity of the cell lines were checked with short tandem repeat analysis of non-coding DNA, while possibilities of microbial contaminations were nullified with mycoplasma check. Genomic DNA and RNA samples were extracted using AllPrep DNA/RNA/Protein Mini Kit (Qiagen, Hilden, Germany) or Puregene DNA Extraction Kit (Qiagen, Hilden, Germany).

Library Preparation and Sequence Annotations for Reduced Representation Bisulfite Sequencing
We prepared libraries for reduced representation bisulfite sequencing (RRBS), following our previously published protocol [30]. Briefly, 100 ng of genomic DNA from MM cell lines or patient-derived samples were digested with MspI overnight to generate~5-kb fragments along the enzyme cut sites across genome. Next, we performed the end-repair and A-tailing of DNA fragments using the KAPA HyperPlus Kit (Roche, Basel, Switzerland) followed by ligation of methylated adaptors (NEBNext Multiplex Oligos for Illumina; Methylated Adaptor, Index Primers Set; New England BioLabs Inc., Ipswich, MA, USA), bisulfite-conversion (EZ DNA Methylation Kit; Zymo Research, Irvine, CA, USA), and amplification of the barcoded libraries (EpiMark Hot Start Taq DNA Polymerase, New England BioLabs Inc., Ipswich, MA, USA). RRBS libraries were then selected based on fragment sizes of 150-250 bp and 250-450 bp. Quality control of the libraries was performed by determining the concentration, using a high-sensitivity DNA assay kit on the Qubit spectrometer (Thermo-Fisher, Waltham, MA, USA), and for size, using the 2200 TapeStation System (Agilent Technologies, Santa Clara, CA, USA), before sequencing with 75 bp singleend reads on the NextSeq 500 System (Illumina, San Diego, CA, USA). Genome sequencing of the raw data (binary alignment/map files) was processed, filtered, and the genome (hg38/GRCh38) annotated using our previously published guidelines [3]. Briefly, differentially methylated regions (DMRs) were annotated by identifying the closest transcription start site (TSS) of a gene based on the National Center for Biotechnology Information Reference Sequence Database. DMRs located from 5 kb upstream to 200 bp downstream of TSS were marked as the promoters, while DMRs spanning up to the 3 -end of the corresponding genes were marked as the body. DMRs having a methylation difference (at a false discovery rate [FDR]-adjusted p value < 0.05) of >20% increase (hypermethylation) or decrease (hypomethylation) between NDMM subgroups and controls) were considered significant for all the downstream analyses. The difference of mean DMR-methylation (MDM) at the promoters and the bodies were determined in NDMM subgroups and MM cell lines in comparison to control.

Gene Expression Analysis
Gene expression profiling was performed using Affymetrix GeneChip Human Genome U133 Plus 2.0 Array (Thermo-Fisher, Waltham, MA, USA). The raw files were processed with Transcriptome Analysis Console v4.0 Software (Thermo-Fisher, Waltham, MA, USA). The signal intensity values of multiple transcripts (whenever applicable) were averaged, while the difference in the median expression of samples between NDMM subgroups and control less than 1-fold with adjusted p < 0.05 (using log2 scale) were filtered out. Venn intersection analyses were performed using Venny v2.1 (BioinfoGP, Madrid, Spain) to determine the genes that have an inverse methylation-expression relationship at the promoters or a direct methylation-expression relationship at the bodies. A regression analysis was performed to determine the coefficient of correlation (R 2 ) between DNA methylation and expression.

Analysis of Chromatin Modifications
We have analyzed the enrichment of four activating histones (H3K27ac, H3K4me1, H3K4me3, and H3K36me3) and one repressive histone (H3K27me3), overlapped to DMRs of the differentially expressed TRGs using chromatin immunoprecipitation sequencing (ChIP-seq) in reference to H929 t(4;14) MM cells. The enrichment of CCCTC binding factor (CTCF)-transcription factor (TF) to the overlapping DMRs was also determined in reference to the ChIP-seq data in H929 cell line. Additionally, open chromatin states on or adjacent to the DMRs were determined with DNAse sequencing (DNAse-seq). Both the ChIP-seq and DNAse-seq data were obtained from the Encyclopedia of DNA Elements (ENCODE) epigenomics data (as uploaded by the Noble Laboratory at University of Washington using avocado:encode2018core) for H929 cells [31]. The chromatin enrichment and position of DMRs were visualized with the Integrated Genomics Viewer (IGV), followed by loading bigwig and customized bed files onto the IGV respectively [32].

Determine G-Quadruplex-DMR Overlaps
The G4FS, overlapped with the DMRs of TRGs, were analyzed using a web-based server that predicts quadruplex forming G-rich sequences (QGRS), the QGRS Mapper [33].

Statistical Analysis
Methylation and expression data were prepared and presented with Prism v8.0 software (GraphPad, San Diego, CA, USA). A two-way ANOVA test was used to compare pair-wise differential methylation at DMRs between NDMM subgroups and controls. A one-way ANOVA with multiple comparison was performed for expression analysis, and an adjusted p value < 0.05 was considered significant.

Selective TRGs Are Differentially Expressed in NDMM
First, we assessed differentially expressed TRGs (NDMM subgroups vs. control), which are involved in a multitude of intrinsic biological pathways including helicase, telomerase, shelterin-complex, DNA-repair, or related biosynthetic activities, as summarized elsewhere [29]. Differential expression of epigenetically regulated TRGs have also been reported in relation to their importance to the development of different types of cancer [34,35]. We aimed to evaluate the impact of 70 known TRGs in MM and observed a subgroup-specific differential expression pattern among the NDMM subgroups. [29]. For instance, we observed 65% (40/62) upregulated and 35% (22/62) Table S2). An unsupervised hierarchical clustering based on TRG-expression in NDMM subgroups revealed closest neighboring clusters between D1 and t(4;14) and between t(14;16) and t (14;20). The D2 remains equidistant to D1 and t(4;14), while t(11;14) remains equidistant to the t(14;16) and t (14;20) subgroups. In contrast, the control group formed a significantly distant cluster to all the NDMM subgroups ( Figure 1B).
Next, we narrowed down a panel of TRGs, based on their significant (p < 0.05) and consistent differential expression across the NDMM subgroups compared to control. For instance, we observed an upregulation of helicase associated DEAD-box helicase 1, DDX1  Next, we narrowed down a panel of TRGs, based on their significant (p < 0.05) and consistent differential expression across the NDMM subgroups compared to control. For instance, we observed an upregulation of helicase associated DEAD-box helicase 1, DDX1 In summary, while TRG expression was typically variable across molecular NDMM subgroups, a cluster of TRGs was consistently differentially expressed in the major molecular subgroups of NDMM discussed here.

Selective TRGs Are Differentially Methylated and Tightly Correlated to Expression in NDMM
We analyzed significant DMRs (p < 0.05 differential methylation greater or less than 20% between NDMM and controls) of 70 known TRGs across their promoters and bodies as detected per NDMM subgroups. The gene bodies were found to have multiple DMRs of a representative gene, while TRGs at the promoters were typically represented with single DMRs. In general, the majority of DMRs in both the genomic regions were found hypomethylated across the subgroups. For instance, we observed 12% hypermethylated (1/8) and 88% hypomethylated (7/8) DMRs representing six unique TRGs at the promoters and 14% hypermethylated (21/146) and 86% hypomethylated (125/146) DMRs at the bodies, representing 30 unique TRGs in NDMM subgroup t(4;14) ( Figure 3A, Table S3). In summary, while TRG expression was typically variable across molecular NDMM subgroups, a cluster of TRGs was consistently differentially expressed in the major molecular subgroups of NDMM discussed here.
Next, we focused on the TRGs that were differentially methylated in all the NDMM subgroups compared to control, and documented five TRGs at the promoters (such as DDX1, RAD54L, TP53BP1, MTR, and TERF1) and 12 TRGs at the bodies (ACYP2, DDB1, MAD1L1, MEN1, PIK3C3, PRMT8, RAP1A, RECQL4, RECQL5, RTEL1, TERT, and ZNF208). A subsequent intersection analysis containing differentially methylated and expressed TRGs among the subgroups revealed five TRGs (DDX1, PIK3C3, RAP1A, TERF1, and ZNF208) (Figure 3G), where expression is dependent on DNA-methylation. Next, we evaluated the correlation between methylation and expression of these genes and observed a strong negative correlation (R 2 ) between methylation and expression for the promoter DMRs, as observed for TRGs such as DDX1 (R 2 = 0.67) and TERF1 (R 2 = 0.49). In contrast, strong positive correlations were observed between methylation and expression for the body DMRs representing TRGs such that RAP1A (R 2 = 0.67) and ZNF208 (R 2 = 0.91) were downregulated in accordance with hypomethylation, while PIK3C3 (R 2 = 0.5) was upregulated in accordance with increase in methylation ( Figure 3H). In summary, we demonstrated a robust correlation between DNA methylation and expression of selective TRGs, evident in the early stage of the disease.

TRG-Methylation Is Fairly Conserved between Patients with NDMM and MM Cell Lines
Epigenetic polymorphisms (including DNA methylation) among patients may significantly influence variations in therapeutic response and survival outcome [36,37]. CpG sites (cytosine and guanine occurring consecutively) that are crucial to the development or maintenance of a disease are typically well-conserved between patients and cell lines [38]. To investigate the importance of the TRG-DMRs, we further compared the DNA methylation level and distribution along DMRs between TC NDMM subgroups and representative MM cell lines. We compared the MDM across NDMM subgroups and MM cell lines compared to control.  Figure 4D). For RAP1A, we observed 7 out of 14 significantly (p < 0.05) methylated DMRs in NDMM subgroups, of which DMR-1 through DMR-4 and DMR-14 were hypomethylated, and DMR-12 was hypermethylated ( Figure 4E). However, we noticed that in DMR-2, only t(4;14) and t(14;20) subgroups were significantly (p < 0.05) hypomethylated compared to control. Another interesting observation from RAP1A was the disparity in methylation between DMRs in NDMM subgroups and MM cell lines, such that DMR-5 was not differentially methylated in NDMM subgroups (MDM = −0.05, p = 0.593), but was hypomethylated in MM cell lines (MDM = −0.21, p = 0.005). Additionally, compared to hypermethylation at DMR-12 in patients, DMR-13 was particularly hypermethylated in cell lines (MDM = 0.4, p = 0.0007). Noteworthy, TRGs were also consistently upregulated or downregulated between patients and cell lines. Cumulatively, we observed considerable similarities in the DNA methylation profile at selective DMRs that were consistent between NDMM subgroups and MM cell lines with few above-mentioned exceptions.  Figure 4D). For RAP1A, we observed 7 out of 14 significantly (p < 0.05) methylated DMRs in NDMM subgroups, of which DMR-1 through DMR-4 and DMR-14 were hypomethylated, and DMR-12 was hypermethylated ( Figure 4E). However, we noticed that in DMR-2, only t(4;14) and t(14;20) subgroups were significantly (p < 0.05) hypomethylated compared to control. Another interesting observation from RAP1A was the disparity in methylation between DMRs in NDMM subgroups and MM cell lines, such that DMR-5 was not differentially methylated in NDMM subgroups (MDM = −0.05, p = 0.593), but was hypomethylated in MM cell lines (MDM = −0.21, p = 0.005). Additionally, compared to hypermethylation at DMR-12 in patients, DMR-13 was particularly hypermethylated in cell lines (MDM = 0.4, p = 0.0007). Noteworthy, TRGs

DNA Methylation Alone or Cooperatively with Chromatin Marks Regulates TRG Expression
We wanted to determine if the expression of TRGs was controlled by DNA methylation alone or in cooperation with other epigenetic factors, such as histone marks, TF binding, or chromatin state. Given the consistencies in DNA methylation across DMRs between NDMM subgroups and MM cell lines, we investigated the chromatin marks in the MM cell line H929, t (4;14) in reference to the DMRs of the NDMM patient-subgroups. We mapped proximities or overlaps between histone marks, CTCF binding sites, and open chromatin state to selective TRG-DMRs. For instance, a DMR (DMR-1 of DDX1) spanning base pairs 15,586,861 to 15,587,037 on chromosome 2 (chr2:15586861 to chr2:15587037) at 4.4 kb upstream of the TSS of DDX1, and a DMR (DMR-1 of TERF1) spanning chr8:73007342 to chr8:73007540 at 1.3 kb upstream of the TSS of TERF1 contained an abundance of activating histone H3K4me3 and relatively lesser enrichment for H3K27ac ( Figure 5A,B) MM cell line H929, t (4;14) in reference to the DMRs of the NDMM patient-subgroups. We mapped proximities or overlaps between histone marks, CTCF binding sites, and open chromatin state to selective TRG-DMRs. For instance, a DMR (DMR-1 of DDX1) spanning base pairs 15,586,861 to 15,587,037 on chromosome 2 (chr2:15586861 to chr2:15587037) at 4.4 kb upstream of the TSS of DDX1, and a DMR (DMR-1 of TERF1) spanning chr8:73007342 to chr8:73007540 at 1.3 kb upstream of the TSS of TERF1 contained an abundance of activating histone H3K4me3 and relatively lesser enrichment for H3K27ac (Figure 5A,B). The chromatin adjacent to DMRs of both the genes was found in an open state, as evidenced from DNAse peaks. Additionally, a sharp CTCF peak was observed at the upstream TSS of DDX1. Hypomethylation at the promoters, coincident with activating histone marks and open chromatin state, explained upregulation of DDX1 (DE: 1.2, log2) and TERF1 (DE: 1.2, log2).

Figure 5.
Proximity and possible overlaps between DMRs and chromatin marks were mapped on selective epigenetically regulated TRGs (n = 5). ChIP-sequencing data from H929, t(4;14) cell line containing four activating histones (H3K4me1, H3K4me3, H3k27ac, and H3K36me3), and one repressive histone (H3K27me3) were mapped to the proximal DMRs. Additionally, DNAse-seq peaks were observed to determine open chromatin, while CTCF binding was also observed. The overlaps between DNA and chromatin level epigenetic marks were visualized at 5 kb upstream TSS of (A) DDX1 and (B) Figure 5. Proximity and possible overlaps between DMRs and chromatin marks were mapped on selective epigenetically regulated TRGs (n = 5). ChIP-sequencing data from H929, t(4;14) cell line containing four activating histones (H3K4me1, H3K4me3, H3k27ac, and H3K36me3), and one repressive histone (H3K27me3) were mapped to the proximal DMRs. Additionally, DNAse-seq peaks were observed to determine open chromatin, while CTCF binding was also observed. The overlaps between DNA and chromatin level epigenetic marks were visualized at 5 kb upstream TSS of (A) DDX1 and (B) TERF1, or across the body region of (C) ZNF208, (D) PIK3C3, and (E) RAP1A. Abbreviations: ChIP, chromatin immunoprecipitation; CTCF, CCCTC binding factor; DMR, differentially methylated regions; DNAse-seq, DNAse sequencing; TRG, telomere-related gene; TSS, transcription start site.
In contrast, hypomethylation observed at DMR-2 spanning chr19:21940986 to chr19:21941044, DMR-5 along the body of ZNF208 spanning chr1:111608116 to chr1:111608192, and DMR-13 spanning chr1:111689837 to chr1:111689970 across RAP1A gene-body were consistent with downregulation of ZNF208 (DE: −0.6, log2); and RAP1A (DE: −0.7, log2) ( Figure 5C,E). While no prominent chromatin marks were found across the promoter or the body of ZNF208, subtle histone enrichments and open chromatin loci were observed in the body of RAP1A; however, were considerably distant (>5 kb) to the DMRs. In these two examples, DNA methylation is likely to control TRG expression, independently or with chromatin modifications. The DMR-2 spanning chr18:42024964 to chr18:42025106) of PIK3C3 was hypermethylated and proximal (<5 kb) to CTCF binding sites, which is consistent with upregulation (DE: 0.61, log2) of the gene ( Figure 5D). Cumulative data demonstrated variable epigenetic regulations, linked to the TRG expression in the disease.

TRG Hypomethylated DMRs Are Stabilized with G-Quadruplexes
Given the fact that aberrantly hypomethylated genomic regions, enriched for G4FS have significant contributions to the genomic instabilities in cancer [39], we investigated overlaps of G4FS with the DMRs of TRGs to understand the possible impact of these geneticepigenetic variants on sustaining stability of these TRG-related epigenetic risk-loci in the disease. We observed three G4FS overlapped with the hypomethylated DMR-1 (177 bp) at the promoter of DDX1 (Figure 6A), which also coincides to the open chromatin region.

Discussion
To the best of our knowledge, this is the first study to report that DNA methylation alone, or in cooperation with chromatin state and histone marks, deregulate aberrant selective TRG expression in early-stage MM.
In general, the majority of TRGs in our NDMM cohort were differentially upregulated (81%) rather than downregulated (19%) across the subgroups. For instance, upregulation of the TRGs, such as DKC1, GAR1, NHP2, and NOP10, were found in several NDMM subgroups (Table S2), which is consistent to the previous finding where expression of these genes was progressively increased along with the development from MGUS to the MM stage [27]. In contrast, WRAP53, which is known to bind the Cajal body box of the telomerase RNA template component, was downregulated (DE: −0.37, log2). Reduced expression of WRAP53 has been previously correlated to attenuated DNA damage response and poor survival in ovarian and lung cancers [40,41], and suggests similar line of dysregulation in MM. Of the shelterin proteins, TERF1 was upregulated, while RAP1A was downregulated, which suggest defunct telomere maintenance, favorable of progression of the disease. In particular, TERF1 overexpression was found strongly correlated with TL and disease progression in gastric and renal carcinomas [42,43]. Another interesting finding from the TRG expression profile was the upregulation of 3 (DCLRE1C, SMC5, and SMC6) of 15 detected repair related TRGs in NDMM subgroups compared to control. DNA-repair genes have previously been found to be upregulated in response to chemotherapy and radiotherapy in a number of cancers including melanoma manifested with metastasis [44,45]. In summary, the TRG expression profile suggests an involvement of a compromised DNA-repair mechanism machinery, coincident with a dysfunctional telomere maintenance in the disease.
Next, we found a strong correlation (R 2 ≥ 0.5) between DNA methylation and expression of five TRGs. In comparison to the previous studies, where DNA methylation changes at individual CpGs were correlated to changes in TRGs expression, we considered all significant DMRs along regulatory regions of the TRGs to determine the methylationexpression relationship in MM [29]. Moreover, we have investigated the role of differential methylation at both the promoters and the bodies of TRGs, considering a unique crucial role of the regulatory regions on gene expression [46,47]. We observed the promoter of DDX1 and TERF1 was hypomethylated, where DMRs of both genes were proximal to H3K4me3 and H3K27ac enriched TSS. This corroborates the fact that hypomethylation at the promoter functions cooperatively with activating chromatin marks to facilitate binding of TFs to upregulate these oncogenes. This seems more reliable given the appearance of DNAse peaks suggesting open chromatin state on the genes. Additionally, we observed CTCF binding to the DDX1 promoter. This could be explained with colocalization of CTCF with the promoter bound H3K4me3 or H3K27ac regarding upregulation of the gene [48]. We also observed a strong positive correlation between DNA methylation and expression at the bodies, such that RAP1A and ZNF208 were downregulated, consistent with hypomethylation across the body DMRs, while PIK3C3 was upregulated, consistent with hypermethylation across the body DMRs. This finding is highly fascinating because aberrant upregulation of several oncogenic TRGs is complemented with epigenetic silencing of two key tumor suppressors, RAP1A and ZNF208 [49][50][51]. Additionally, gene-body DNA methylation was found as a principal determinant of expression because DMRs were not overlapped or proximal (>5 kb) to any regulatory histone marks or CTCF binding sites at open chromatin regions.
Our study provides a potential mechanistic insight of robust epigenetic control of TRG expression in patients with NDMM, which could be considered as potential earlydetection biomarkers and therapeutic targets for the disease. The DMRs, as observed in the NDMM subgroups might also have a profound role in the progression of the disease, as we found several of these DMRs to be conserved between NDMM subgroups and MM cell lines, which represent a terminal differentiation stage of the disease. However, we found a relatively higher degree of hypermethylation and hypomethylation in the cell lines compared to the patient-derived samples. We also found a general trend of hypomethylation (except for PIK3C3) in TRG-DMRs, which were frequently overlapped with G4FS. Our findings are in alignment with previous findings that G4FS may stabilize and conserve CpG sites at hypomethylated states by inhibiting DNA methyltransferase 1 enzymatic activity, which in turn facilitates enhanced DNA breakpoints or chromatin accessibility [39,52,53].
In summary, we have identified a number of crucial epigenetic risk-loci that might play crucial roles in TRG dysfunctions and overall development of the disease. These epigenetic risk-loci are consistently deregulated in NDMM subgroups and detectable at the early stage of the disease. However, future study needs further expansion with inclusion of subjects from different disease stages (such as premalignant MGUS, transient SMM, and MM with patients from relapsed and refractory cohorts) to compare dynamics of epigenetic alterations at TRG loci. Finally, given the advancement in sequencing techniques to accurately match the DNA methylation profile of cell-free DNA to the tissues and cell types, it could serve as a noninvasive screening platform for epigenetic biomarkers in patients with MM [54].

Conclusions
We have identified a 5 TRG panels, tightly regulated with epigenetic modifications such as DNA methylation, chromatin modifications, or G4FS in patients with NDMM. Using clustered regularly interspaced short palindromic repeats (CRISPR) tools, it would be highly interesting to investigate whether conditional alterations in DNA methylation or chromatin marks on these TRGs could alter gene expression, TL and inhibit overall progression of the disease. Investigations of epigenetically deregulated TRG pathways are also expected to provide valuable information for developing novel chemotherapies and immunotherapies in MM.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/cancers13246348/s1, Table S1: Demographic features of patients with NDMM, Table S2: Median expression of TRG transcripts in major molecular subgroups of patients with NDMM, Table S3: Differential methylation across the promoter and the body of TRGs in t(4;14) subgroup, Table S4: Differential methylation across promoter and body of TRGs in t(11;14) subgroup, Table  S5: Differential methylation across promoter and body of TRGs in t(14;16) subgroup, Table S6: Differential methylation across promoter and body of TRGs in t(14;20) subgroup, Table S7: Differential methylation across promoter and body of TRGs in D1 subgroup, Table S8: Differential methylation across promoter and body of TRGs in D2 subgroup.
Author Contributions: S.R.C. conceptualized the work and wrote the manuscript. S.R.C. performed and analyzed DNA methylation. The omics data was annotated and further analyzed by C.A. and S.R.C., F.Z. and F.v.R. helped in identifying, characterizing, and selecting NDMM patient samples. All co-authors contributed to the final version of the manuscript. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The datasets used and/or analyzed during the current study are available from the corresponding author and with permission of University of Arkansas for Medical Sciences on reasonable request.

Acknowledgments:
The authors acknowledge the participants who generously gave their time to the study. UAMS Myeloma Center foundation grant is acknowledged for the generous fund support. Authors are thankful to Yan Xiao and Geethamary Gerad in the accessioning lab for their instrumental support in the acquisition and processing the patients' samples, Ruslana Tytarenko and Purvi Patel, who processed part of the methylation and gene expression array, and Yan Wang, who helped CA in annotating RRBS data for the present study.

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