Revealing Genomic Profile That Underlies Tropism of Myeloma Cells Using Whole Exome Sequencing

Background. Previously we established two cell lines (SNU_MM1393_BM and SNU_MM1393_SC) from different tissues (bone marrow and subcutis) of mice which were injected with single patient's myeloma sample. We tried to define genetic changes specific for each cell line using whole exome sequencing (WES). Materials and Methods. We extracted DNA from SNU_MM1393_BM and SNU_MM1393_SC and performed WES. For single nucleotide variants (SNV) calling, we used Varscan2. Annotation of mutation was performed using ANNOVAR. Results. When calling of somatic mutations was performed, 68 genes were nonsynonymously mutated only in SNU_MM1393_SC, while 136 genes were nonsynonymously mutated only in SNU_MM1393_BM. KIAA1199, FRY, AP3B2, and OPTC were representative genes specifically mutated in SNU_MM1393_SC. When comparison analysis was performed using TCGA data, mutational pattern of SNU_MM1393_SC resembled that of melanoma mostly. Pathway analysis using KEGG database showed that mutated genes specific of SNU_MM1393_BM were related to differentiation, while those of SNU_MM1393_SC were related to tumorigenesis. Conclusion. We found out genetic changes that underlie tropism of myeloma cells using WES. Genetic signature of cutaneous plasmacytoma shares that of melanoma implying common mechanism for skin tropism. KIAA1199, FRY, AP3B2, and OPTC are candidate genes for skin tropism of cancers.


Introduction
Multiple myeloma is a malignant B-cell disorder characterized by proliferation of atypical plasma cells in bone marrow [1] with or without the presence of monoclonal immunoglobulin protein in serum and/or urine [2]. Multiple myeloma has correlation with plasmacytoma, which is a mass of plasma cells found outside of bone marrow [3] that needs medical intervention with radiotherapy [4] or chemotherapy. While multiple myeloma frequently accompanies plasmacytoma at the time of diagnosis, plasmacytoma precedes multiple myeloma in some cases. The disease entity called solitary plasmacytoma exists in 4% of plasma cell tumors [5,6], and approximately 40-50% of patients with solitary plasmacytoma will develop multiple myeloma [7]. Hence, plasmacytoma is an early form or an accompanying disease of myeloma, and the data regarding the "clinical behavior" of plasmacytoma are quite accumulated. However, not much is known about the cellular biology of plasmacytoma per se.
Fortunately, we previously succeeded in establishing two cell lines with different tropism from a single patient [8]. Briefly, 14 weeks after injecting mononuclear cells obtained from multiple myeloma patient's bone marrow in a NRG/ SCID mouse via tail vein, tumor developed at subcutis of the mouse. The engraftment of myeloma cells into mouse bone marrow (BM) was also observed. After separation and culturing cells from subcutis and BM, we succeeded in the establishing of two cell lines (SNU MM1393 SC and SNU MM1393 BM) from a single patient with different tropism. These two cell lines showed different biologic behavior. In contrast to SNU MM1393 BM, cell proliferation of SNU MM1393 SC was IL-6 independent. SNU MM1393 BM 2 International Journal of Genomics and SNU MM1393 SC showed a high degree of resistance against bortezomib compared to U266 cell line. SNU MM1393 BM had the greater lethality compared to SNU MM1393 SC.
From genetic perspective, we believe that unique somatic mutations may allow tumor cells to adapt and survive in tumor microenvironment. In other words, there may be specific genetic changes that contribute to tumor tropism. Hence, in this study, we tried to characterize genomic profile specific for SNU MM1393 BM and SNU MM1393 SC to understand genetic background of difference in these cell lines. We were particularly interested in genetic changes specific for SNU MM1393 SC, because we thought these SNU MM1393 SC specific genetic changes would reveal genetic background for plasmacytoma which exhibit tropism for extramedullary space. To find these genetic changes, we performed whole exome sequencing (WES) using DNA of both cell lines. As it is well known, WES allows comprehensive characterization of genomic changes in individual tumors [9].

