Abstract

   Multiple myeloma (MM) is a hematological malignancy characterized by
   expanding clonal plasma cells in the bone marrow (BM) that produce
   monoclonal immunoglobulin. It is an incurable disease, accounting for
   about 10% of blood malignancies and the second most common hematologic
   malignancy. Therefore, in-depth research into the molecular mechanisms
   and therapeutic targets of the disease is crucial. For the first time,
   we performed high-throughput chromosome conformation capture (Hi–C)
   analysis of plasma cells in five multiple myeloma patients, and
   integrated it with genome resequencing and transcriptomic associated
   with genomic variation and gene expression. As a result, 19 specific
   TAD (Topologically Associating Domain) boundaries in MM samples related
   to the immune response and Wnt signaling pathways were identified.
   Additionally, Loop structures were also analyzed, revealing that
   promoter-enhancer-associated loops were the most prevalent. Genomic
   characteristics of MM patients were explored, identifying SNPs, InDels,
   and CNVs, with variations in the CDS region potentially affecting gene
   function. Transcriptome analysis showed differentially expressed genes
   in MM patients, mainly involved in p53 signaling and cell adhesion.
   Multi-omics analysis identified overlapping genes related to MM,
   including those involved in MHC class II protein complex assembly and
   antigen presentation. The study provides insights into the complex
   genomic and transcriptomic changes in MM plasma cells, potentially
   aiding in identifying therapeutic targets.

Supplementary Information

   The online version contains supplementary material available at
   10.1038/s41598-025-03132-2.

   Keywords: 3D genome, Genome resequencing, Transcriptome, Plasma cell,
   Multiple myeloma

   Subject terms: Cancer, Molecular biology, Medical research, Molecular
   medicine

Introduction

   Multiple myeloma (MM) is a form of cancer that predominantly affects
   adults and is centered around the excessive multiplication of plasma
   cells within the bone marrow^[46]1. It is an incurable disease, and
   most patients experience relapse/refractory^[47]2, posing significant
   challenges for patients and healthcare providers. However, targeted
   immunotherapeutic strategies such as CAR-T, CAR-natural killer (NK),
   and TCR-T have emerged as a beacon of hope in the treatment of MM due
   to their potential to minimize harm to healthy tissues^[48]3,[49]4.

   Recent applications of chromosome conformation capture techniques,
   particularly Hi–C, have revolutionized our understanding of the 3D
   cancer genome^[50]5. Studies using Hi–C data from different species
   have shown that the switching of chromatin compartments between active
   (A) and inactive (B) states is closely linked to changes in gene
   expression^[51]6–[52]8. In prostate cancer, topologically associating
   domains (TADs) are reduced in size and exhibit alterations at key tumor
   suppressor loci, such as TP53^[53]9. Previous studies have shown that
   the destruction of TAD boundaries will lead to the regulation of the
   genes and the occurrence of diseases^[54]10,[55]11. Therefore, the TAD
   boundary plays a very important role. Deleting the boundary section
   will cause a disorder of gene regulation, resulting in the
   transcription of originally silenced genes, while those that should
   have been transcribed were silenced^[56]12. The boundaries of TAD can
   not only guide the folding of chromosomes into higher structures, but
   also correctly guide the remote transcriptional regulation^[57]13.
   Moreover, the spatial architecture of the genome actively shapes the
   nature of cancer-associated genomic alterations^[58]14. These insights
   reveal a bidirectional relationship between the 3D organization of the
   genome and the genetic changes that drive cancer^[59]15. Understanding
   this interplay is crucial for unraveling the molecular mechanisms of
   cancer and for developing targeted therapies that consider the spatial
   context of genomic alterations^[60]16. A previous study found that in
   the context of multiple myeloma (MM), the three-dimensional (3D)
   organization of the genome undergoes significant changes when compared
   to normal B cells^[61]5. These alterations are characterized by an
   increase in the number of topologically associating domains (TADs), a
   reduction in the average size of these domains, and a shift in the
   chromatin state of certain genomic regions^[62]5,[63]11. Multi-omics
   analysis is an important method to explore the formation of traits and
   the development of diseases^[64]5,[65]17. However, limited multiple
   myeloma study was explored using 3D genomics methods based on the
   different patients.

   Despite considerable progress, the development of extensive 3D datasets
   that capture the cancer genome’s structure is still in its infancy.
   Therefore, in-depth research into the disease’s molecular mechanisms
   and therapeutic targets is crucial. Plasma cells are a specialized type
   of white blood cell that plays a vital role in the immune system. They
   are derived from B lymphocytes, a type of immune cell that produces
   antibodies. 3D genomics is an emerging field of study that focuses on
   the spatial organization of the genome. However, the 3D genomic study
   of plasma cells is limited. Characterizing the cancer genome’s spatial
   disarray and its functional implications is paramount, especially
   considering the prevalence of genomic changes such as mutations in
   cancer. In our study aimed at unraveling the molecular characteristics
   of MM cancers, we employed a multi-faceted analysis. We utilized Hi–C
   contact maps, WGS, and RNA-seq data from 5 patients and 1 control
   normal person (Table [66]S1) to scrutinize potential biases in Hi–C
   data. Our findings unveiled a significant link between the
   architectural features of the 3D genome, genomic variations, gene
   expression, and MM cancer. This integrated analysis enhances our
   comprehension of the 3D cancer genome’s role in MM and opens new
   avenues for therapeutic intervention by targeting genomic disarray.

Results

