Transcriptome Analysis of Mesenchymal Stem Cells from Multiple Myeloma Patients Reveals Downregulation of Genes Involved in Cell Cycle Progression, Immune Response, and Bone Metabolism

A growing body of evidence suggests a key role of tumor microenvironment, especially for bone marrow mesenchymal stem cells (MSC), in the maintenance and progression of multiple myeloma (MM), through direct and indirect interactions with tumor plasma cells. Thus, this study aimed to investigate the gene expression and functional alterations of MSC from MM patients (MM-MSC) in comparison with their normal counterparts from normal donors (ND-MSC). Gene expression analysis (Affymetrix) was performed in MM-MSC and ND-MSC after in vitro expansion. To validate these findings, some genes were selected to be evaluated by quantitative real time PCR (RT-qPCR), and also functional in vitro analyses were performed. We demonstrated that MM-MSC have a distinct gene expression profile than ND-MSC, with 485 differentially expressed genes (DEG) - 280 upregulated and 205 downregulated. Bioinformatics analyses revealed that the main enriched functions among downregulated DEG were related to cell cycle progression, immune response activation and bone metabolism. Four genes were validated by qPCR - ZNF521 and SEMA3A, which are involved in bone metabolism, and HLA-DRA and CHIRL1, which are implicated in the activation of immune response. Taken together, our results suggest that MM-MSC have constitutive abnormalities that remain present even in the absence of tumors cells. The alterations found in cell cycle progression, immune system activation, and osteoblastogenesis suggest, respectively, that MM-MSC are permanently dependent of tumor cells, might contribute to immune evasion and play an essential role in bone lesions frequently found in MM patients.

previous treatment for the disease, i.e., no chemotherapy, no corticosteroids, no immunomodulators, no proteasome inhibitors, or bisphosphonates, were successfully enrolled in this study and allocated in the case group. The clinical laboratory characteristics of MM patients at the diagnosis are shown in Table 1. Seven BM normal donors for allogeneic stem cell transplantation, not matched by age or gender, were also included in the study and allocated in the control group. As an additional control, the HS-5 human bone marrow normal stromal cell line (ATCC, Manassas, VA, USA) was also used in some experiments.
Isolation, expansion and characterization of MsC. BM samples from normal donors (n = 7) and newly diagnosed MM patients (n = 19) were harvested from patients' iliac crest. Then, BM mononuclear cells were isolated using Ficoll-Paque PLUS (GE Healthcare, Little Chalfont, Bucks, GBR), according to the manufacturer's instructions. Finally, MSC were sorted by Magnetic-Activated Cell Sorting (MACS) methodology, using CD105 + as a positive marker (Miltenyi Biotec, Bergisch Gladbach, DEU). MSC expansion was performed on αMEM with GlutaMAX and nucleoside medium 36 , supplemented with penicillin (100 U/mL)/streptomycin (100 μg/mL), fungizone (2.5 μg/mL) (All Gibco, Carlsbad, CA, USA) and 10% fetal bovine serum (FBS) (Vitrocell, Campinas, SP, BRA). MSC were incubated at 37 °C, 5% CO 2 and high humidity. During passage zero, the cells were fed twice a week, replacing only 50% of culture medium volume, until the cells reached 80% confluence or up to 21 days. Then, MSC were detached from the culture plastic with the aid of 0.25% trypsin-EDTA reagent (Gibco, Carlsbad, CA, USA), counted by tripan blue exclusion method and seeded again. From the first The osteoblastic differentiation of MM-MSC (n = 4) and ND-MSC (n = 4) was performed in technical duplicates in 12-well microplates, using StemPro Osteogenesis Differentiation kit (Gibco, Carlsbad, CA, USA) and following the manufacturer's instructions. Cells were fed twice a week, and on days 7, 14 and 21, the differentiation medium was removed and frozen at −80 °C for further analysis. In order to confirm the cell differentiation by a quantitative methodology, osteocalcin measurement was performed on the differentiation media collected from the cases and controls on days 7, 14 and 21, using the Human Osteocalcin Quantikine ELISA kit (R&D Systems, Minneapolis, MN, USA), according to the manufacturer's instructions.
Microarray hybridization and data acquisition. For gene expression analysis, RNA samples were extracted from MM-MSC (n = 4) and ND-MSC (n = 4) using RNeasy Mini kit (Qiagen, Valencia, CA, USA). For each sample, three independent RNA extractions were performed. RNA quantification and purity analysis were performed on the NanoDrop ® ND-8000 UV spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). In addition, RNA integrity was verified by 1% agarose gel electrophoresis, stained with ethidium bromide.
The microarray platform used was GeneChip ® Human Exon 1.0 ST Array (Affymetrix, Santa Clara, CA, USA). On this platform, there are about 1.4 million probe sets, being approximately 4 probes per exon and 40 probes per gene. Microarray images were obtained using GeneChip Scanner 3000 7 G, and data were quantified using Affymetrix GeneChip Command Console ® Software (both Affymetrix, Santa Clara, CA, USA), generating CEL files containing the raw data.
Microarray data preprocessing and analysis. Microarray raw data preprocessing and identification of differentially expressed genes (DEG) were performed using the AltAnalyze software www.altanalyze.org 38 . Once the analysis parameters were set, the raw data, in CEL format, was preprocessed by the Robust Multi-array Analysis (RMA) method 39 , which includes background correction, quantile normalization and summarization of the probes into specific probe sets. The probe sets were defined as differentially expressed between the groups when the p-values, corrected by the False Discovery Rate (FDR) method 40 , were less than 0.05 and the fold-changes (difference in expression in the case group versus control group) were greater than 1.5, in module. Differentially expressed probe sets were annotated for the purpose of identifying which genes they represent. To ensure that there was no great variability among within-condition samples, the coefficients of variation (CV), of the normalized gene expression values in log2, were calculated and, arbitrarily, the CV cut-off criteria less than 15% was established to consider a gene consistent. The microarray data, discussed in this article, have been deposited in NCBI's Gene Expression Omnibus, and can be accessed through GEO Series accession number (ref GSE113736).

