Expression patterns of LILRBs in normal tissues and cancer cell lines
In our first attempt, the expression patterns of LILRBs in different human tissues were determined based on RPKM values using GTEx [39] (http://www.GTExportal.org/home/). We note that the highest expression of LILIBs were observed in spleen, followed by blood and the lung tissue, while in other tissues, they were only weakly expressed (Fig. 1a and Supplementary Figure S1). Importantly, the preferential enrichment of LILRBs in spleen was further validated in the FANTOM5 and HPA (Human protein atlas) dataset (Supplementary Figure S2 and S3). Next, we explored the expression profiles of LILRBs in cancer cell lines from Cancer Cell Line Encyclopedia (CCLE). As shown in Fig. 1b and Supplementary Figure S4, LILRBs showed relatively high expression in cell lines of malignant hematological cell lines, such as acute myeloid leukemia (AML), acute lymphocytic leukemia (ALL), lymphomas, and multiple myeloma (MM). We also analyzed the expression levels of LILRBs based on RNA-seq in 64 different cell lines, as part of the Cell Atlas in the HPA. This analysis confirmed that LILRBs was highly expressed in cell lines of the myeloid and lymphocytic origin (Supplementary Figure S5). Notably, LILRB1 and 2 showed higher expression in monocytes and the THP1 monocyte cell line. Together, these findings indicated a cellular-, tissue-, and disease- specificity of LILRBs expression.
Analysis of LILRB family gene expression levels in tumor and non-tumor tissues
Previous studies revealed that members of the LILRB family were upregulated in a number of cancers [40, 41]. Here, using pan-cancer datasets from the TCGA project, we found that the expression levels of LILRB1-LILRB5 were positively correlated with each other (Fig. 1c), suggesting they may share some common features or functions; of them LILRB5 showed relatively weak correlation with the other four genes. Combining the normal tissue of the GTEx dataset as controls, we then systematically compared LILRBs expression between tumor and adjacent normal tissue across 28 cancer types (9465 tumor and 7831 normal samples). Surprisingly, LILRBs were significantly dysregulated in almost all cancer types (Fig. 1d and Supplementary Figure S6). For LILRB1, LILRB2, and LILRB4, increased expression in tumors was more commonly seen; whereas LILRB3 and LILRB5 were significantly down-regulated in the majority of cancer types (Fig. 1d and e). For LILRB1-LILRB4, the most remarkable difference was observed between AML and normal its normal counterparts (Fig. 1d and Supplementary Figure S6). We also found that LILRB1-LIRB4 were highly expressed in glioblastoma multiforme (GBM),head and neck squamous cell carcinoma (HNSC), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), pancreatic adenocarcinoma (PAAD), and skin cutaneous melanoma (SKCM), whereas they were markedly decreased in adrenocortical carcinoma (ACC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), and thymoma (THYM), as compared with normal controls (Fig. 1d and Supplementary Figure S6).
LILRBs were significantly up‑regulated in AML and dynamically expressed during normal hematopoiesis
Our observations, together with previous findings [18], reflect an AML-specific expression patterns of LILRBs. We thus further explored the differential expression of LILRBs in two independent AML datasets with normal controls available. For these two datasets, only LILRB1 and LILRB3 were consistently up-regulated in AML (Fig. 1f and g), whereas LILRB5 was down-regulated (Supplementary Figure S7a and b). Interestingly, LILRB4, which was previous shown to be increased in AML [18], exhibited no difference in expression in the two datasets (Fig. 1f and g).
We next sought to determine the expression patterns of LILRBs during normal hematopoiesis using a published data (GSE42519) [42]. We found that the transcripts of LILRBs was relatively high in bone marrow (BM) hematopoietic stem cells (HSCs), with a gradual diminishment when HSC differentiated to committed myeloid progenitors. The expression remained low in megakaryocyte-erythroid progenitor (MEP) and then steadily increased during myeloid maturation, reaching the highest level of expression in the mature polymorphonuclear cell (PMN) (Figure S7c). This dynamic regulation of LILRBs transcripts suggests a potential role of LILRBs in HSC biology.
Landscape of genetic alterations of LILRBs in tumors
We also investigated genetic alterations (including mutations, amplifications, and deletions) frequencies of LILRBs across pan-cancers. The average alteration frequencies of five genes were summarized in Fig. 4a and the oncoprint was present in Supplementary Figure S8a. The highest mutation loads of LILRBs was observed in SKCM (31.53%), followed by uterine corpus endometrial carcinoma (UCEC) (14.56%), LUAD (14.49%), and LUSC (13.14%) (Fig. 4a and c). Overall, LILRB1 was the most highly mutated and LILRB3 the least; the most frequent genomic variants were missense mutations for five genes (Fig. 4b and Supplementary Figure S8b). For diffuse large B-Cell lymphoma (DLBC), mesothelioma (MESO), and THYM, copy number alterations (CNAs) were the only genetic events (Fig. 4a). Amplifications were more commonly seen in cancers such as ACC, uterine carcinosarcoma (UCS), and bladder urothelial carcinoma (BLCA), while deletions were mostly found in cancers like brain lower grade glioma (LGG), ovarian serous cystadenocarcinoma (OV), and testicular germ cell tumor (TGCT) (Fig. 4a and d). In AML where LILRBs expression were markedly changed, however, genetic alteration frequencies of these genes were extremely low, suggesting other mechanisms might contribute to the abnormal LILRBs expression in AML (Fig. 4a, c, and d).
LILRBs were hypomethylated in AML
We next asked whether DNA methylation regulates the expression of LILRBs in cancers. We first retrieved the methylome data of LILRBs across 30 cancer types with matched controls through the human disease methylation database Diseasemeth version 2.0 (http://bio-bigdata.hrbmu.edu.cn/diseasemeth/). Surprisingly, we found that LILRB members were significantly hypomethylated in almost all cancer types analyzed as compared to normal samples (Fig. 3a). We then investigated the correlation between methylation and expression levels of LILRBs using GSCALite [31]. We found that methylation level of LILRB1 was negatively associated with its mRNA expression in most cancer types (Fig. 3b). For other LILRB members, methylation and expression levels were also mainly negatively correlated, with only a few positive correlations (Fig. 3b). Importantly, consistent with the increased expression of LILRBs in AML, significantly hypomethylated promoters of LILRBs were observed in both the Diseasemeth (AML, n = 271; normal, n = 10) and GSE63409 dataset (AML, n = 44; normal, n = 30) (Fig. 3c and d). We then studied the correlation between promoter methylation and expression of LILRBs in TCGA AML dataset. Interestingly, expression of LILRB2, LILRB3, and LILRB4 correlated negatively with promoter methylation and the most significant correlation was observed for LILRB4 (Fig. 3b and e). This observation is consistent with a previous report that decitabine (DAC, a demethylating agent) treatment with AML cells remarkably promoted expression of LILRB family members, especially LILRB4 [43]. Collectively, these results suggest that DNA hypomethylation might contribute to the activation of LILRB members in AML. Further, analyzing the relation between methylation and survival revealed that hypomethylation of LILRBs predicted worse survival in most cancers (Fig. 3f). In AML, hypomethylation of LILRB4 was associated with adverse outcome and hypomethylation of LILRB5 with favorable outcome (Fig. 3f).
Adverse prognostic impact of LILRBs in AML
Next, we used Cox regression analyses to explore the association between LILRBs expression and overall survival (OS) in TCGA pan-cancer datasets. Overall, we found the significance and direction of the prognostic significances varied, depending on the cancer types analyzed. For example, increased expression of LILRB family members were generally associated with worse OS in KIRC, LAML, LGG, TGCT, THYM, and uveal melanoma (UVM). While in SKCM and metastatic SKCM, the reverse was observed (Fig. 4a). Interestingly, a previous study has also performed Cox analyses in TCGA data, showing that LILRB1-LILRB4 negatively impacts the survival of AML patients. It is of particular interest to validate the prognostic value of LILRBs using Kaplan-Meier methods in larger patient cohorts of AML. To this end, we collected five independent datasets from GEO; X-tile was used to determine the optimal thresholds for each LILRB members in TCGA and GEO datasets. First, we were able to validate the adverse prognostic impact for LILRB1-LILRB4 in the TCGA cohorts (Fig. 4b), whereas high LILRB5 was associated with favorable outcome (Supplementary Figure S9). Importantly, the prognostic value of LILRB1-LILRB4 also extended to the event-free survival (EFS) endpoint and cytogenetically normal (CN)-AML subsets (Supplementary Figure S10a-c). Furthermore, the adverse prognostic impact of LILRBs was validated in TCGA microarray data (n = 183) (Supplementary Figure S10d) and other five independent cohorts of AML patients (GSE10358, n = 304; GSE37642 [U133A], n = 422; GSE37642 [U133plus2], n = 140; GSE106291, n = 250; GSE71014, n = 104) (Fig. 5a-e), although in some cases only a trend for shorter OS was observed. For subsequent analyses, we will focus on the role of LILRB1-LILRB5 in AML, due to frequent alterations and high prognostic values of LILRBs in this malignancy.
LILRB4 is aberrantly overexpressed in MLL-rearranged AML and may be a target of MLL fusion proteins
We next asked whether LILRBs expression could be associated with specific molecular subtypes in AML. To this end, we examined the expression differences of LILRBs across published transcriptomic subtypes in the Hemap dataset (including AML, pre-B-ALL, DLBCL, and MM) [27]. As expected, all five LILRB members were more highly expressed in monocyte-like AML, while their expressions were relatively weak in the other three malignancies (Fig. 6a). One unanticipated exception to this overall trend was the strong enrichment of LILRB4 in MLL-rearranged AML (Monocyte − like − MLL) and ALL (KMT2A) (Fig. 6a). To confirm this observation, we subsequently analyzed the transcript levels of LILRB4 in 15 leukemia cell lines with or without MLL-rearrangements from the CCLE database. Leukemia cell lines with presence of MLL-fusion genes exhibited markedly higher LILRB4 expression than those lack MLL-fusion genes, whether LILRB4 expression was detected by RNA-seq (Fig. 6b) or Affymetrix microarray (Supplementary Figure S11a). Accordingly, analysis of three large primary patient datasets (BeatAML, TCGA, and GSE13159) revealed consistently highest LILRB4 expression in MLL-rearranged AML as compared to other cytogenetic/clinicopathologic leukemia entities (Fig. 6c, Supplementary Figure S11b and c). To further confirm the relevance of LILRB4 expression in MLL-rearranged AML, we collected four MLL-rearrangement-related gene signatures from MSigDB and computed ssGSEA scores of these signatures for each sample in the TCGA dataset. Then, we compared the ssGSEA scores computed for high LILRB4-expressing samples with those in low LILRB4-expressing samples. We found gene-sets down-regulated in MLL-rearranged AML (MULLIGHAN_MLL_SIGNATURE_1_DN) showed significantly lower ssGSEA scores in LILRB4-high patients than in LILRB4-low patients; whereas for gene-sets up-regulated in MLL-rearranged AML (MULLIGHAN_MLL_SIGNATURE_1_UP), the opposite was seen (Fig. 6d). Also, the ssGSEA scores of two MLL-rearranged-governed signatures (ROSS_AML_WITH_MLL_FUSIONS and VALK_AML_WITH_11Q23_REARRANGED) were significantly up-regulated in high LILRB4 expressers (Supplementary Figure S11d).
It has been shown that target genes of MLL-fusions were often hypomethylated [44, 45], which is consistent with our previous observations. Also, promoters of these genes were often enriched with transcription activation-associated histone marks (H3K79me2, H3K27ac, and H3K4me3) [32]. To determine whether LILRB4 expression could be directly regulated by MLL-fusion gene, we analyzed a published ChIP-seq dataset (GSE79899) of MLL-fusion proteins, H3K79me2, H3K27ac, and H3K4me3 for MV4-11 (MLL-AF4) and THP-1 (MLL-AF9) cell lines. We found a significant enrichment of MLL-N proteins in the promoter regions of LILRB4 gene for both cell lines, while punctuated binding peaks of H3K79me2, H3K27ac, and H3K4me3 were observed in both the promoter and gene body of LILRB4 (Fig. 6e). Importantly, a similar enrichment of the three epigenetic marks was seen in five other ChIP-seq datasets (H3K79me2 from GSE82116 and GSE71779; H3K27ac from GSE89336 and GSE71776; H3K4me3 from GSE61785 and GSE82116) (Fig. 6f). Overall, these results suggest that MLL fusion proteins may be a direct regulator of LILRB4 expression.
LILRB1 expression correlates with distinct genomic alterations in AML
We then examined the associations between LILRBs expression and the clinical and genetic characteristics in the TCGA AML cohort. We found an association between LILRBs expression and the French-American-British (FAB) classification of AML: a higher percentage of myelomonocytic or monocytic morphology (M4/M5 subtypes) and a lower percentage of FAB M2/M3 was observed in patients with high LILRBs expression (Fig. 7a). Moreover, high LILRBs expressers were more likely to be > 60-year-old and less likely to present with favorable cytogenetics (Fig. 7a).
To determine whether LILRB1-LILRB5 correlated with distinct mutational profiles characterized for AML, we identified significantly mutated genes occurred in patients with high and low LILRB1-LILRB4 expression (as stratified by the median expression value of respective genes), using curated mutational data from TCGA. Overall, we found LILRB1 and LILRB5 expression was associated with more mutational events than the other three genes (Fig. 7b). As shown in Fig. 7c, patients with high LILRB1 expression had higher frequency of mutations in U2AF1 (7% vs 1%) and RUNX1 (14% vs 4%), while IDH1 (14% vs 4%) was more frequently mutated in those with low LILRB1 expression. High LILRB5 expression was positively correlated with TP53 mutations and negatively correlated with FLT3 and WT1 mutations. For other three genes, LILRB2 was associated with mutations in IDH1 and STAG2, LILRB3 with WT1, and LILRB4 with RUNX1 (Fig. 7b).
To further explore the association between LILRBs expression and copy number variation (CNV), we performed GISTIC2.0 analysis of TCGA copy number data and assessed copy number alterations between in two patient groups. We focused on LILRB1, as it was consistently dysregulated and showed the greatest mutational events in AML patients. Interestingly, LILRB1-low patients had no somatic copy number alterations (Supplementary Figure S12), whereas LILRB1-high patients had 14 significantly deleted regions and four significantly amplified region (FDR = 0.25) (Fig. 7d). Interestingly, the majority of genes deleted in LILRB1-high patients were involved in inflammatory responses (including cytokines and genes essential for microbial killing and antigen processing and presentation, see Supplementary data 1 for detail). Also, a number of genes belongs to the cadherin (CDH) and protocadherin (PCDH) family, which often exerts tumor-suppressive functions [46], were significantly deleted. In contrast, LILRB-high AML patients had recurrent amplification at loci essential in AML pathogenesis, including KMT2A and ERG [47, 48] (Fig. 7d).
Correlations Between LILRBs and Tumor Immune Infiltrating Cells (TIICs) in AML
Considering that LILRBs might play important roles in the TME, we further explored the correlations between LILRBs and the level of immune cell infiltration in TCGA AML cohort. It is noteworthy that, among the 22 cell types, monocytes had the highest positive correlations with LILRB1-LILRB4 (Fig. 8a), consistent with previous finding that LILRBs were preferentially expressed in monocytic AML [19, 49]. This monocytic preference was also confirmed in two recently published scRNA-seq datasets of AML (Van Galen AML scRNA, Fig. 8b and FIMM AML scRNA, Supplementary Figure S13a). Interestingly, LILRB4 was exclusively correlated with M2 macrophages (Fig. 8a), a high immunosuppressive component in the TME. By contrast, LILRB1-LIRB4 were negatively correlated with the infiltrating levels of tumor-suppressive immune cells, such as resting T cells CD4 memory cell, CD8 T cells, memory B cells, plasma cells, and resting NK cells (Fig. 8a). Similar results were found by analyzing the CIBERSORT estimates in the GSE10358 and GSE6891 dataset (Supplementary Figure S13b and c). Importantly, when other methods were used for calculating the relative fractions of TIICs, positive associations between LIRB1-LILRB4 and monocytes were consistently seen, while negative associations between LILRB1-LILRB4 and CD8 T cells were proved for most-if not all-methods in all three datasets (Supplementary Figure S13d-f). Further analysis of normal cell populations from the Hemap dataset revealed that LILRBs were highly expressed in myeloid lineage immune cells (monocytes, macrophages, dendric cells, myeloid progenitors, and neutrophils), with consistent low expression in T cells (CD4 + T cells and T/NK cells) (Fig. 8c). Collectively, these findings further confirmed the immunosuppressive roles of LILRBs in cancer TME.
Correlation between LILRBs and immune checkpoints in AML
Given that immune checkpoints have been proved to be promising therapeutic target for cancer treatment, we therefore evaluated the relationship between LILRBs and a collection of checkpoint genes describe by De Simone et al. [37]. Results from Spearman correlation analyses are given in Supplementary Data 2. As shown in the Circos plots, LILRB1-LILRB3 all showed strong positive correlations with CD86, VISTA, and HAVCR2 (Fig. 8d-f), indicating a possible synergistic effect between these genes. In contrast, no significant correlations were observed between LILRB4/5 and these checkpoints (Supplementary Figure S14a and b). These results further highlight LILRBs potentially as major signaling pathways involved in immunosuppression in the AML microenvironment.
The biological significance of LILRBs expression in AML
We then sought to investigate the biological features associated with LILRBs in AML. Since the expression of five LILRB members were highly correlated, a comparison of gene expression profiles of patients with high and low LILRB1 expression (as determined by the median expression value) was performed. Overall, 799 genes (490 up- and 309 downregulated; adjusted p < 0.05; log2FC ≤ − 1.5 or log2FC ≥ 1.5) were differentially expressed in LILRB1high versus LILRB1low patients (Fig. 9a and Supplementary Data 3). Among the genes positively correlated with LILRB1 were, as expected, the other members of the LILRB family (Fig. 9a). Also, genes associated with presence of monocytes/macrophages (CD14, CD68) or M2 macrophage polarization (MSR1, MRC1, CD163) were significantly up-regulated in high LILRB1 expressers (Fig. 9a), in line with our previous findings. Next, we used STRING database to construct a protein-protein interaction (PPI) network of the differentially expressed genes (DEGs), with a confidence score > 0.90. Genes interacted with LILRB1 and their sub-networks were shown through Cytoscape software (Fig. 9b). We found 12 genes directly interacting with LILRB1: PILRA, TLR8, SIGLEC7, CD300C, FCGR2A, FCGR2B, FCGR3A, CD86, FGR, HCK, IL10, ITGAX. Among them, CD300C, FCGR2A, FCGR2B, and FCGR3A also had connections with the other four LILRB members (Fig. 9b). GeneMANIA results also revealed that genes of the FCGR and CD300 family were closely correlated with LILRBs. These genes were mainly involved in negative regulation of leukocyte mediated immunity and negative regulation of immune system process (Supplementary Figure S15a).
We then performed GO analysis using these DEGs and the top 10 significant terms of BP, MF and CC enrichment analysis were shown (Fig. 9c). Notably, in terms of BP, immune response-related processes were significantly enriched, such as inflammatory response, immune system process, and immune response. KEGG and Reactome Pathway analyses also revealed immune response pathways, including cytokine − cytokine receptor interaction, cytokine signaling in immune system, innate immune system, antigen processing − cross presentation, and adaptive immune system were mainly enriched (Fig. 11d and Supplementary Figure S15b).
Finally, GSEA was conducted in the LILRB1high and LILRB1low cohorts. For the C2 collection of curated gene sets from the MSigDB, the VALK_AML_CLUSTER_5 gene set (96% of the samples are FAB M4 or M5 subtype) was predominantly enriched in LILRB1high group. Also enriched were gene sets of MLL-fusion and NPM1-mutation, two distinct entities often associated with monocytic features of AML (Fig. 10a). For the C7 immunologic collection, the LILRB1high group had principal enrichment in genes up-regulated in monocytes compared to other immune cells (Fig. 10b), and multiple immune activities were enriched in the LILRB1high group for HALLMARK gene sets (Fig. 10c).