Clinical information about the patients

   The cohort of this study consists of individuals aged between 44 and
   80, with both male and female participants (Table [67]S1). The group
   includes five MM patients, each with distinct treatment regimens. The
   treatments range from combinations of bortezomib, thalidomide, and
   dexamethasone to more complex therapies involving pomalidomide,
   bortezomib, dexamethasone, daratumumab, and selinexor. The time of
   diagnosis varies, with the earliest in August 2018 and the latest in
   July 2023. Genetically, the patients exhibit a range of chromosomal
   abnormalities. For instance, several individuals have FISH results
   showing 1q21 amplifications and deletions of 13/13q. One patient has a
   complex karyotype with multiple abnormalities, including deletions,
   translocations, and aneuploidy, while another shows a highly abnormal
   karyotype with numerous structural and numerical changes. In contrast,
   one patient has normal FISH and karyotype results. Overall, this group
   is characterized by significant genetic heterogeneity and diverse
   treatment approaches, reflecting the complexity of their conditions. In
   addition, all MM patients were seriously ill and died within 2 years
   after diagnosis, these samples had certain characteristics and were
   difficult to collect. With this clinical context in mind, we next
   explored the 3D genomic architecture of MM plasma cells to understand
   how chromatin organization might contribute to disease progression.

Dynamic changes in compartmentalization and local accessibility

   To elucidate the multiscale rewiring of chromatin architecture and its
   influence on MM plasma cells, we used in situ Hi–C to map chromatin
   contacts for plasma cells of five patients (MMC1, MMC2, MMC3, MMC5, and
   MMC6) and one control sample (Control). As a result, the average clean
   data of the sample is 180.37 Gb. We then generated a total of 1.45
   billion valid contacts (with an average of 242.28 million contacts per
   sample (Table [68]S2) and reached a maximum resolution of 8.15 kb
   (Table [69]S2). Most (57.18%) contacts occurred within chromosomes and
   consisted of the dominant (87.61%) long-range interactions (Table
   [70]S2).

   We then constructed genome-wide contact maps by dividing the genome
   into 500 kb regions (Figs. [71]1A and [72]S1A). Inter‑chromosomal
   interactome indicated that the chromosomes with similar lengths have a
   similar likelihood to mutually contact each other, which was consistent
   with the previous study^[73]5. All samples showed a strong decrease in
   contact probability with an increase in the distance between loci
   (Fig. [74]1B). Compartment A, rich in actively transcribed genes,
   features open chromatin with histone marks for active transcription and
   is centrally located within the nucleus. In contrast, Compartment B,
   with fewer genes and inactive transcription, has closed chromatin with
   marks for gene silencing and is peripherally situated in the
   nucleus^[75]8. These compartments are crucial for genome organization
   and gene regulation^[76]18. We identified the compartment state of each
   sample (Fig. [77]1C) and identified substantial levels of compartmental
   switching in plasma cells across control and 5 MM patients
   (Fig. [78]1D). In these regions, 408 bins were switching from B to A
   (which contained 259 genes) and 321 bins from A to B (contained 274
   genes) (Fig. [79]1E).

Fig. 1.

   [80]Fig. 1
   [81]Open in a new tab

   Overview of 3D organization in MM plasma cells. (A) Genome-wide contact
   maps of MMC1 patients. (B) P(s) curves (at 100 kb resolution) averaged
   across all autosomes in the genome of each sample. (C) Compartment
   state of chromosome 1 indicated by PC1 values of MMC1 patients. Red
   represents A compartment and cyan represents B compartment. (D)
   Compared with Control, the MM patients showed a few compartment changes
   in common. (E) The most enriched terms for genes with A to B or B to A
   switching. F. Enrichment analysis of genes with A-to-B switching event.

   Functional enrichment analysis demonstrated that genes embedded in
   regions experience the A-to-B switching event and were primarily
   involved in signaling receptor activity and G protein-coupled receptor
   activity (Fig. [82]1F). This includes molecular transducer activity,
   transmembrane signaling receptor activity, interleukin-1receptor
   activity and detection of chemical stimulus involved in sensory
   perception (Fig. [83]1F). Nonetheless, genes located in regions that
   were subject to B-to-A switching events were primarily involved in
   DNA-binding transcription factor activity and RNA polymerase II
   transcription regulatory region sequence-specific DNA binding molecular
   function, it is related to the development of organs and system of
   animals. In addition, it is also enriched in type III interferon
   signaling (Fig. [84]S1B).

Most TADs were highly stable

   Having established compartment-level changes in chromatin organization,
   we next examined the stability and alterations of topologically
   associating domains (TADs) in MM plasma cells to identify specific
   regulatory units associated with disease development. At the
   sub-megabase scale, the local chromatin architecture can be
   characterized by TAD^[85]19. A topologically associating domain (TAD)
   is a highly self-correlated continuous region where the interactions
   between fragments tend to be more within the TAD than between TADs, and
   it is separated from its neighbors by distinct boundaries to form an
   independent regulatory unit that presents a square structure on the
   heat map diagonal. As a regulatory unit, genes in TAD share common
   regulatory elements, so there are cooperative expression
   characteristics of genes in TAD (providing basis for co-expression of
   adjacent genes on chromosomes). According to the Inclusion ratio (IR),
   we identified 3197–3988 TADs in the six samples (Fig. [86]2A), with an
   average of 702.45 Kb length (Fig. [87]2B, Table [88]S2).

Fig. 2.

   [89]Fig. 2
   [90]Open in a new tab

   TAD and loop characteristics. (A) Number of TADs in each sample. (B)
   Boxplots showing TAD sizes in each sample. (C,D) GO enrichment
   biological process (BP) (C) and KEGG (D) analysis of genes with TAD
   boundary of difference. (E) Hi–C contact matrix and TADs on Chr 1:
   1–10 Mb in control and MMC1 patient.

   We next calculated the genome-wide DI (Directionality index score)
   values of the samples to explore the differential TADs among all
   samples. Then, the TAD boundaries of all samples are merged. As a
   result, compared with the control, we found 19 specific TAD boundaries
   (which embedded 43 genes) in the MM samples. These genes were mainly
   involved in “immune response”, “neutrophil degranulation”, “leukocyte
   degranulation”, “leukocyte activation” and “immune system process”
   biological process terms (Fig. [91]2C); “Human immunodeficiency virus 1
   infection”, “Wnt signaling pathway”, and “Cytokine − cytokine receptor
   interaction” KEGG pathways (Fig. [92]2D); “molecular function
   regulator”, “catalytic activity, acting on a protein”, and
   “guanyl − nucleotide exchange factor activity” molecular function terms
   (Fig. [93]S2A). Figure [94]2E showed an example fragment of 10 Mb on
   chromosome 1 with resolution of 40 Kb. These results indicate altered
   immune response due to formation of new and disappearance of the
   original TAD boundaries associated with multiple myeloma development.
   Building on the TAD analysis, we next explored chromatin loop
   structures to understand how promoter-enhancer interactions might
   contribute to gene expression changes in MM.

Global rewiring of loops in multiple myeloma

   If the interaction frequency of a pair of chromosomal sites is higher
   than the interaction frequency of the adjacent chromosome segments on
   the linear line, then the pair is called a significant interaction site
   (which is also called loop)^[95]20, that shows the presence of strong
   interaction signal points at non-diagonal locations on the heat map.
   The existence of loop structure is the biological reason for the
   emergence of significant interaction sites, so we can identify
   significant interaction sites and identify loop structure through Hi–C
   data. The ends of the loop are referred to as the anchor points of the
   loop, which include common promoter-enhancer interaction (PEI) sites.
   We next identified significant interacting sites through Hi–C data and
   identifying loop structures. As a result, we identified 1069 loops from
   MMC6 to 6929 loops in control (Table [96]S2). If a loop where one
   anchor is in the promoter region (the 2 kb upstream of the gene
   transcription start site serves as the promoter region), and the other
   anchor is located in a non-promoter region (potential enhancer-like
   region) is referred to as a promoter-enhancer associated loop (PEL).
   Finally, we identified the greatest number of loop types was
   promoter-enhancer associated loops. Most loops were not anchored
   enhancers or promoters (Fig. [97]3A–C), and the number of each loop
   type shows a significant difference; given that patients MMC1 and MMC2
   (normal karyotypes) have different karyotypes compared to patients
   MCC4, MMC5, MMC6 (complex karyotypes, Table [98]S1), it is presumed
   that MM patients with complex genomic variants lose more
   regio-interactions in the genome.

Fig. 3.

   [99]Fig. 3
   [100]Open in a new tab

   Loop characteristics in MM. (A) The identified Loop is classified into
   three types. E-E, E-P, and P-P are the number of enhancer-enhancer
   interaction, enhancer-promoter interaction, and promoter-promoter
   interaction loops, respectively. Loop classification according to the
   number of enhancers (B) and promoters (C) of each loop contained. (D)
   Distribution of chromatin loop in each chromosome in control (red) and
   MM (green) groups. GO enrichment biological process (BP) (E) and KEGG
   (F) analysis of genes with loops of difference.

   The differential loop analysis showed 10 specific loops in the 5 MM
   samples compared with the control (Fig. [101]3D). These differential
   loop-related genes (if a bin size region upstream and downstream of the
   differentially loop-related boundary intersects with the promoter
   region of the gene) were further used to perform GO/KEGG enrichment
   analysis. These genes were mainly involved in “anterior/posterior
   pattern specification” biological process terms (Fig. [102]3E);
   “Epstein − Barr virus infection”, “Autophagy – animal”, “Relaxin
   signaling pathway”, and “ECM − receptor interaction” KEGG pathways
   (Fig. [103]3F); “regulatory region nucleic acid binding”, and
   “chromatin DNA binding” molecular function terms (Fig. [104]S2B). These
   findings indicate the potential role of chromatin loops in autophagy
   events in MM patients.

Genomic characteristics of the MM

   With insights from chromatin loops, we then turned to the genomic
   characteristics of MM to identify specific mutations and structural
   variations that might drive disease progression. To explore the genomic
   characteristics of MM patients, we determined the whole genomes of 5
   patients, compared them with one control person, and identified the
   SNPs and other variations (Fig. [105]4A). As a result, we totally
   identified 6.58 M, 4.40 M, and 2.18 M SNPs, transitions, and
   transversions, respectively (Table [106]S3). In total, there’re
   3,636,606 SNPs were identified in the intergenic region (55.28%), and
   only 45,438 SNPs (0.69%) were identified in CDS regions (Table
   [107]S4). Of these SNPs in CDS region, most of them were nonsynonymous
   (49.62%) and synonymous coding SNPs (49.23%) (Fig. [108]4B). Similar to
   SNPs, small InDel annotation results also showed that most InDels were
   in intergenic regions (51.17%) (Table [109]S5). We also detected the
   structure variation (SV) of the MM samples. The greatest number of SV
   (69,169) was found in MMC2 (Table [110]S6), most of them were inversion
   (65,308). Copy number variation (CNV) detection results showed that the
   control has the least number of CNV (1067). These variations in the CDS
   region might cause the function change of the genes, compared with the
   reference genome, we found there were 6752, 927, and 1211 genes with
   non-synonymous SNP, InDel, and SV, respectively, for MMC1 patients
   (Table [111]S7). There are 7765, 12,497, 7957, 8266, and 8034 genes
   with variations (including Non-synonymous SNP, InDel and SV) in MMC1,
   MMC2, MMC4, MMC5, and MMC6 respectively. Among these genes, 4691 of
   them were shared among these five MM patients (Fig. [112]4C).

Fig. 4.

   [113]Fig. 4
   [114]Open in a new tab

   Genomic variations of MM. (A) Distribution of variations on the genome.
   From outside to inside were chromosome position, gene density, SNP
   density and InDel density respectively. (B) Distribution of SNPs of
   each type based on all samples. (C) Veen diagram showing number of
   genes with variations (including Non-synonymous SNP, InDel and SV)
   shared among the five MM patients. KEGG (D) and Reactome (E) enrichment
   analysis of the shared genes with variations in the CDS region.

   These overlapped variated genes were further used to perform GO/KEGG
   enrichment analysis. These genes were mainly involved in “multicellular
   organismal process”, “response to stimulus”, and “developmental
   process” biological process terms (Fig. [115]S3A); “ECM proteoglycans”,
   “Collagen formation”, “Collagen chain trimerization”, and “Diseases of
   glycosylation” Reactome terms (Fig. [116]4D); “ECM-receptor
   interaction”, “Olfactory transduction”, “Graft-versus-host disease”,
   and “Taste transduction” KEGG pathways (Fig. [117]4E); “small molecule
   binding”, “ion binding”, and “protein binding” molecular function terms
   (Fig. [118]S3B). These results implicate SNPs, InDels and SVs in the
   development of MM were associated with binding processes.

Transcriptome changes in MM plasma cells

   To further link genomic variations to gene expression changes, we
   performed transcriptome analysis to identify differentially expressed
   genes and pathways in MM plasma cells. We sequenced five samples
   (including Control, MMC1, MMC2, MMC4 and MMC6). RNA-seq results showed
   that after sequencing quality control, a total of 34.01 Gb of clean
   data was obtained (Table [119]S8). According to the expression of the
   samples (Fig. [120]5A), the correlation analysis showed that these MM
   patients showed more similar expression pattern compared with control
   (Fig. [121]5B). Compared with control, we identified 1,619 (1,209 up
   and 410 down regulated genes) differentially expressed genes (DEGs) in
   the four MM patients (Fig. [122]5C,D). These DEGs were mainly involved
   in “p53 signaling pathway”, “Various types of N-glycan biosynthesis”,
   “Pathways in cancer”, “Cell adhesion molecules” and “N-Glycan
   biosynthesis” KEGG pathways (Fig. [123]5E); “protein disulfide
   isomerase activity”, “ADP binding” and “AMP binding” GO molecular
   function terms (Fig. [124]S4A); “negative regulation of cell adhesion”,
   “negative regulation of cell communication”, and “negative regulation
   of signaling” GO biological process terms (Fig. [125]S4B). These
   results indicate a large proportion of perturbations in the
   transcriptome of MM blood with the majority of dysregulated
   protein-coding genes associated with adenylate binding.

Fig. 5.

   [126]Fig. 5
   [127]Open in a new tab

   Analysis on DEGs. (A) FPKM density distribution of each sample. (B)
   Correlation heatmap between samples. (C) MA plot of differentially
   expressed genes between MM patients and Control. (D) Hierarchical
   clustering of differentially expressed genes. (E) Bubble chart of KEGG
   pathway enrichment on DEGs. Top 20 enriched pathways (with smallest
   Q-value) were shown.

   Tumors are intricate conglomerates of both malignant and non-malignant
   cells. The tumor purity—defined as the proportion of cancer cells
   within a sample—can complicate integrative analyses by introducing
   variability^[128]21. Conversely, it also presents an opportunity to
   study tumor heterogeneity, offering insights into the complex interplay
   between cancer cells and their surrounding microenvironment^[129]22. We
   then used PUREE^[130]23 to calculate the tumor purity of these MM
   samples. As a result, the MM samples had tumor purity values ranging
   from 0.46 (MMC1) to 0.58 (MMC4), which belonged to mid-high (0.38–0.97)
   purity range samples^[131]23. The mixture of the genomes of tumor cells
   and normal cells, which may be responsible for the large number of
   mutations in the genome.

Multi-omics analysis

   Integrating the multi-omics data, we explored the overlap between
   genomic variations, chromatin organization, and gene expression to
   identify key genes and pathways that may drive MM development. We first
   checked the overlapped genes among differentiated compartments, TAD,
   and loop (Fig. [132]6A). Three hundred and fifty overlapped genes were
   found in DEGs and genomic CDS variated genes (Fig. [133]6B). These
   genes were mainly involved in “MHC class II protein complex assembly”
   and “Antigen processing and presentation of exogenous peptide antigen”
   GO biological process terms (Fig. [134]6C); “fibronectin binding”, and
   “MHC class II protein complex binding” GO molecular function terms; and
   “Autoimmune thyroid disease”, “Cell adhesion molecules”, “N-Glycan
   biosynthesis”, “Intestinal immune network for IgA production” and
   “Hematopoietic cell lineage” KEGG pathways (Fig. [135]6C). In addition,
   we found that genes such as NSG2, ALDH1A2, and CALCRL showed a higher
   expression in the MM patients. ARHGAP24 showed a higher expression in
   the control (Fig. [136]6D). NSG2 is a member of the neuron-specific
   gene (NSG) family, which is specifically expressed in neurons^[137]24.

Fig. 6.

   [138]Fig. 6
   [139]Open in a new tab

   Multi-omics analysis. (A) Venn diagram showing the overlapped genes for
   differentiated compartment change, TADs and loops. (B) Venn diagram
   showing overlap between genes with CDS variation and DEGs. (C)
   Enrichment analysis on overlapping genes between genes with CDS
   variation and DEGs. (D) FPKM of ARHGAP24, ALDH1A2, CALCRL, and NSG2
   genes in four MM patients and Control.

   Cytogenetic analysis has unveiled a complex landscape of genetic
   mutations in multiple myeloma (MM), with the majority of these
   alterations being concentrated in structural rearrangements and copy
   number variations (CNVs)^[140]25. Among these, the most commonly
   observed CNVs include gains of chromosome 1 and loss of chromosome
   17^[141]26,[142]27. Here, consistent with previous studies, we also
   found that the most commonly observed CNVs include gains of chromosome
   1 and losses of chromosome 13 were detected (Table [143]S9). We also
   found that CNV had an important effect on the loops difference in MM
   patients. Among the 161 loops differentiated between control and MM
   patients, 141 of them were related to genes (39% of these loops were
   related to CNV genes). These genes include LRRC63, GPR183, OBI1,
   POU4F1, ZIC5, COG3, UBAC2, TM9SF2, GPR18, and ZIC2. For most of these
   genes, decreasing number of interactions were detected, the
   interactions were lower in MM than the control, with an average − 5.73
   of log2FC.

   These results indicate the functional relevance of each level of the 3D
   genome hierarchy (such as compartment, TAD and loop), and the
   variations of genome and establishment of 3D genome provides multiple
   regulatory layers to alter gene expression in MM development.

Discussion

   Chromatin conformation capture techniques, such as 3C, 4C, 5C, Hi–C,
   and ChIA-PET, have recently been developed to explore the
   three-dimensional (3D) genome organization of genomes at
   high-resolution^[144]7,[145]19 and reveal gene regulation
   mechanisms^[146]11. Using these techniques, many studies have found
   that the mammalian and bird’ genomes are organized into gene-dense and
   transcriptionally active compartment A and gene-sparse and
   transcriptionally inactive compartment B at the megabase
   scale^[147]6,[148]11. Topologically associating domains (TADs) are
   formed at the sub-megabase scale, which functions as units for gene
   regulation^[149]28. Chromatin loops facilitate long-range interactions
   between enhancers and promoters for gene regulation within
   TADs^[150]20. The 3D organization of the genome is dynamically
   regulated in key biological processes such as stem cell
   differentiation^[151]29, cell division^[152]30, and B-cell
   activation^[153]31.

   This study comprehensively analyzes Hi–C, genomic resequencing, and
   RNA-seq across five multiple myeloma (MM) and control plasma cells.
   Consistent with previous study^[154]5, we found that the 3D cancer
   genome is influenced by cancer-specific genome alterations and
   differential gene expression events. For the majority of these CNV
   related genes, a reduction in the number of detected interactions was
   observed. The level of interactions was found to be lower in multiple
   myeloma (MM) compared to the control groups. Many genes related to CNV
   were also found with loop difference, such as GPR18 and GPR183. These
   two genes both showed CNV loss in MM patients. The GPR18 gene is
   related to the number of B cells and the expression of B cell-related
   genes^[155]32, which may have clinical value for the prognosis of
   various cancers, including multiple myeloma. A recent study indicates
   that GPR183 is highly expressed in cell lines resistant to HHT (an
   anti-tumor drug), suggesting that it may be associated with drug
   resistance in tumor cells^[156]33.

   At the compartment level, we found that more genes were activated, as
   there were 408 genes changed from compartment B state to A, and only
   321 genes changed from compartment A to B. GO enrichment analysis (BP)
   of genes with A-to-B switching event showed that these inactive genes
   were involved in “homophilic cell adhesion via plasma membrane adhesion
   molecules” term, which indicated that there was a decrease function for
   “Cell–cell adhesion via plasma-membrane” term. The genes in this term
   include PCDH9, LRFN3, IL1RAP, PCDHA6, PCDHA9, CDHA8, PCDHA7, PCDHA5,
   PCDHA4, PCDHA2, PCDHA1, PCDHA13, PCDHA11, PCDHA10, PCDHA12, and PCDHA3.
   A previous study investigated the interaction between B9/BM1 cells and
   osteoclasts and showed the possibility of tumor metastasis in bone
   marrow^[157]34. CML hematopoietic stem cells expressing IL1RAP can be
   targeted by chimeric antigen receptor-engineered T cells^[158]35, the
   inactivation of this gene might affect its function.

   At the TAD scale, consistent with the previous study^[159]5, we also
   found that MM genomes contain more TADs (~ 3627 in MM and 3197 in
   normal cells), and the average TAD size is smaller than in normal
   plasma cells (~ 0.70 Mb in MM and 0.72 Mb in normal cells). It is
   recommended that heterogeneity of cancer cells may contribute to more
   diverse 3D genomes within a cell population, increasing the detected
   TAD numbers^[160]5,[161]36. These 43 genes were found in MM-specific
   TADs. The genes were mainly related to “Human immunodeficiency virus 1
   infection” (including genes TNFRSF1B, FBXW11, NFATC1 and AP1S3), “Wnt
   signaling pathway” (including genes FBXW11, NFATC1, RSPO3, and FZD9),
   and “Cytokine–cytokine receptor interaction” (including genes TNFRSF1B,
   TNFRSF8, IL19, IL20, and IL24) KEGG pathways. Receptor activator of
   nuclear factor (NF)-κΒ ligand stimulation in multiple myeloma-derived
   osteoclasts induced elevated NFATC1, and selenoprotein W into the
   nucleus as compared to that in the control cells^[162]37. Our results
   showed that genes in the cytokine-cytokine receptor interaction KEGG
   pathway, including IL19, IL20, and IL24 changed their 3D structures,
   which could regulate the immune system and have inflammation
   effects^[163]38. Interleukin-20 (IL-20) is a pro-inflammatory cytokine
   with diverse angiogenic properties, serum IL-20 concentrations were
   found to participate actively in the pathophysiology of MM
   progression^[164]39. The presence of IL24 might affect tumor growth
   using a mouse model, but not as much as the therapeutic effect of
   HSV-tk, injection of IL24 reduced the tumor size^[165]40. At the loop
   level, these differentiated loops involved 113 genes, mainly in
   “AGE-RAGE signaling pathway in diabetic complications” and
   “Autophagy–animal” KEGG pathways. The genes were mainly MAPK10, EGFLAM,
   COL4A5, GAD1, ARSB, MNX1, GUCY2F, PABPC3, CIR1, CPLX1, WDR36, PITX2,
   NKX3-1, GATA4, and POM121L2. The identified chromatin interactions and
   differential gene expression patterns provide insights into the
   molecular mechanisms of MM, highlighting pathways such as immune
   response and cell adhesion that could be targeted therapeutically.
   These findings suggest that drugs interfering with these pathways or
   modulating the expression of key genes involved in MM pathogenesis
   might serve as potential therapeutic targets.

   In addition, by integrating 3D genome data, we found that genes such as
   NSG2, ALDH1A2, and CALCRL showed a higher expression in the MM
   patients. NSG2 is a member of the neuron-specific gene (NSG) family,
   which is specifically expressed in neurons^[166]17 and localized to the
   plasma membrane, the trans-Golgi network, and multiple endolysosomal
   compartments^[167]41. The aldehyde dehydrogenase 1 (ALDH1) family
   contains major enzymes that produce retinoic acid by the oxidation of
   all-trans-retinal and 9-cis-retinal, which mainly participates in
   biological functions such as cell differentiation, apoptosis, cell
   cycle arrest, and eventually^[168]42–[169]44. A previous study showed
   that ALDH1A2 expression is regulated by the epigenetic regulation of
   DNMTs, and subsequently might act as a tumor suppressor in ovarian
   cancer^[170]45. CALCRL is a G protein-coupled receptor that regulates
   the concentration of calcium ions in cells^[171]46. It could inhibit
   cell proliferation and angiogenesis^[172]47. CALCRL also contributes to
   the drug resistance in AML by controlling the ADM-CALCRL
   axis^[173]48,[174]49. ARHGAP24 showed a higher expression in the
   control. The Rho GTPase activating protein 24 (ARHGAP24) has been
   reported as a tumor suppressor in multiple cancers^[175]50. These
   results indicated that the NSG2, ALDH1A2, CALCRL, and ARHGAP24 genes
   might also have important potential functions during the multiple
   myeloma process. In addition, this study’s findings may be limited by a
   relatively small sample size of multiple myeloma patients with the same
   symptom, age and sex, which could affect the generalizability of the
   results. Meanwhile, recent studies^[176]51,[177]52 have shown that
   single-cell RNA sequencing is a crucial technique for investigating
   cellular heterogeneity, highlighting the need for its incorporation in
   future research. Further studies should consider there are more than 2
   MM patients and controls with the same sex, age and symptoms for
   analysis.

   Combining 3D genome, genome, and transcriptome analyses, we reveal that
   during MM development, multiple levels of alterations, such as spatial
   genome reorganization, occur accompanied by gene expression. The 3D
   genomic architecture identified in this study provides a foundation for
   exploring the spatial organization of genes and their regulatory
   elements in MM, which could lead to the development of targeted
   therapies that disrupt specific chromatin interactions crucial for the
   disease’s progression. For instance, the identification of key genes
   such as NSG2, ALDH1A2, and CALCRL, which show differential expression
   in MM patients, suggests potential avenues for therapeutic
   intervention, possibly through small molecule inhibitors or gene
   editing technologies. Furthermore, understanding the 3D genome’s role
   in antigen presentation, as indicated by the involvement of MHC class
   II protein complex assembly, could pave the way for novel
   immunotherapies that enhance the immune system’s ability to recognize
   and attack myeloma cells. However, since cancer types are diverse and
   alterations are heterogeneous, the phenomena observed in one may not
   hold in other cancer types. In the future, we need to investigate
   whether these observations are universal phenomena across cancer types.

Materials and methods

Plasma cell identification

   In this study, MACSprep Multiple Myeloma CD138 MicroBeads (Order no.
   130-111-744) have been developed to positively select CD138 + cells
   directly from whole blood. The operation method is modified according
   to the manufacturer’s manual, and the details are described in the
   following process.

   Adjust Cell Concentration: Ensure the specimen is qualified, with no
   hemolysis or coagulation; Use an automated hematology analyzer to
   measure the white blood cell count and calculate the suitable cell
   number: Suitable Cell Number = White Blood Cell Concentration
   (WBC) × Specimen Volume; Calculate the volume of cells to be added per
   tube (µL) = 1 × 10^6/Cell Concentration.

   Cell Surface Marker Staining: Label the tubes according to the
   experiment number and the antibody combination to be tested; Add the
   pre-prepared cocktail reagent to each tube for the chosen antibody
   combination; Mix the specimen well (at least 5 times by inversion), and
   add the calculated volume of cells to the bottom of the tube, mix well,
   and incubate in the dark for 15–20 min; Add 500 µl of lysing solution,
   place in the dark for 10 min until complete hemolysis occurs;
   Centrifuge at 2000 rpm for 3 min, discard the supernatant; Add 2 mL of
   1% newborn calf serum in PBS buffer to resuspend the cells, mix well,
   centrifuge at 2000 rpm for 3 min, and discard the supernatant; Add an
   appropriate amount (about 200–600 µL) of 1% paraformaldehyde fixative,
   resuspend the cells; Filter through a regular mesh with a pore size
   greater than 200 mesh, observe with the naked eye until no visible
   flocculent material or precipitate is present, and re-filter if
   necessary before preparing for machine detection.

   Intracellular Staining of Samples: Add 0.5 mL of 1 × FACS
   permeabilizing solution, mix well, and incubate at room temperature in
   the dark for 5 min; Add 2 mL of PBS, mix well, centrifuge at 2000 rpm
   for 3 min, and discard the supernatant; Mix the cells, add the standard
   amount of fluorescent McAb against intracellular antigens, mix well,
   and incubate at room temperature in the dark for 30 min; Add 2 mL of
   PBS, mix well, centrifuge at 2000 rpm for 3 min, and discard the
   supernatant; Add an appropriate amount (about 200–600 µL) of 1%
   paraformaldehyde fixative, mix well, and store in the dark at 2–8 °C,
   analyze within 24 h; Filter through a regular mesh with a pore size
   greater than 200 mesh, observe with the naked eye until no visible
   flocculent material or precipitate is present, and re-filter if
   necessary before preparing for machine detection.

   Plasma Cell Enrichment: CD138 Magnetic Bead Method: Use 1–2 mL of
   Running buffer to moisten the purple blood filter plug; Take 2–4 ml of
   blood into the filtration set and mark the blood level in a 15 mL tube;
   After centrifugation at 1400 r/min for 10 min, discard the supernatant
   and part of the settled red blood cells; Add Running buffer up to the
   mark, then add 100–150 µL of magnetic beads, mix slightly, and stand
   for 15 min; Use the magnetic bead column, first moisten with Running
   buffer, then filter the blood sample; Take Running buffer to the mark
   in a large tube, wash three times; Remove the magnetic bead column, add
   5 mL of Elution buffer, and quickly press down with a gas plug;
   Centrifuge at 1400 r/min for 10 min and discard the supernatant.

   Fish Detection: the enriched plasma cells were exposed to 2 × SSC,
   2 × SSC, 70% ethanol, 80% ethanol and 100% ethanol for 3 min, and then
   dried. Adding Probes, Hybridizing: Denaturing at 88 °C for 2 min,
   hybridizing at 45 °C for 2-16 h; 2 × SSC at room temperature for 1 min
   and in 0.3%NP-40/0.4 × SSC solution preheated at 68 °C for 2 min;
   Preheated distilled water at 37 °C for 1 min and naturally dried in the
   dark; Add 10μL of hybrid blue staining solution to the target area of
   the slide, and cover the slide. Select the appropriate filter to
   observe under the fluorescence microscope.

   Plasma Cell Freezing and Thawing, Freezing: Place the collection tube
   containing plasma cells into liquid nitrogen for rapid freezing and
   store it in a − 80 °C refrigerator for future use. Thawing: Place the
   plasma cell freezing tube in the − 80 °C refrigerator on ice.
   Continuously shake in a 37 °C constant-temperature water bath until the
   cells are thawed and ready for use.

Hi–C library construction, quality control and sequencing

   Cell cross-linking: the sample was fixed with formaldehyde, and the
   intracellular protein was cross-linked with DNA to preserve the
   interaction and maintain the 3D structure of the cell. Endonuclease
   digestion: DNA was digested by restriction endonuclease to produce
   sticky ends on both sides of the cross-linking. The restriction enzyme
   used in this project is DpnII. Terminal repair: Biotin-labeled bases
   are introduced to facilitate the purification and capture of subsequent
   DNA using the mechanism of terminal repair. Cyclization: cyclization of
   the DNA repaired at the end, and cyclization between the DNA fragments
   containing interactions, to ensure that the location of the interaction
   DNA is determined in the process of subsequent sequencing and analysis.
   DNA purification and capture: the DNA was de-crosslinked, and the
   purified DNA was broken into 300–700 bp fragments. The DNA fragments
   containing the interaction were captured by streptavidin magnetic beads
   and the library was constructed. After the construction of the library,
   the concentration of the library and the size of the inserted fragment
   (Insert Size) were detected by Qubit 2.0 and Agilent 2100,
   respectively, and the effective concentration of the library was
   quantified accurately by the Q-PCR method to ensure the quality of the
   library. After passing the library inspection, high-throughput
   sequencing was carried out on the Illumina platform, and the sequencing
   read length was PE150.

   Raw data (raw reads) of fastq format were first processed through
   in-house Perl scripts. In this step, clean data were obtained by
   removing reads containing adapter, reads containing ploy-N, and
   low-quality reads from raw data. At the same time, Q20, Q30, and GC
   content of the clean data were calculated. All the downstream analyses
   were based on clean, high-quality data. Clean Reads used BWA
   (Burrows-Wheeler Aligner, v0.7.17)^[178]53 to compare the two-terminal
   sequencing data with the reference genome sequences, respectively. The
   reads that can be compared are called Mapped Reads. Alignment
   efficiency refers to the percentage of Mapped Reads in Clean Reads.
   Using BWA^[179]53 with the command “mem -t 10-k 32,” to align the
   two-terminal sequencing data with the assembled genome sequence to
   obtain the only Read. Then using HiC-Pro (An optimized and flexible
   pipeline for Hi–C data processing, v2.10.0)^[180]54 to analyze the
   alignment results, identify the Valid Interaction Pairs and Invalid
   Interaction Pairs. We used HiC-Pro v2.10.0^[181]54 to obtain the
   corresponding standardized interaction matrix at various resolutions
   (10-, 20-, 100-kb, and 500-kb), and then calculate the Pearson
   correlation coefficient among the five samples.

Resolution analysis

   Sequencing depth of data determines the resolution of Hi–C data (the
   size of bin). Different resolutions should be used to study different
   sequencing depths of data. The resolution is calculated based on the
   method in reference^[182]20. That is, the number of interactive reads
   between each bin is arranged from high to low according to the
   different sizes of the bin. When 80% of the bin is arranged, the bin
   still covers the size of the smallest bin with more than 1000 reads,
   satisfying this condition, that is, the resolution of the Hi–C library.
   The Hi–C data are standardized by using HiC-Pro v2.10.0^[183]54
   software. Then the relationship between interaction frequency and
   genome linear distance is calculated by HiCdat^[184]55. The slope of
   the model is the corresponding IDEs (Interaction decay exponents). In
   the measured samples, it can be observed that the interactive signal
   decreases rapidly with the increase in distance. The mds algorithm of
   Pastis software is used to simulate the three-dimensional location of
   chromatin^[185]56. The software generates pdb (protein data bank) files
   based on Poisson distribution and uses pymol software to visualize the
   three-dimensional structure.

   The PCA (principal component analysis) was calculated by HiTC v1.24.0
   software^[186]57 with a bin size of 100 Kb. The PCA value was positive
   for the A compartment with high gene density and negative for the B
   compartment with low gene density. If there is a biological repetition
   for each sample, we merge the data from each library and identify A and
   B.

TAD and loop analysis

   The identification of topology-related domains (TAD) was carried out
   according to bin = 40 Kb using TadLib^[187]58 software. At the same
   time, the software HOMER (Hypergeometric Optimization of Motif
   EnRichment, v5.1)^[188]59 was used to calculate the ratio of IR
   (Inclusion Ratio), retain the TAD with IR > 1, and filter out the TAD
   with a length of less than 5 bins to obtain the final TAD. The DI
   values of each sample were calculated. Finally, the DI delta value of
   each TAD is further calculated (that is, the average difference of DI
   of 4 bins upstream and downstream of the TAD boundary). Then, the limma
   package was used to calculate the difference between the control and MM
   samples. Only when the FDR value of difference significance was less
   than 0.01, and the DI delta score of the two groups of samples was not
   all greater than 200 could it be considered as a TAD boundary of
   difference.

   Loop identification was using Juicer (v2.0)^[189]60 for analysis with
   bin size = 25 Kb, FDR ≤ 0.01. To compare the loop structures between
   different samples, we first merge the loops from each sample^[190]61,
   remove redundancy from the merged loops of all samples, and then
   calculate the original interaction values for each sample’s
   non-redundant loops as input for edgeR (v3.8.6)^[191]62. Loops with an
   FDR value less than 0.01, a significant P value less than 0.05, and an
   FC (Fold Change) greater than 1.5 are considered as differential loops.

Gene ontology analysis

   The gene ontology (GO) enrichment analysis was carried out by using
   clusterProfiler R software package, and the hypergeometric test was
   used to find the GO items which were significantly enriched compared
   with the whole genome background. KEGG (Kyoto Encyclopedia of Genes and
   Genomes)^[192]63,[193]64 pathway enrichment analysis was also
   conducted. We used KOBAS (KEGG Orthology Based Annotation System,
   v3.0)^[194]65 software to test the statistical enrichment of
   differential expression genes in KEGG pathways. We used clusterProfiler
   R packages to find KEGG pathways that are significantly enriched
   compared to the entire genome background.

Whole-genome sequencing experiments and analysis

   The whole genome DNA of the five MM patients was extracted and
   sequenced at an average 37.5 × depth through XTen (Illumina). The
   sequenced reads were mapped to the human reference genome (hg19) by the
   bwa-mem software^[195]66. Only uniquely mapped reads were used for
   downstream analysis. The Picard software
   ([196]http://broadinstitute.github.io/picard) was used to remove PCR
   duplicates. The total mapping rate is above 90% for each sample. The
   detection of SNP (Single Nucleotide Polymorphism) and small INDEL
   (small Insertion and Deletion) is mainly implemented by GATK software
   (The Genome Analysis Toolkit, v4.0)^[197]67. We next used Manta
   (v1.6.0)^[198]68 to identify genomic structural variation (SV).

RNA-seq experiments and analysis

   We use the TRIzol Reagent (Life Technologies, California, USA)
   instruction manual to extract total RNA from cells. Use the NanoDrop
   2000 (Thermo Fisher Scientific, Wilmington, DE) to measure the
   concentration and purity of RNA. Evaluate the integrity of the RNA
   using the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA,
   USA) with the RNA Nano 6000 Kit. The RIN (RNA integrity number) values
   ranged from 8.0 to 9.4. Indicating that all sample’s RNA were
   qualified, the quality meets the requirements of database construction,
   the total amount meets two or more conventional quantitative databases.
   Magnetic beads coated with oligo(dT) are used to selectively capture
   mRNA from total RNA. First-strand cDNA synthesis is then performed,
   followed by the synthesis of the second strand of cDNA. The overhangs
   at the ends of the cDNA are repaired to create blunt ends through the
   activity of exonucleases and polymerases. After the 3′ ends of the DNA
   fragments are adenylated, NEBNext adapters with a hairpin loop
   structure are ligated to them. The library fragments are then purified
   using AMPure XP magnetic beads from Beckman Coulter in Beverly, USA.
   Subsequently, 3 μl of USER Enzyme from NEB (USA) is added, and the
   mixture is incubated at 37 °C for 15 min. Before proceeding to PCR, the
   reaction is heated at 95 °C for 5 min to denature any remaining
   secondary structures. The PCR is then carried out using a high-fidelity
   DNA polymerase, along with universal PCR primers and index (X) primers.
   Finally, the PCR products are purified once more using AMPure XP
   magnetic beads, and the quality of the library is assessed using an
   Agilent Bioanalyzer 2100 system. The sequenced reads were mapped to the
   human reference genome (hg19) by TopHat (v2.1.1)^[199]69 and gene
   expressions were quantified by Cufflinks (v2.2.1)^[200]70. We used the
   RStudio software for the downstream statistical analyses. PUREE^[201]23
   was used to test the tumor purity of these MM samples.

Conclusion

   In summary, we have investigated the 3D genome of MM and analyzed the
   differences at compartment, TADs, loops, genomes, and gene expression
   levels to identify MM-related pathways and key genes. These findings
   extend our understanding of MM, which may implicate in clinical
   treatment and drug development for MM.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [202]Supplementary Material 1^ (558.6KB, pdf)
   [203]Supplementary Material 2^ (169.9KB, xlsx)

Author contributions

   Conceptualization, K.Z.,Y.L., and X.F.; writing—original draft
   preparation, K.Z., M.C.(Mengsi Chen), M.C.(Ming Chen), Y.W. and H.L.;
   writing—review and editing, Y.L., X.G., L.L.,L.T.,X.L., D.H., and X.F.;
   Supervision, X.L., D.H., and X.F.; funding acquisition, K.Z.; All
   authors have read and agreed to the published version of the
   manuscript.

Funding

   This work was supported by The Science and Technology Department of
   Sichuan Province (No.2020YJ0438).

Data availability

   Data availability: The accession number for the WGS, RNA-seq, and HiC
   data sets generated in the study is HRA007587via GSA
   (https://ngdc.cncb.ac.cn/gsa/).

Declarations

Competing interests

   The authors declare no competing interests.

Ethics statement

   All methods carried out in accordance with the relevant guidelines and
   regulations. This study has been approved by the Ethics Committee of
   Chengdu First People’s Hospital, ethics number ZXKY No.002, in 2020.
   All patients and blood donors in this study have signed informed
   consent forms by themselves or their guardian(s).

Footnotes

   Publisher’s note

   Springer Nature remains neutral with regard to jurisdictional claims in
   published maps and institutional affiliations.

References