DNA Preparation and Whole Exome
Sequencing. DNA was extracted from two cell lines (SNU MM1393 BM and SNU MM1393 SC). QuickGene DNA whole blood kit S (Kurabo Industries Ltd., Japan) was used to extract DNA according to the manufacturer's recommendations. For WES, we sequenced exome using the Solexa sequencing technology platform (HiSeq2000, Illumina, San Diego, CA, USA) following the manufacturer's instructions. We randomly sheared 3 ug of genomic DNA using Covaris System to generate about 150 bp inserts. The fragmented DNA was end-repaired using T4 DNA polymerase and Klenow polymerase, and Illumina paired-end adaptor-oligonucleotides were ligated to the sticky ends. We analyzed the ligation mixture by electrophoresis on an agarose gel, sliced and purified fragments with 200-250 bp sizes. Purified DNA library was hybridized with SureSelect Human All Exon V4 probes set (Agilent, Santa Clara, CA, USA) to capture 50Mb targeted exons following manufacturer's instruction. We prepared the HiSeq2000 paired-end flow cell to the manufacturer's protocol using captured exome library. Clusters of PCR colonies were then sequenced on the HiSeq2000 platform using recommended protocols from the manufacturer.

Alignment of FASTQ File
. FASTQ files were aligned to human reference (human g1k v 37.fasta) by using the Burrows-Wheeler aligner (BWA-0.7.5) [10] to make SAM file. SortSam in Picard-tools-1.68 was used to convert to BAM file and sort with chromosome and went through a PCR duplicate marking process, which enables the Genome Analysis Toolkit (GATK-1.6.5) to ignore duplicates in subsequent processing [11]. We performed a local realignment prior to recalibration, which gives the most accurate quality scores for each sample. Recalibration was performed to increase recalibration accuracy. A default value was set for all processes.