Bioinformatics analyses workflow.
After identification of DEG, we performed the bioinformatics analyses in order to extract relevant biological information among these genes. Gene Co-Expression Network Analysis. Gene co-expression network construction and additional analyses were performed using Cytoscape 3.5.1 software 41 , and three of its plug-ins. First, the GeneMANIA plug-in 42 was used to generate the network, through the prediction of interactions among DEG, based exclusively on data published in the literature concerning co-expression. Then, another plug-in, CentiScaPe 43 was used to calculate centrality measures of the genes (nodes) belonging to the constructed network. In our study, the calculated centrality measures were degree and betweenness, which represent, respectively, the number of connections of a node, i.e., the number of interactions of a gene with other genes in the network, and the number of shortest paths that pass through a node to connect other pairs of nodes. Lastly, GLay plug-in 44 was used to find modules, also known as communities or clusters, which means groups of highly interconnected genes in the network.
Identification of high-hubs, hubs and bottlenecks. The calculated degree and betweenness values were used to construct a scatter plot, using GraphPad Prism 7.0 statistical software (GraphPad Software, San Diego, CA, USA). The scatter plot allows categorization of nodes in high hubs, hubs, and bottlenecks, as previously described by Azevedo et al. 45 . In summary, by dividing the plot into quadrants, the genes located in the upper right quadrant represent the high hubs (high degree and betweenness values), whereas the genes located in the lower right quadrant represent the hubs (genes with high degree and low betweenness values), and, finally, genes located in the upper left quadrant represent bottlenecks (genes with high betweenness and low degree values  . Then, all DNA samples were diluted to a final concentration of 50 ηg/μL. MSC telomere length was determined by multiplex real time qPCR, as previously described by Cawthon 56,57 , with minor modifications. In summary, this methodology consists in determining the relative ratio (T/S) between the telomere region copy number (T) and a single copy gene (S), using a relative standard curve. In our study, we chose the ALB gene as the single copy gene. T/S ratio for each sample is proportional to the mean telomere length. All experiments were performed in triplicate and our CV inter-assay was around 13.04%.
Cell cycle analysis. MM-MSC and ND-MSC frequencies distribution among cell cycle phases were evaluated in the BD FACSCanto II flow cytometer, using propidium iodide reagent (both Becton, Dickinson and Company, Franklin Lakes, NJ, USA). The results were analyzed using ModFit LT software (Verity Software House, Topsham, ME, USA). statistical analyses. All statistical analyses were performed on IBM SPSS Statistics 20.0 software (IBM Corporation, Armonk, NY, USA), adopting α = 5% significance level. All graphs were plotted in GraphPad Prism 7 software (GraphPad Software, San Diego, CA, USA) and the results are shown as mean and standard deviation (SD). In order to evaluate the group effect (MM-MSC versus ND-MSC) over time (7, 14 and 21 days) on the measurements of the continuous variable osteocalcin, we used the Generalized Estimating Equation (GEE) with gamma distribution. Mann-Whitney U test was used to perform comparison among groups regarding relative gene expression by RT-qPCR. Additionally, to evaluate group effect on the continuous dependent variable mean telomere length (T/S), we used the independent t-test, as the probabilistic distribution of this variable was considered normal (p = 0.01, Kolmogorov-Smirnov test). We also assumed the homogeneous variance distribution between groups, since Levene's test showed no significant difference between group variances (F = 0.053 and p = 0.819). Lastly, to investigate the existence of an association between the group (MM-MSC versus ND-MSC) and the relative frequency of cells in the different cell cycle phases (G0/G1, S and G2/M), we performed the Fisher's exact two-tailed test, since some expected frequencies were less than five. Principal component (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analyses were implemented in the R software in order to perform dimensionality reduction and assess how the samples group to each other.

Distinct gene expression profiling between MM-MsC and ND-MsC. After microarray data
pre-processing and establishment of cut-off criteria (adjusted p-value < 0.05 and fold-change >1.5), we found 485 DEG between MM-MSC and ND-MSC, including 280 upregulated and 205 downregulated genes. About 50% of them correspond to protein-coding genes, while the other half comprised genes encoding for long non-coding RNAs, small nuclear RNAs, and small nucleolar RNAs. Arbitrarily, we considered that genes with the CV values less than 15%, i.e., with low variability within-condition, were consistent. Only 30 genes in the MM-MSC group and 28 genes in the ND-MSC were considered not consistent out of the 485 DEG, with an overlap of 12 non-consistent genes between the groups. The majority of the non-consistent genes were non-protein coding genes. PCA and t-SNE analyses were performed to assess how the samples grouped to each other. Interestingly, both analyses were able to separate the samples into the studied groups, showing a good degree of similarity between the samples belonging to the same group (Figs S1 and S2). The first two PCA components (PC1 and PC2) were able to explain 23.4% and 11.9% of the data, respectively (Fig. S3).
Network and functional enrichment analyses reveal the functional gene signature of MM-MsC. The 485 DEG were used to build a gene co-expression network, where genes are represented as nodes, and their relationship, in this case co-expression, is represented as edges connecting the genes. The Cytoscape plug-in GeneMANIA was used to construct the network, using only interactions from the co-expression category. After filtering out the genes that were not connected in the network, we observed 195 genes (nodes) and 1515 interactions (edges) in the final network, with 31 genes upregulated, and 164 downregulated.
The degree and betweenness values of the 195 nodes were calculated, using the CentiScaPe plug-in, and a scatter plot was built to visualize the relationship between degree (x-axis) and betweenness (y-axis) for each node. The plot was divided into four quadrants by applying an arbitrary cut-off of degree ≥29 and betweenness ≥671.5 (Fig. 2). The 20 DEG with the highest degree (hubs) values and the 20 DEGs with the highest betweenness (bottlenecks) values were highlighted (total of 36 downregulated genes). Among the hubs, half of them participate directly or indirectly in cell cycle progression. Additionally, some other genes are involved in functions related to the immune response, such as antigen processing and presentation via MHC class II, and regulation of the complement system activation ( Table 2). In parallel, bottleneck genes participate in cell cycle regulation or immune response. Besides, some genes are involved in osteoblastogenesis (SEMA3A, BICC1, CHRLD1, and ZNF521) and HAS1 is involved in Waldenstrom's macroglobulinemia (Table 3), a post-follicular B-cell lymphoproliferative disorder also associated with M-protein production (IgM). Finally, three genes with high degree and betweenness (high hubs) are involved in cell cycle progression, while one high-hub participates in the regulation of complement system activation (Table 4). We also performed the identification of functional modules within the gene co-expression network, i.e., the localization of groups of highly connected nodes inside the original network (Fig. 3). Enrichment analysis of the functional modules was performed using the Enrichr and the DAVID softwares. Since the results obtained from both softwares were very similar, i.e., there was a great overlap between the terms of GO and the biological pathways of KEGG significantly enriched, we decided to show only the results generated by Enrichr. Overrepresented GO biological processes and KEGG pathways are represented in Figs 4 and 5 for module 1, and in Figs 6 and 7 for module 3. We did not find any biological processes nor biological Real-time Rt-qpCR validation. The microarray data were validated at the mRNA level by RT-qPCR for SEMA3A, ZNF521, HLA-DRA, CHI3L1, NEIL3, and GPC6. The genes HLA-DRA, CHI3L1, and ZNF521 were downregulated in 69%, 69% and 62% of MM-MSC, whereas genes SEMA3A, NEIL3 and GPC6 were downregulated in 54%, 38%, and 31% of MM-MSC, respectively (Fig. 8). Although all genes were downregulated in tumor samples when compared with normal controls, we only found a statistically significant difference for the ZNF521 gene (p = 0.046) (Fig. 9).

Multiple myeloma effects on MM-MSC telomere length. Telomere length comparison between
MM-MSC and ND-MSC are showed in Fig. 10. As expected, MM-MSC presented lower telomeric length (mean = 0.97, SD = 0.11) than ND-MSC (mean = 1.04, SD = 0.11). However, the independent t test for the two samples showed that the difference found was not statistically significant (t (24) = 1.578, p = 0.128).

Discussion
Gene expression profiling analysis of MM-MSCs compared to ND-MSCs revealed 485 DEG, being 280 upregulated and 205 downregulated. When we built the co-expression network, there was a reversal and the downregulated genes became more represented in the network than the upregulated (164 versus 31), which was expected, since among the 280 upregulated genes, only 49 were protein-coding genes, whereas, of the 205 downregulated ones, 170 were protein-coding genes. Further exploration of DEG with different bioinformatics tools, showed that the most relevant enriched pathways and functions were among the downregulated genes, especially for those involved in cell cycle progression, immune response activation, and osteoblastic function and maturation. We have validated part of these findings through real-time quantitative PCR methodology, and additional in vitro functional assays. Six downregulated target genes, whose functions were mainly related to cell cycle progression, immune response activation, and osteoblastic function and maturation, were selected for validation at the mRNA level. Of the six target genes selected, a statistically significant difference was detected only for the ZNF521 (p = 0.048), which was downregulated in 62% of the MM-MSC evaluated samples, validating the microarray findings. In line with these results, a tendency to down expression was also observed for the HLA-DRA and CHI3L1 genes (both 69%), and for SEMA3A gene (54%). Regarding the functional in vitro assays, we chose to validate the cell cycle progression, since this function has already been reported by other authors as being impaired in MM-MSC, leading these cells to enter into an early cellular senescence process when compared . However, the difference found was not statistically significant, suggesting that the alterations in mRNA level related to the cell cycle are not related to a shortening of telomeres, or the casuistic was insufficient to detect this difference due to the biological variability among individuals. Our results suggest that MM-MSCs have a distinct gene expression profiling in comparison with ND-MSCs. Interestingly, these cells showed hundreds of DEG even after in vitro expansion in the absence of tumor cells, confirming previous findings [32][33][34]58 (Table 5). However, as discussed by André et al. 33 , the variability of the results found in these studies is relatively high, possibly due to several factors, including but not limited to: methodological differences prior to microarray execution -from the local and the form of MSC isolation, as well as its method of cultivation -up to differences in the microarray platform chosen, as well as the statistical and bioinformatics approaches used for the pre-processing of the raw data and identification of the DEG. However, despite the differences, it is still possible to make relevant comparisons among these studies, and to raise quite pertinent hypotheses.
The first work to evaluate the overall gene expression profiling of MM-MSC, in comparison with ND-MSC, was published in 2007 by Corre et al. 32 . Among the 183 DEG/Probesets found, 59 were classified by the authors as belonging to the category of tumor microenvironment, comprising functions such as cellular communication, receptor signaling molecules, extracellular matrix, and secretory molecules. Additionally, they highlighted 40 genes (20 upregulated and 20 downregulated) as being essential for MM, of which four were also found in our study, one gene upregulated -ANGPTL4 -and three with diminished expression -NPR3, TNFRSF19, and FBLN1. The gene ANGPTL4 was also found downregulated in MM-MSC in three previous independent studies [32][33][34] . More recent publications have demonstrated its multiple roles in osteolytic lesions 59 , and MM bone disease 60 , for example, through the promotion of osteoclast-mediated bone resorption, cartilage degradation, and angiogenesis. The NPR3 and FBLN1 genes, among other functions, participate in bone formation 61,62 , i.e., the reduction of their expression can potentially contribute to the development of osteolytic lesions frequently found in patients with MM. Finally, the TNFRSF19 gene, a member of the TNF receptor superfamily, appears to mediate caspase-independent cell death (Gene database, NCBI).  In line with our results, André et al. 33 also found enriched categories related to the cell cycle, such as, "M-phase", "DNA replication", "cell cycle regulation", etc, among downregulated genes, like those results found by Wagner et al. (2008) in a study with ND-MSC in replicative senescence 63 . In our study, one of the functional modules detected through Cytoscape was composed mainly of genes involved in different pathways and biological processes related to cell cycle progression. In addition, André et al 33 . also demonstrated that MM-MSCs have a lower proliferative rate, higher cell size, β-galactosidase increased activity, retention of cells in the S phase of the cell cycle, and secrete a senescence-associated molecule profile. Thus, the authors hypothesized that MM-MSCs appear to become senescent earlier than ND-MSCs. It is important to highlight that in this study, MSCs were expanded in vitro in monoculture, i.e. in the absence of other cells including MM plasma cells. Berenstein et al. 64 demonstrated that MM-MSC, when cultured in the absence of tumor cells, accumulate in the S-phase of the cell cycle, increase β-galactosidase activity, and increase the expression of microRNAs associated with senescence. Moreover, they also demonstrated that MM-MSC co-cultivation with a MM cell line (KMS12-PE) is able to reverse, at least partially, the senescence phenotype of MM-MSC 64 . Contributing to these findings, the study of Garcia-Gomez et al. 34 , which evaluated the gene expression profiling of MM-MSC expanded in monoculture and co-culture with the MM cell line MM.1 S, observed that, after co-cultivation, some genes related to cell cycle progression become upregulated.
With regard to the other functional module identified through Cytoscape after performing functional enrichment analysis, we observed that the majority of the enriched categories was related to different aspects of the immune response, including antigen processing and presentation via MHC classes I and II, T cells activation, and immune response triggered by inflammatory cytokines. All genes belonging to these categories were downregulated in MM-MSC, suggesting that these biological pathways and processes could figure as a possible mechanism of immune escape. The immune system of MM patients is highly impaired 65    showed that in BM of MM patients there is an increased expression of Treg cell markers, FOXP3 and CTLA4 genes, suggesting their accumulation in the BM microenvironment of MM patients, and a possible mechanism of tumor evasion from the immune system 66    Regarding bone formation and resorption, the ZNF521 gene, which was downregulated in MM-MSCs, plays a key role in bone metabolism -this gene acts inhibiting the differentiation of osteoblast progenitors, through binding to RUNX2 pro-differentiation transcription factor, and simultaneously promoting maturation and correct function of mature osteoblasts 75 . In agreement with these data, the RUNX2 gene was also downregulated in MM-MSC (FC = −1.65, gene expression evaluated only by microarray), possibly contributing to the imbalance of the bone metabolism observed in patients with MM. In addition, the downregulation of SEMA3A, which was classified as a bottleneck in the co-expression network, may also contribute to bone metabolism imbalance, since it plays an important role in osteoblastogenesis, inhibiting osteoclastic differentiation and stimulating osteoblastic differentiation 76 . This gene also acts as an inhibitor of angiogenesis in endothelial cells, and a study reported that the loss of its inhibitory capacity may contribute to the transition from MGUS to the active form of MM 77 . Finally, the BICC1 gene, also downregulated in the MM-MSC and classified as a bottleneck, is a genetic determinant of osteoblastogenesis and mineral bone density 78 .
The main limitation of our study is the small number of MM patients enrolled (n = 19). All of them were diagnosed in the same public hospital which receives most of MM new cases in stage III (84% in this study) and medical emergency, such as spinal cord compression, renal insufficiency, hypercalcemia or bone fractures, when immediate therapeutic interventions with corticosteroids and bisphosphonates are necessary, making patients ineligible for gene expression studies. Another limitation was that the controls were not age-matched to cases. In general, transplant normal donors were younger than MM patients. This situation could raise the hypothesis that the early senescence profile of MM-MSC, or the other differences detected through bioinformatics analysis, could be artificially created by the lack of age matching. Magalhães et al. 79 conducted a meta-analysis of microarray studies that evaluated aging-related genes, and they identified 74 genes with higher levels of evidence. Of the 485 DEG in the MM-MSC compared to the ND-MSCs identified in this study, there was an overlap of only two genes within the 74 reported by Magalhães et al. 79 . The genes were SERPING1 and S100A6, which are involved with the regulation of complement cascade and cell cycle progression and differentiation, respectively. Thus, the absence of age-matched cases and controls probably did not significantly affect our data. Additionally, MSC from both groups were expanded in vitro, which might introduce artifacts. However, unfortunately, MSC from normal donors and patients with MM are found in very low number in the BM. Therefore, in order to obtain the appropriate number of cells to perform experiments, these cells must be expanded in vitro previously. Finally, MM-MSCs were expanded in vitro in the absence of MM tumor cells. However, despite this limitation, the comparisons among different studies, including those that were carried out in co-culture, allow the researchers to  generate quite interesting hypotheses, which can be tested through comparisons between the monoculture and the co-culture of the MM-MSC with or without MM cells.
In summary, our study demonstrated that MM-MSCs have a distinct gene expression profile when compared to the ND-MSCs, corroborating previous studies. The functional enrichment analysis of the gene co-expression network revealed that the main deregulated functions in MM-MSC are related to cell cycle progression, activation of the immune response, and to bone metabolism, which may contribute directly or indirectly to MM physiopathology. Due to the essential role of these cells in the maintenance and progression of MM, potential therapeutic targets and new drugs capable of disrupting the interactions between MM-MSC and MM cells are welcome.

Casuistic
Comparison design