Aberrant metabolic processes promote the immunosuppressive microenvironment in multiple myeloma

Introduction Multiple myeloma (MM) is still an incurable plasma cell malignancy. The efficacy of immunotherapy on MM remains unsatisfactory, and the underlying molecular mechanisms still are not fully understood. Methods In this study, we delineated the dynamic features of immune cell in MM bone marrow (BM) along with elevated tumor cell infiltration by single-cell RNA sequencing (scRNA-seq), and investigated the underlying mechanisms on dysfunction of immune cells associated with myelomagenesis. Results We found that immune cells were activated in those patients with low infiltration of tumor cells, meanwhile suppressed with elevated infiltration of MM cells, which facilitated MM escaping from immune surveillance. Besides PD-1, abnormal expression of PIM kinases, KLRB1 and KLRC1 were involved in the defect of immune cells in MM patients. Importantly, we found aberrant metabolic processes were associated with the immunosuppressive microenvironment in MM patients. Disordered amino acid metabolism promoted the dysfunction of cytotoxicity CD8 T cells as well as lipid metabolism disorder was associated with the dysregulation of NK and DCs in MM. As metabolic checkpoints, PIM kinases would be potential effective strategies for MM immunotherapy. Discussion In summary, redressing the disordered metabolism should be the key points to get promising effects in immune-based therapies.

cells was then assessed based on four metrics step by step: (1) the number of detected genes per cell; (2) the number of detected UMI per cell; (3) the proportion of mitochondrial gene counts; and (4) the proportion of rRNA genes counts (RNA18S5 or RNA28S5). The following criteria were then applied to filter low-quality cells: gene number < 200 or > 6,000, UMI > 1000, ribosomal gene proportion > 0.4 or mitochondrial gene proportion > 0.3. Finally, a total of 42,936 cells were incorporated into further analysis. For the integration of the cells from different samples, the Gene-cell matrix of all samples was integrated with Seurat to remove batch effects across different samples. In parameter settings, the first 30 dimensions of canonical correlation analysis (CCA) and principal component analysis (PCA) were used.

Dimensionality reduction, clustering of cells, and visualization
The filtered gene-cell matrix was first normalized using "LogNormalize" methods in Seurat v.3 with default parameters. The top 2,000 variable genes were then identified using the "vst" method in the Seurat Find Variable Features function. PCA was performed using the top 2,000 variable genes. Graph-based clustering was performed on the PCA-reduced data for clustering analysis with Seurat v.3. The resolution was set to 0.5 to obtain a more refined result. Briefly, the first 50 PCs of the integrated gene-cell matrix were used to construct a shared nearest-neighbor graph (SNN; FindNeighbors() in Seurat), and this SNN was used to cluster the dataset (FindClusters()) using a graph-based modularityoptimization algorithm of the Louvain method for community detection. Then UMAP was performed on the top 30 principal components for visualizing the cells.
Cell cluster annotation with specific maker genes expression Find All Markers in Seurat (Wilcoxon rank-sum test) was used to perform differential gene expression analysis. For each cluster, marker genes were generated relative to all other cells. Cellular identity was determined by comparing cluster-specific markers of each cluster to known cell-type-specific genes from previous studies. Cluster annotation was confirmed using the R package SingleR, which compares the transcriptome of every single cell to reference datasets to determine cellular identity. DEGs identification and functional enrichment analysis Differential expression genes (DEGs) among different sample groups within a cluster were identified using FindMarkers in Seurat (wilcox.test), with default parameters. A gene was considered significantly differentially expressed if the false discovery rate (FDR) < 0.05 and expression fold change (FC) > 1.3. The heat map was then generated using the pheatmap R package for filtered DEGs. The gene set variation analysis (GSVA) was applied to the scRNA-seq data, and average GSVA scores were calculated for each cell using the GSVA function in the GSVA software package. Differential pathway analysis between clusters was done with the limma R software package. Significantly enriched pathways were identified with an FDR value < 0.05. Gene ontology (GO) enrichment analysis on DEGs was performed using cluster profiler4.

Cell function analysis based on scRNA-seq
The cytotoxic score and exhausted score for T cells and active score for dendritic cells (DCs) were defined by AddModuleScore. CellPhoneDB were used to analyze cell-cell interactions among immune cells. The interaction strength refers to the total average of the mean expression value of a single ligand-receptor partner in the corresponding interacting cell type. The expression of any complex output by CellPhoneDB was calculated as the sum of the expression of the component genes. Mouse model and flow cytometry analysis C57BL/KaLwRij mice (purchased from Harlan Laboratories Inc., Netherlands) and housed in our lab were utilized in the present study, according to the protocol reported by our previous study. 1×10 6 5TGM1-GFP cells were injected into C57BL/KaLwRij mouse (female) via tail vein. Control mice received the PBS injection of equal volume. Bone marrow cells were collected 5 weeks after MM cell injection, and FACS was performed to analyze the composition in bone marrow cells. Fresh bone marrow aspirates obtained from MM patients after informed consent were placed in ethylenediaminetetraacetic acid (EDTA)-containing tubes and immediately transported to the lab. Fresh BMMCs were isolated by Ficoll density-gradient centrifugation and stained with fluorochrome-conjugated antibodies for 15 minutes at room temperature. Flow cytometry was performed on CantoⅡ flow cytom eter (BD Biosciences), and the data were analyzed by Flowjo V10 software (Treestar). The detailed information with the antibodies utilized is listed in suppl. Table 1.

Supplemental Figure 3. NK cell sub-clusters in MM patients
A. Heatmap shows the expression profile of top 10 signature genes for the definition of NK/NKT sub-clusters. The top bars label the cluster numbers corresponding to the subcluster number in Fig.4A. B. Bar charts show the proportions of NK/NKT cell sub-clusters from each HD and MM patients. The sub-cluster numbers in right correspond to the ones in Fig. 4A. C. GO Enrichment of DEGs in NK-S100A8 between high and low tumor burden groups of MM patients. Each dot in the graphs represents a single gene from DEGs. Upregulated genes are indicated as red dots and downregulated genes are indicated as blue dots. The color bar indicates the z-score of each pathway.

Supplemental Figure 4. Myeloid sub-clusters in MM patients
A. Heatmap shows the expression profile of top 10 signature genes for the definition of myeloid cell sub-clusters. The top bars label the cluster numbers corresponding to the sub-cluster number in Fig.5A. B. Heatmap shows the DEGs in macro-IL1B among HD and MM patients in different infiltration groups (HD: n=6; Low: n=6; High: n=6). C. Heatmap shows the DEGs in mono-FCGR3A among HD and MM patients in different infiltration groups (HD: n=6; Low: n=6; High: n=6).