Variant Calling and Annotation.
For single nucleotide variant (SNV) and small indel calling, we used Varscan2 (http://genome.wustl.edu/) [12]. Because germline control was absent, we used pilup2snp command for single sample calling using human g1k v 37 as reference genome. We set the following options: (1) minimum coverage above value of 20, (2) minimum variant allele frequency above value of 0.1, and (3) default value below 0.05 other option value set as default values. To select unique mutation, we performed comparison between two calling results. For functional annotation and prediction of variant effect, we used ANNOVAR [13] with Polyphen [14] database version 2.2.2.
The rate of transversion and transition in the coding region was different between the two cell lines. While transversion was dominant event in SNU MM1393 BM cell line, transition was dominant event in SNU MM1393 SC. Absolute transversion rate was much higher in SNU MM1393 BM (65.5%) than SNU MM1393 SC (34.0%) (Figure 1(d)).   (Table 1). Then, the frequencies of these SNVs were investigated in open source data of various cancers (http://www.cBioportal.org). When this analysis was performed, SNVs which are unique for SNU MM1393 SC were frequently detected in melanoma (52.8%) and uterine cancer (49.0%). On the other hand, SNVs those which are unique for SNU MM1393 BM were frequently found in ovarian cancer (74.7%) and bladder cancer (72.8%) (Figure 2).

Monte Carlo Simulation and Pathway Analysis. Using
Monte Carlo Simulation (Ratio of dNs/dS), we examined distribution statistics of SNVs found in two cell lines, respectively. It is believed that nonrandomness of SNV distribution is related to the functional importance of those SNVs. When this analysis was performed, ratio for SNU MM1393 BM was 2.8 ( = 0.14), while it was 1.1 for SNU MM1393 SC ( = 0.07). Hence, SNV distribution in both cell lines was random with cut-off value of 0.05. Our results indicated that unique nonsynonymous mutations of SNU MM1393 SC seemed biologically more neutral than those of SNU MM1393 BM although they were statistically insignificant.
In KEGG pathway analysis of unique somatic mutation from both cell lines, the result of pathway analysis showed both similarities and differences. Chemokine signaling pathway and chemokine-chemokine interaction pathway were showed in two cell lines at the same time. While tight junction and adherens junction were shown to be involved only in SNU MM1393 SC, SNU MM1393 BM cell line showed diverse disorder in pathway such as cell cycle, Wnt signaling pathway, and MAPK signaling pathway (Figure 3).

Discussion
In this report, we analyzed genomic signature of two cell lines derived from single patient. These cell lines have different tropism, and SNU MM1393 SC is plasmacytoma which has tropism for skin. Thus, we thought this study would reveal genetic changes that are related to skin tropism and plasmacytoma. In fact, recent study reported the process of segregation of genetic changes in tumor cells during clonal expansion [16]. As far as we know, this is the first WES study comparing the genomics of myeloma cells in bone marrow and plasmacytoma. Varscan2 which was used in our study for somatic variant calling is one of the most commonly used tools to detect somatic mutation along with MuTect [12,17]. In fact, many papers already published in the field of cancer genomics already used Varscan2 [18][19][20].
As expected, most of SNVs (96.5%) overlapped between the cell lines, suggesting their common originality from a single patient. And the most coding sequence SNVs in both SNU MM1393 BM and SNU MM1393 SC were neutral with respect to adaptation and cancer cell growth. This is also supported by the outcome of Monte Carlo Simulation, where distribution of SNVs did not show significant nonrandomness. It has been suggested that bone marrow microenvironment strictly regulates the growth of cell via dynamic interplay among hematopoietic cells [21]. And our results indicate that only a small subset of nonsynonymous SNVs in cancer are affected by selection, making it possible to interpret SNV trends as reflection of underlying mutational processes.
The most interesting finding in our data is related to the genetic changes unique for SNU MM1393 SC. We conjectured in the planning of this study that genetic change unique for SNU MM1393 SC would be related to skin tropism and formation of plasmacytoma. And, as expected, genomic signature that is unique to SNU MM1393 SC had more than 50% of overlapping with melanoma which is a primary skin tumor. This finding highly coincides with our conjecture and SNVs such as KIAA1199, FRY, AP3B2, and OPTC may be the very gene related to skin tropism of cancers. These 4 are the genes that are mutated in SNU MM1393 SC and are frequently found in melanoma samples. In fact, it has been known that there are common genetics between melanoma and plasmacytoma such as CDKN2A germline mutation in prenext generation sequencing (NGS) era [22].
One more noticeable finding in our study is that the number of SNVs was higher in SNU MM1393 BM than in SNU MM1393 SC and transversion rate was higher in SNU MM1393 BM than in SNU MM1393 SC. Also dynamic gene to gene interaction in SNU MM1393 BM was more complex than that in SNU MM1393 SC. Along with this phenomenon and from the previous report that the number and pattern of somatic SNVs determine pathway underling cancer [23], we think that it would be biologically simpler to form a tumor in subcutis than to form a tumor in bone marrow in animal xenograft model used in our experiment. We think this is related to the fact that growth and differentiation of cancer cells in bone marrow are strictly regulated by dynamic interplay among various hematopoietic cells, compared to subcutis.   ARSA  AP3B2  CTNS  EPB41L2  EPB41L3  NKX2-2  FCGR3B  HSD3B1  PSTK  XCL2  SMAD6  ALDOA  PTPRM  OR5AP2  S1PR5  GABRE  PRKX  PPM1D  CCR5  CDC27  ITGA6  CACNA1I  ARHGAP5  OR2T33   AK2  HLA-DRB1  SLC27A6  IARS2  COQ2  AOX1  ABCA6  ALDH5A1  DHRS4  B3GAT1  PTPLA  EGLN1  PANK2  PIWIL4  CNOT6  B4GALNT1   ACER2  OXA1L  DHODH  SPRED1  CHAT  TAF6L  NRP1  KISS1R  ZIC2  IGFBP3  DNAJC25  ITGAE  IL17RA  In fact, we had difficulties in the analysis of WES data due to the lack of germline reference DNA in this study. Because we used public germline database as reference in the calling of SNVs, the number of SNVs was very large compared to the previous multiple myeloma genomic studies. Moreover, it is well known that the majority of mutations observed in cancer sequencing studies are believed to be passenger mutations having little impact on the cancer cell [24]. To overcome this issue, we did not focus on the specific genetic changes found in our study. Rather, we analyzed the "pattern" of genetic changes found in our study and compared them with public database so as to find out genetic clues underlying tropism for International Journal of Genomics 7 skin and plasmacytoma. And this is also the reason why we focused on the "unique" SNVs for both SNU MM1393 BM and SNU MM1393 SC instead of whole genomic picture of SNU MM1393 BM and SNU MM1393 SC.