Pan-cancer analysis reveals distinct clinical, genomic, and immunological features of the LILRB immune checkpoint family in acute myeloid leukemia.

Leukocyte immunoglobulin (Ig)-like receptor Bs (LILRBs), a family of type I transmembrane glycoproteins, are known to inhibit immune activation. Here, we comprehensively evaluated the molecular, prognostic, and immunological characteristics of LILRB members in a broad spectrum of cancer types, focusing on their roles in acute myeloid leukemia (AML). We showed that LILRBs were significantly dysregulated in a number of cancers and were associated with immune-inhibitory phenotypes. Clinically, high expression of LILRB1-LILRB4 predicted poor survival in six independent AML cohorts. Genetically, LILRB1 was associated with more mutational events than other LILRB members, and multiple genes involved in immune activation were deleted in LILRB1 high patients. Epigenetically, LILRB4 was significantly hypomethylated and marked by MLL-associated histone modifications in AML. Immunologically, LILRBs were positively associated with monocytic cells, including M2 macrophages, but were negatively associated with tumor-suppressive CD8 T cells. Importantly, patients with higher LILRB expression generally showed a better response to immune checkpoint blockade (ICB) in five independent immunotherapy cohorts. Our findings reveal critical immunological and clinical implications of LILRBs in AML and indicate that LILRBs may represent promising targets for immunotherapy of AML.

Introduction

Acute myeloid leukemia (AML) is a highly fatal hematopoietic malignancy marked by various cytogenetic and molecular abnormalities and variable responses to treatment.1, 2, 3 Currently, the mainstay of treatment for AML is cytotoxic chemotherapy, yet chemoresistance and relapse are commonly seen in clinical practice. Some novel regimens, such as hypomethylating agents (HMAs), Bcl-2 inhibitors, Fms-like tyrosine kinase 3 (FLT3), and isocitrate dehydrogenase (IDH) inhibition, have shown promising results in certain subsets of patients with AML., Since 2017, the US Food and Drug Administration (FDA) has approved some new agents, such as enasidenib for patients with relapsed/refractory IDH2-mutated AML, gilteritinib for FLT3-mutated AML, and liposomal cytarabine-daunorubicin CPX-351 for therapy- and myelodysplastic syndrome (MDS)-related AML. However, there still remains an urgent need to develop novel effective therapies for various subsets of AML.

Of note, immune checkpoint inhibitors (e.g., anti-PD-1 and anti-PD-L1 antibodies) have revolutionized cancer treatment during the past decade in treating cancers such as non-small cell lung carcinoma and melanoma; however, the transfer of immunotherapy to AML has been less successful than to other cancers. Indeed, the AML microenvironment is predominantly immunosuppressive. For example, we have previously demonstrated that M2 macrophages, a classical immunosuppressive component, were preferentially enriched in AML than other hematological malignancies and normal controls. Also, a recent single-cell RNA sequencing (RNA-seq) study has identified proportionally fewer T cells and cytotoxic T lymphocytes (CTLs) in AML than normal controls, and the function of these T cells is profoundly impaired, probably mediated by CD14+ monocyte-like cells., Moreover, Noviello et al. found that bone marrow (BM) T cells at AML relapse showed an exhausted phenotype, which was absent in patients maintaining long-term complete response. These findings suggest encouraging therapeutic opportunities by modulating the immune environment in AML.

Leukocyte immunoglobulin (Ig)-like receptor subfamily B (LILRB) proteins are a group of type I transmembrane glycoproteins with extracellular Ig-like domains that bind ligands and intracellular immunoreceptor tyrosine-based inhibitory motifs (ITIMs). This group of receptors contains 5 members (LILRB1–LILRB5) mainly expressed in hematopoietic cells and also various types of tumors. As these proteins negatively regulate immune activation,15, 16, 17 they are often considered as immunosuppressive components in the tumor microenvironment (TME). In AML, the TME-modulating role of LILRBs has recently come into focus, especially for LILRB4. Gui et al. demonstrated that LILRB4 facilitates tissue infiltration of AML cells by substantially suppressing T cell activities, while blocking LILRB4 activity efficiently inhibited AML development in vitro and in vivo., In addition, LILRB1 was found to be more highly expressed in dysfunctional CD8+ T cells from AML than T cells from healthy controls. Interestingly, a non-immunological AML-promoting role was reported for LILRB2, which binds ANGPLT2 to maintain stemness of normal stem cells and support leukemia development by inhibiting differentiation of AML cells. Despite the functional importance of LILRBs in AML, there lacked a systematic study to explore the expression patterns, clinical implications, and immunological features of all LILRB members in AML. Therefore, in this study, drawing on rich multi-omics data in the public domain, we comprehensively evaluated the transcriptional levels and prognostic significances of LILRB members in a broad spectrum of cancer types, focusing on its role in AML. In addition, we systematically characterized the genomic and immune landscape in patients with AML with altered LILRB expression.

Results

Landscape of genetic and expression alterations of LILRBs across cancer types

We first determined the expression patterns of LILRBs in different human tissues based on reads per kilobase of transcript per million mapped reads (RPKM) values using Genotype-Tissue Expression (GTEx; http://www.GTExportal.org/home/). We observe that LILIBs were highest expressed in the spleen, followed by blood and the lung tissue, while weakly expressed in other tissues (Figure S1). Importantly, the preferential enrichment of LILRBs in spleen was further validated in the FANTOM5 and Human Protein Atlas (HPA) dataset (Figures S2 and S3). Next, using Cancer Cell Line Encyclopedia (CCLE), we showed that LILRBs were relatively highly expressed in malignant hematological cell lines from AML, acute lymphocytic leukemia (ALL), lymphomas, and multiple myeloma (MM) (Figure S4). Moreover, we observed a strong protein expression of LILRB1–LILRB4 in monocytes via Human Proteome Map (https://www.humanproteomemap.org/) (Figure S5). Together, these findings indicated cellular-, tissue-, and disease-specific LILRB expression. Combining the normal tissue of the GTEx dataset as controls, we then systematically compared LILRB expression between tumor and adjacent normal tissue across 28 cancer types (9,465 tumor and 7,831 normal samples). Surprisingly, LILRBs were significantly dysregulated in almost all cancer types (Figures 1A and 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 (Figures 1A and 1B). For LILRB1–LILRB4, the most remarkable difference was observed between AML and normal counterparts (Figures 1A and S6). We also investigated genetic alteration (including mutations, amplifications, and deletions) frequencies of LILRBs across pan-cancers. The average alteration frequencies of five genes are summarized in Figure S7, and the oncoprint is present in Figure S8A. The highest mutation loads of LILRBs were observed in skin cutaneous melanoma (SKCM) (Figure 1C). Overall, LILRB1 was the most highly mutated and LILRB3 the least; the most frequent genomic variants were missense mutations for five genes (Figures 1D and S8B). Amplifications were more commonly seen in cancers such as adrenocortical carcinoma (ACC) and uterine carcinosarcoma (UCS), while deletions were mostly found in brain lower grade glioma (LGG) (Figures 1E and S7). By analyzing the methylome data of LILRBs across 30 cancer types with matched controls through the human disease methylation database Diseasemeth v.2.0 (http://bio-bigdata.hrbmu.edu.cn/diseasemeth/), we found that LILRB members were significantly hypomethylated in almost all cancer types analyzed compared with normal samples (Figure 1F). Furthermore, the level of methylation was negatively associated with the level of mRNA expression in most cancer types (Figure S9A). Analyzing the relation between methylation and survival revealed that hypomethylation of LILRBs predicted worse survival in most cancers (Figure S9B). Finally, Cox regression analyses were used to explore the association between LILRB expression and overall survival (OS) in TCGA pan-cancer datasets. Overall, we found that the significance and direction of the prognostic significances varied, depending on the cancer types analyzed. For example, increased expression of LILRB family members was generally associated with worse OS in kidney renal clear cell carcinoma (KIRC), LAML, LGG, thymoma (THYM), and uveal melanoma (UVM), while in SKCM, the reverse was observed (Figure 1G).

Figure 1

Landscape of genetic and expression alterations of LILRBs across cancer types

(A) Heatmap of differential expression profiles of LILRBs between tumor and normal tissue samples, combining data from TCGA and GTEx databases. The color depicts the log2-transformed fold change (Log2FC) between tumor and normal tissues. (B) Bar plot showing genes significantly up- and down-regulated (p < 0.05) across different cancer types. Red, up-regulated expression; blue, down-regulated expression. (C) Heatmap showing mutation frequencies of LILRBs across different cancer types. Numbers on the cells represent mutation percentages. (D) Bar plot showing the percentages of various mutation types for five LILRB genes. (E) Pie plots showing the percentages of various copy-number alteration (CNA) types for five LILRB genes. (F) Heatmap of differential methylation profiles of LILRBs between tumor and normal samples, using data from the Diseasemeth database. The color depicts methylation differences between tumor (T) and normal (N) tissues. (G) Association between LILRB expression and patient prognosis across 33 cancer types as determined by the Cox regression model. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001. See also Figures S1–S9.

Association between LILRB expression and immune responses in cancers

LILRB family genes have been known for their immune inhibitory functions in cancers. For example, LILRB4 has been shown to suppress T cell activation and support tissue infiltration of AML cells. We hypothesized that they might be associated with specific immunologic programs in cancers. Previously, Thorsson et al. have identified six immune subtypes across cancers: C1 (wound healing), C2 (INF-r dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (tumor growth factor beta [TGF-β] dominant). We found that all LILRB members exhibit the highest expression in C6 (Figure 2A), a highly immunosuppressive subtype displaying increased M2 macrophage infiltration, and conferred the worst prognosis in respective tumors. Recent research has developed four distinct TME subtypes conserved across a broad array of cancers: (1) immune enriched, fibrotic (IE/F), (2) immune enriched, non-fibrotic (IE), (3) fibrotic (F), and (4) immune depleted (D). We found that patients with high LILRB expression possessed primarily subtypes IE/F and IE, whereas patients with low LILRB expression were mainly concentrated in the D subtype (Figure 2B). Next, we investigated the association between LILRBs and 29 TME signature scores calculated using TCGA pan-cancer data. We observed strong positive correlations between LILRB expression with both anti-tumor and tumor-promoting immune processes (especially checkpoint inhibition) but weak correlations with stromal components and cancer cell properties (Figure 2C).

Figure 2

Association of LILRB family gene expression with tumor microenvironment factors

(A) The expressions of LILRBs within different immune infiltrate subtypes across pan-cancers tested with ANOVA. C1: wound healing, C2: INF-r dominant, C3: inflammatory, C4: lymphocyte depleted, C5: immunologically quiet, and C6: TGF-β-dominant.(B) Sankey diagram showing the links between LILRB expression and TME subtypes across pan-cancers. IE/F, immune-enriched, fibrotic; IE, immune-enriched, non-fibrotic; F, fibrotic; D, immune-depleted. LILRB expression statuses were stratified by the median expression value of respective genes. (C) Correlation matrix plots showing the association between expressions of LILRBs and 29 TME signature scores calculated using TCGA pan-cancer data. Spearman correlation was used for testing.

Validation of the prognostic significance of LILRBs in AML

Our data, along with previous studies, reflect an AML-specific expression pattern of LILRBs., In this study, we focus on LILRBs in AML. Cox analyses in TCGA data showed that LILRB1–LILRB4 negatively impact the survival of patients with AML. 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 member in TCGA and GEO datasets. First, we were able to validate the adverse prognostic impact for LILRB1–LILRB4 in the TCGA cohorts, whereas high LILRB5 was associated with a favorable outcome (Figure S10). Importantly, the prognostic value of LILRB1–LILRB4 also extended to the event-free survival (EFS) endpoint and cytogenetically normal (CN) AML subsets (Figures S11A–S11C). Furthermore, the adverse prognostic impact of LILRBs was validated in TCGA microarray data (n = 183) (Figure S11D) and five other independent cohorts of patients with AML (GEO: GSE10358, n = 304; GSE37642 [U133A], n = 422; GSE37642 [U133plus2], n = 140; GSE106291, n = 250; GSE71014, n = 104) (Figures 3A–3E and S12A), although in some cases, only a trend for shorter OS was observed. However, LILRB5 showed opposite prognostic effects in the GEO: GSE37642 (U133A) and GSE71014 datasets compared with that of TCGA, and no statistically significant associations were detected in the other three datasets (Figures S12A and S12B).

Figure 3

Independent validation of the prognostic significance of LILRBs in five GEO datasets

(A–E) Kaplan-Meier curves representing OS of five AML cohorts from GEO (GEO: GSE10358, n = 304; GSE37642 [U133A], n = 422; GSE37642 [U133plus2], n = 140; GSE106291, n = 250; GSE71014, n = 104) based on the expression of indicated LILRB members (LILRB1–LILRB4). The optimal cutoff of each gene was determined by the X-tile method. See also Figures S10–S12.

LILRB expression correlates with distinct genomic alterations in AML

We then examined the associations between LILRB expression and the clinical and genetic characteristics in the TCGA AML cohort. We found an association between LILRB 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 were observed in patients with high LILRB expression (Figure 4A). Moreover, high LILRB expressers were more likely to be >60 years old and less likely to present with favorable cytogenetics (Figure 4A).

Figure 4

LILRB expression correlates with distinct genomic alterations in AML

(A) Heatmap showing association between LILRB expression and clinical characteristics in the TCGA AML cohort (top panel). Pie charts showing the chi-squared test of clinicopathologic factors for LILRB1 status using the median expression as the cut off (bottom panel). (B) Bubble plot showing associations between the expression of LILRB1–LILRB5 (as stratified by the median expression value of respective genes) and common mutational events in the TCGA dataset (top panel). Bubble size indicates -log10 (Fisher test p value). Color signifies log10 (odds ratio [OR]), positive association is indicated with red circles, negative with blue circles, and non-association with gray circles. LILRB1 expression stratified by RUNX1 mutation status, and LILRB5 expression stratified by TP53 mutation status (bottom panel). (C) Co-bar plots showing the comparison of mutational profiles between patients with high and low LILRB1 expression (as stratified by the median expression value) in the TCGA dataset. (D) GISTIC analyses identified recurrent CNAs in patients with AML with high LILRB1 expression. (E) The CNV variation frequency of indicated oncogenes and tumor-suppressor genes in patients with AML with high LILRB1 expression. The height of the column represented the alteration frequency. Yellow dot indicates the amplification frequency; blue dot indicates the deletion frequency. See also Figure S13.

We hypothesized that altered LILRB expression would have an impact on the mutation landscape of AML patients. To determine whether LILRB1–LILRB5 correlated with distinct mutational profiles characterized for AML, we identified significantly mutated genes that 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 (Figure 4B). As shown in Figure 4C, patients with high LILRB1 expression had a higher frequency of mutations in U2AF1 (7% versus 1%) and RUNX1 (14% versus 4%), while IDH1 (14% versus 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 (Figure S13A). In addition, RUNX1-mutated AML highly expressed the LILRB1 gene, and TP53 mutations were linked to higher LILRB5 expression (Figure 4B). For the other three genes, LILRB2 was associated with mutations in IDH1 and STAG2, LILRB3 with WT1, and LILRB4 with RUNX1 (Figure 4B).

We also considered the possibility that specific regions of the genome may be preferentially focally amplified or deleted in patients with high or low LILRB expression. We therefore performed GISTIC2.0 analysis of TCGA copy-number data and assessed copy-number variations (CNVs) in two patient groups. We focused on LILRB1, as it was consistently dysregulated and showed the greatest mutational events in patients with AML. Interestingly, LILRB1low patients had no somatic copy-number alterations (Figure S13B), whereas LILRB1high patients had 14 significantly deleted regions and four significantly amplified regions (false discovery rate [FDR] = 0.25) (Figure 4D). Interestingly, the majority of genes deleted in LILRB1high patients were involved in inflammatory responses (including cytokines and genes essential for microbial killing and antigen processing and presentation; see Table S1 for details). Also, a number of genes belong to the cadherin (CDH), protocadherin (PCDH) family (e.g., CDH1, PCDH1, and PCDH12), and cyclin-dependent kinase (CDK) inhibitors (e.g., CDKN1B, CDKN2A, and CDKN2B), which often exert tumor-suppressive functions,, were significantly deleted. In contrast, LILRB-high patients with AML had recurrent amplification at loci essential in AML pathogenesis, including KMT2A and ERG (Figure 4E).,

LILRB4 is aberrantly overexpressed in MLL-rearranged AML and may be a target of MLL fusion proteins

We next asked whether LILRB 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, diffuse large B cell lymphoma [DLBCL], and MM). 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 (Figure 5A). One exception to this overall trend was the strong enrichment of LILRB4 in MLL-rearranged AML (monocyte-like MLL) and ALL (KMT2A) (Figure 5A). This agreed favorably with previous findings that LILRB4 was correlated with MLL-rearranged leukemia., To further 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 the presence of MLL fusion genes exhibited markedly higher LILRB4 expression than those lack MLL fusion genes, whether LILRB4 expression was detected by RNA-seq (Figure 4B) or Affymetrix microarray (Figure S14A). Accordingly, analysis of three large primary patient datasets (BeatAML, TCGA, and GSE13159) revealed consistently higher LILRB4 expression in MLL-rearranged AML compared with other cytogenetic/clinicopathologic leukemia entities (Figures 5C, S14B, and S14C). 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 (Figure 5D). 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 (Figure S14D).

Figure 5

LILRB4 is aberrantly overexpressed in MLL-rearranged AML and is likely a direct target of MLL fusion proteins

(A) Expression differences of LILRB genes in molecular subtypes of AML and pre-B-ALL, DLBCL, and MM. The expression FC between each subtype and the remaining samples in the same disease were compared using the Wilcoxon rank-sum test. The color of the dots indicates FCs (log2), and size indicates the FDR values. The FDR values were categorized into six groups based on significance cutoffs for visualization (0.05, 0.01, 0.001, 1 × 10-5, 1 × 10-16). (B) Bar plot showing LIRB4 expression (RNA-seq) in non-MLL-rearranged (HEL, MEG01, KASUMI1, KG1, NB4, K562, HL60, U937) and MLL-rearranged (SEM, MONOMAC6, NOMO1, RS411, MOLM13, THP1, MV411) cell lines from the CCLE database. The dotted line represents the mean expression of LILRB4 across all cell lines analyzed. (C) Comparison of LILRB4 expression among human primary AML cases with MLL rearrangements and those without MLL rearrangements in the BeatAML dataset. (D) Box plots showing ssGSEA scores of two MLL-related gene signatures (MULLIGHAN_MLL_SIGNATURE_1_DN and MULLIGHAN_MLL_SIGNATURE_1_UP) between patients (TCGA dataset) with high and low LILRB4 expression (as stratified by the median expression value). (E) ChIP-seq tracks for MLL fusion proteins, H3K79me2, H3K27ac, and H3K4me3 at LILRB4 gene loci in MV4-11 and THP-1 cell lines. ChIP-seq data were obtained from GEO: GSE79899. (F) ChIP-seq tracks for H3K79me2, H3K27ac, and H3K4me3 at LILRB4 gene loci in MV4-11- and MLL-AF9-transformed blast cells. ChIP-seq data were obtained from GEO: GSE82116, GSE71779, GSE89336, GSE71776, and GSE61785. See also Figures S14 and S15.

It has been shown that target genes of MLL fusions were often hypomethylated., Consistently, 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) (Figures S15A and S15B). Moreover, the expression of LILRB2, LILRB3, and LILRB4 correlated negatively with promoter methylation, and the most significant correlation was observed for LILRB4 (Figure S15C). This observation is consistent with a previous report that decitabine (DAC; a demethylating agent) treatment with AML cells remarkably promoted the expression of LILRB family members, especially LILRB4. Also, promoters of MLL fusion target genes were often enriched with transcription activation-associated histone markers (H3K79me2, H3K27ac, and H3K4me3). To determine whether LILRB4 expression could be directly regulated by the MLL fusion gene, we analyzed a published chromatin immunoprecipitation (ChIP)-seq dataset (GEO: 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 (Figure 5E). Importantly, a similar enrichment of the three epigenetic marks was seen in five other ChIP-seq datasets (H3K79me2 from GEO: GSE82116 and GSE71779; H3K27ac from GEO: GSE89336 and GSE71776; H3K4me3 from GEO: GSE61785 and GSE82116) (Figure 5F). Overall, these results suggest that MLL fusion proteins may be a direct regulator of LILRB4 expression.

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 the TCGA AML cohort. It is noteworthy that, among the 22 cell types, monocytes had the highest positive correlations with LILRB1–LILRB4, while only a weak correlation was observed between LILRB5 and monocytes (Figure 6A), consistent with the previous findings that LILRBs were preferentially expressed in monocytic AML., This monocytic preference was also confirmed in two recently published single-cell RNA-seq (scRNA-seq) datasets of AML (Van Galen AML scRNA, Figure 6B, and FIMM AML scRNA, Figure S16A). Interestingly, LILRB4 was exclusively correlated with M2 macrophages (Figure 6A), 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 memory CD4 T cells, CD8 T cells, memory B cells, plasma cells, and resting natural killer (NK) cells (Figure 6A). Similar results were found by analyzing the CIBERSORT estimates in the GEO: GSE10358 and GSE6891 datasets (Figures S16B and S16C). 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 (Figures S16D–S16F). 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) (Figure 6C). Collectively, these findings further confirmed the immunosuppressive roles of LILRBs in cancer TME.

Figure 6

The relation between LILRBs expression with immune cell infiltration and immune checkpoints

(A) Correlation matrix plot showing correlations between LILRBs and tumor immune infiltrating cells (TIICs). The overall immune cell compositions were estimated by CIBERSORT in the TCGA dataset. (B) Dot plot showing expression pattens of LILRBs and selected immune checkpoint genes in annotated cell types from 16 AML scRNA-seq samples (Van Galen AML scRNA). The color of the dots indicates average expression, and size indicates percentage of cells with detectable expression. (C) Heatmap showing LILRB expression in normal cell populations from the Hemap dataset. (D) Correlogram showing correlations between the expression of LILRBs and selected immune checkpoint genes in the TCGA dataset. Positive correlation was marked with orange and negative correlation with blue. See also Figure S16.

Correlation between LILRBs and immune checkpoints in AML

Given that immune checkpoints have been proven to be promising therapeutic targets for cancer treatment, we therefore evaluated the relationship between LILRBs and a collection of checkpoint genes describe by De Simone et al. Results from Spearman correlation analyses are given in Table S2. As shown in the correlogram, LILRB1–LILRB3 all showed strong positive correlations with CD48, CD86, PD-L2, TIM-3, and VISTA (Figure 6D), while relatively weaker associations were observed between LILRB4/5 and these checkpoints. Moreover, analysis at the single-cell level revealed that CD86 and VISTA, which, like LILRBs, were preferentially expressed in monocytes (Figures 6B and S16A). In contrast, LILRBs did not show any correlations with PD-1, and only weak correlations between CTLA-4 and LILRBs (except for LILRB4) were observed (Figure 6D). These results further highlight LILRBs potentially as major signaling pathways involved in immunosuppression in the AML microenvironment.

LILRB expression predicts responses to immunotherapy

Considering the strong connection between LILRB expression and immune response, we next asked whether LILRB expression can be utilized as a tool to predict response to immune checkpoint blockade (ICB). We first used Tumor Immune Dysfunction and Exclusion (TIDE; http://tide.dfci.harvard.edu/) to assess the potential of LILRBs as new biomarkers by comparing their predictive power with that of existing biomarkers. Surprisingly, we found that LILRBs had an area under the curve (AUC) of >0.5 in 17 of the 22 (77%) ICB subcohorts, comparable to the predictive performance of TIDE (18 out of 25, 72%). It also showed a higher predictive value than tumor mutational burden (TMB), T clonality, and B clonality (Figure S17). Moreover, our analyses revealed that LILRBs could predict patients' survival in five independent immunotherapy cohorts, including two melanoma cohorts treated with anti-PD-1 therapy (Liu2019_PD1_Melanoma and Gide2019_PD1_Melanoma), two melanoma cohorts treated with anti-CTLA-4 therapy (Nathanson2017_CTLA4_Melanoma_Post and VanAllen2015_CTLA4_Melanoma), and one clear cell renal cell carcinoma (ccRCC) cohort with anti-PD-1 monotherapy (Miao2018_ICB_Kidney_Clear). Remarkably, LILRBs exhibited the highest predictive value in three of the five datasets (Figure 7A). In the Gide2019_PD1_Melanoma and Nathanson2017_CTLA4_Melanoma_Post cohorts, the percentage of responders (complete response [CR] or partial response [PR]) to ICB was generally higher in patients with high LILRB expression than those with low LILRB expression (Figure 7B). Similar findings could be extended to other cancer types (lung cancer and gastric cancer) with ICB treatment (Jung2019_PD1/PDL1_ Lung and Kim2018_PD1_Gastric) and melanoma patients treated with adoptive T cell therapy (ACT) (Lauss2017_ACT_Melanoma) (Figure S18). In addition, high LILRB expressions were generally correlated with PD-1/CTLA-4 up-regulation in cohorts treated with respective antibodies (Figure S19A). To test the potential of LILRBs in predicting ICB response in patients with AML, we checked the relationship of LILRBs with expression signatures for predicting ICB response in the TCGA AML dataset. Surprisingly, we found a negative correlation of LILRBs with T cell-exclusion signatures, including myeloid-derived suppressor cells (MDSCs), M2 subtype of tumor-associated macrophages (TAMs), exclusion, and TIDE (except for LILRB5) score but a positive correlation with the T cell dysfunction score, interferon gamma (IFNG), and merck18 signatures (Figure 7C). These observations indicate that LILRBs might contribute to immune evasion through the induction of T cell dysfunction. In the TCGA AML cohort, we found that LILRBs (except for LILRB5) showed significantly higher expression in predicted responders than non-responders (Figure 7D), suggesting that AML with high LILRB expression may benefit more from ICB treatment. While no differences in PD-1 expression between low and high LILRBs expressers were observed, patients with high LILRBs showed an obviously high expression of CTLA-4 in the TCGA AML cohort (Figure S19B).

Figure 7

LILRB expression predicts responses to immunotherapy

(A) Scatterplot comparing predict performance of LILRBs to that of standardized cancer immune evasion biomarkers for OS among indicated ICB cohorts. The x axis denotes the Z score on Cox-PH regression, and the y axis indicates its significance level (two-sided Wald test). The red horizontal line in each plot indicates threshold for significance (p = 0.05). (B) Percentages of responders (complete response [CR] or partial response [PR]) and non-responders (stable disease [SD] or progressive disease [PD]) to ICB among indicated ICB cohorts between patients with high and low LILRB expression (as stratified by the median expression value of respective genes). (C) Correlogram showing the association between LILRB expression with T cell-dysfunction and T cell-exclusion signatures in TCGA AML cohort, as determined using the tumor immune dysfunction and evasion (TIDE) method. (D) Violin plots comparing the expression of LILRBs between patients who benefitted and did not benefit from immunotherapy in AML, as predicted by the TIDE algorithm. Significances were calculated by Wilcoxon rank-sum tests. See also Figures S17–S19.

The biological significance of LILRB expression in AML

We then sought to investigate the biological features associated with LILRBs in AML. Since the expressions 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 down-regulated; adjusted p < 0.05; log2 fold change [FC] ≤ -1.5 or ≥ 1.5) were differentially expressed in LILRB1high versus LILRB1low patients (Figure 8A; Table S3). Among the genes positively correlated with LILRB1 were, as expected, the other members of the LILRB family (Figure 8A). Also, genes associated with the presence of monocytes/macrophages (CD14, CD68) or M2 macrophage polarization (MSR1, MRC1, CD163) were significantly up-regulated in high LILRB1 expressers (Figure 8A), in line with our previous findings. Next, we used the STRING database to construct a protein-protein interaction (PPI) network of the differentially expressed genes (DEGs), with a confidence score >0.90. Genes interacting with LILRB1 and their subnetworks were shown through Cytoscape software (Figure 8B). We found 12 genes directly interacting with LILRB1: PILRA, TLR8, SIGLEC7, CD300C, FCGR2A, FCGR2B, FCGR3A, CD86, FGR, HCK, IL10, and ITGAX. Among them, CD300C, FCGR2A, FCGR2B, and FCGR3A also had connections with the other four LILRB members (Figure 8B). 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 the immune-system process (Figure S20A).

Figure 8

The biological significance of LILRBs expression in AML

(A) Volcano plot showing differentially expressed genes (DEGs) between high and low LILRB1 expressers (as stratified by the median expression value). (B) Cytoscape analysis of LILRB1-related network using PPI information obtained from STRING database (http://stringdb.org/). Red nodes represent up-regulated genes and blue represent down-regulated genes. (C and D) GO (C) and KEGG (D) analysis of DEGs. See also Figures S20 and S21.

We then performed Gene Ontology (GO) analysis using these DEGs, and the top 10 significant terms of biological process (BP), molecular function (MF), and cellular component (CC) enrichment analysis were shown (Figure 8C). Notably, in terms of BP, immune-response-related processes were significantly enriched, such as inflammatory response, immune-system process, and immune response. Kyoto Encyclopedia of Genes and Genomes (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 (Figures 8D and S20B).

Finally, gene set enrichment analysis (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 the LILRB1high group. Also enriched were gene sets of MLL fusion and NPM1 mutation, two distinct entities often associated with monocytic features of AML (Figure S21A). For the C7 immunologic collection, the LILRB1high group had principal enrichment in genes up-regulated in monocytes compared with other immune cells (Figure S21B), and multiple immune activities were enriched in the LILRB1high group for HALLMARK gene sets (Figure S21C).

Discussion

The LILRB family members LILRB1–LILRB5 are a group of proteins containing the immune-inhibitory ITIM motifs that negatively regulate immune cell activation. Here, using RNA-seq data of normal tissues from GTEx, FANTOM5, and HPA, we showed that LILRB members were predominantly enriched in the spleen, consistent with their immune-modulatory functions. In cancer cell lines, LILRBs showed relatively high expression in cell lines of malignant hematological origin, in line with the selective expression of LILRBs in hematopoietic lineage cells. Indeed, abnormal expression of LILRBs has been documented in various cancers, such as lung cancer, hepatocellular carcinoma (HCC), and certain types of subtypes of adenocarcinoma. In this study, based on combined datasets from TCGA and GTEx, we comprehensively analyzed LILRB expression between tumor and adjacent normal tissue across 28 cancer types (9,465 tumor and 7,831 normal samples). Our data showed that LILRBs were significantly dysregulated in the majority of tumor types. For LILRB1–LILRB4, the most striking difference was seen between AML and its normal counterparts. We also observe a strong enrichment for LILRB1–LILRB4 in the monocytic lineage; this observation was confirmed in mass spectrometry proteomic data, single-cell transcriptomics of immune cells, immune cell abundances estimated using bulk TCGA samples, and GSEA of monocyte-related gene sets, in agreement with previous reports.,,, One limitation is that many of the findings were based on correlation analyses; the results could, therefore, be biased by normalization methods and statistical analyses along the way. Future functional immunological data and prospective validation will still be required before these in silico approaches can be used in a clinical setting.

Despite being positively correlated with monocytes, LILRB1–LILRB4 were negatively correlated with the density of CD8+ T and NK cells, which are considered essential for effective anti-tumor immunity. It has been shown that activated LILRB4 on monocytic AML cells recruits SHP-2 and upregulates nuclear factor κB (NF-κB), leading to increased ARG1 and uPAR accompanied by a concomitant suppression of T cell activity., This might provide a potential mechanistic explanation to our observations. It should be noted that BM T cells in AML are often functionally impaired,11, 12, 13, possibly mediated by malignant monocyte-like cells from AML.,,, Further research aimed at unraveling the underlying molecular mechanisms is clearly warranted, as this may provide opportunities for the identification of new drug targets and therapeutics that can circumvent the T cell-suppression state in AML.

Immunosuppressive factors, such as indoleamine 2,3-dioxygenase 1 (IDO1), CD200, and TIM-3 were reported to be closely associated with a poor outcome in AML.43, 44, 45 In a preliminary analysis, Deng et al. studied the prognostic relevance of several co-stimulating and co-inhibitory receptors in the TCGA AML dataset, including LILRB1–LILRB4. Here, we independently validated the prognostic significances of LILRB members in five independent datasets. Strikingly, we showed that LILRB1–LILRB4 adversely impacted survival in almost all analyzed datasets. Of interest, we also noticed that LILRB4 was significantly associated with M2 macrophage abundances. This observation raises the possibility that LILRB4 might contribute to leukemogenesis through M2 macrophages. Our group has recently reported that M2 macrophage fractions were more selectively up-regulated in AML than the other four hematological malignancies and normal controls. Importantly, we also demonstrated superior predictive performance of the M2 marker CD206 (MRC1) than classical prognosticators in AML. Interestingly, in this study, we found that CD206 was significantly up-regulated in high LILRB1 expressers. As CD206+ and/or LILRB4+ monocytes could suppress T cell proliferation and create an immunosuppressive microenvironment in AML,, it could be hypothesized that at least part of the prognostic value of LILRBs could be attributed to the immune-suppressive TME it contributed. Acute monocytic leukemia often harbors mixed-lineage leukemia (MLL) rearrangements, an aggressive phenotype with limited treatment options and poor survival rates, which might also explain the observed result. Indeed, we demonstrated that LILRB4 was aberrantly overexpressed in MLL-rearranged AML and might be a direct target of the MLL fusion proteins.

In a recent pan-hematological-malignancies study, the authors found that LILRB2 could distinguish lymphoma and leukemia subtypes with high immune infiltration from those harboring lower cytolytic score. We consistently found multiple genes involved in immune activation (including cytokines and genes essential for microbial killing and antigen processing and presentation) were deleted in LILRB1high patients, indicating a delicate balance between immune activation and suppression in the TME.

Indeed, an integrated analysis of transcriptomic and proteomic data has uncovered and ranked LILRBs among the top potential chimeric antigen receptor (CAR) targets in AML. Gui et al. also found that blocking LILRB4 activation effectively reversed T cell suppression and inhibited AML cell infiltration. Given that LILRBs are selectively dysregulated in AML, it is tempting to speculate that AML positive for these proteins might be good candidates for immunotherapy. Indeed, we found that LILRBs showed comparable or even superior predictive power for ICB response than other biomarkers reported in the literature, as quantified both by AUC and Z score in the Cox proportional hazard (PH) regression. This agreed favorably with our observation that patients with high LILRB expression possessed primarily immune-enriched subtypes, which were characterized by highly abundant lymphocytes, increased PD-L1 expression, and amplifications of PD-L1 and PD-L2 genes, suggesting an environment with dampened anti-tumor immunity that could benefit from immunotherapy. Also, we found high LILRB expressers often exhibited increased PD-1/CTLA-4 expression in ICB cohorts. Our findings were in line with previous observations that LILRB2 was co-expressed with CTLA-4 and that patients responding to anti-PD-1 ICB showed an enriched expression of LILRB4., Importantly, Chen et al. reported that LILRB2 blockade substantially enhanced the efficacy of anti-PD-1 treatment in a human lung cancer model. We also demonstrated the predictive potential of LILRBs in AML immunotherapy. Although we observed no correlations between LILRBs and PD-1 expression in AML, LILRB agonism could reprogram TAMs into a non-suppressive and immunostimulatory phenotype, thereby enhancing the therapeutic efficacy of ICB treatment. Future cancer immunotherapy clinical trials will be critical to further validate these findings.

In this study, we provided a comprehensive analysis of the expression patterns and clinical significances of LILRBs across pan-cancers, focusing on their role in AML. We also analyzed the association of LILRB expression with genomic features and tumor immunity in AML. Our data revealed up-regulated expression of LILRBs in AML and that higher expression levels of these genes predicted worse outcomes. In addition, LILRBs were associated with an immune-suppressive TME in AML. Overall, these findings suggest important immunological and clinical implications of LILRBs in AML, which warrants further clinical investigation with immunotherapy specifically targeting AML with LILRB dysregulations.

Materials and methods

Analysis of gene-expression data

Briefly, the mRNA expression data of the LILRB family in normal tissues were obtained from the GTEx project (www.gtexportal.org/). Datasets used to assess the expression patterns of LILRBs in normal tissues and cell lines are described in detail in the supplemental methods. To determine the expression patterns of LILRBs between tumor and adjacent normal tissues across a broad range of cancer types, we systematically analyzed the gene-expression data of 9,465 tumor and 7,831 normal samples based on RNA-seq data from the TCGA and GTEx projects. All these datasets were downloaded from the UCSC Xena project and were normalized between arrays using the limma package. These studies were approved by the respective institutional review boards with written informed consent obtained from all patients.

Analysis of AML scRNA-seq data

For scRNA data analysis, previously published scRNA-seq data from 16 AML samples at diagnosis consisting of 30,712 BM cells (Van Galen AML scRNA) were downloaded from GEO (GEO: GSE116256). Another scRNA-seq data for 8 patients consisting of 30,579 AML BM cells (FIMM AML scRNA) were retrieved via the Synapse Web Portal (https://www.synapse.org and https://doi.org/10.7303/syn21991014). Data were processed and visualized using custom scripts provided by Dufva et al.

Analysis of genetic alteration data

The genetic alterations of LILRBs from TCGA PanCancer Atlas studies (10,967 patients), including somatic mutations, amplification, and deep deletion, were assessed through the cbioportal for Cancer Genomics (http://www.cbioportal.org). Procedure details are provided in the supplemental methods.

Analysis of gene-methylation data

For comparison of methylation status of LILRBs between tumor and normal samples, beta values of Illumina 450k probes at the promoter region of five genes were retrieved by the DiseaseMeth v.2.0 web portal (http://bio-bigdata.hrbmu.edu.cn/diseasemeth/analyze.html). Procedure details are provided in the supplemental methods.

Analysis of the association between LILRB4 and MLL rearrangement

We used the Hemap dataset to analyze the association between LILRB expression and common molecular subtypes. Datasets used to determine the association of LILRB4 expression with MLL rearrangement and analysis of ChIP-seq data are described in detail in the supplemental methods.

Survival analysis

We investigate the association between the expression of LILRB members and clinical outcomes across 33 cancer types. The association between transcript levels of LILRB members and OS across cancers were assessed by univariate Cox regression. To confirm the prognostic value of LILRBs in AML, we further obtained five independent GEO datasets (GEO: GSE10358, n = 304; GSE37642 [U133A], n = 422; GSE37642 [U133plus2], n = 140; GSE106291, n = 250; GSE71014, n = 104) with available survival information. Patients with AML from these datasets and the TCGA dataset were divided into those with high and low gene expression, according to the optimal cutoff determined by the X-tile method. We then performed Kaplan-Meier analysis (log rank test) to compare the survival differences of two groups regarding OS (six datasets) and EFS (only in TCGA dataset).

Immune-response analysis

The relative abundances of 22 immune cell populations in patients with AML were estimated using the CIBERSORT algorithm, as previously described.10 As CIBERSORT may not be suitable for the use of the RNA-seq data, this algorithm was exclusively applied to the TCGA LAML microarray dataset. For validation purposes, the relative fractions of immune cells were also estimated in two relatively large GEO datasets, GEO: GSE10358 and GSE6891. In addition, we used other deconvolution methods to quantify the proportions of monocytes (quanTIseq, MCP-counter, CIBERSORT abs, and xCell) and CD8 T cells (EPIC, TIMER, quanTIseq, MCP-counter, CIBERSORT abs, and xCell). These methods have been integrated as a unified interface by Sturm et al. and are freely available through the TIMER 2.0 web portal (http://timer.comp-genomics.org/). We evaluated the relationship between LILRBs and several notable immune checkpoint genes.36 Spearman correlation analysis was used to test the association between LILRB expression and these parameter estimates. Immunotherapy-associated dataset collection and analyses are provided in the supplemental methods.

Differential gene-expression analysis and functional enrichment analysis

Differential gene-expression analysis for RNA-seq data was performed using the raw read counts with the R/Bioconductor package “DESeq2,” controlled for the FDR by the Benjamini–Hochberg procedure. GO analysis and KEGG pathway analysis of LILRB1-co-expressed genes were performed using the STRING database (http://www.string-db.org/). GO and KEGG terms with FDR-corrected p values less than 0.05 were considered significantly enriched. For displaying purposes, the top 10 GO terms of each three GO categories—BP, CC, and MF—and the top 10 KEGG pathway terms were visualized as bar plots.

PPI-network analysis

We applied STRING (http://string.embl.de/) to construct a PPI network of the DEGs. We chose a confidence score >0.9 as the judgment criterion. Cytoscape visualization software (v.3.6.1) was used to present the LILRB1-related subnetwork.

GSEA

GSEA was performed on the TCGA dataset using GSEA v.4.1.0 software (http://www.broad.mit.edu/gsea). Procedure details are provided in the supplemental methods.

Statistical analysis and visualization

All statistical analyses and visualizations were performed using either indicated web servers or R v.4.0.4. For details, see the supplemental methods.

Data and code availability

The datasets analyzed in this study are available in the following open access repositories: GTEx, www.gtexportal.org/; HPA, https://www.proteinatlas.org/; CCLE, https://www.broadinstitute.org/ccle; Human Proteome Map, https://www.humanproteomemap.org/; TCGA, https://portal.gdc.cancer.gov/; UCSC Xena, https://xena.ucsc.edu; cBioPortal, http://www.cbioportal.org; GEO, https://www.ncbi.nlm.nih.gov/geo/ (GEO: GSE13159, GSE116256, GSE63409, GSE79899, GSE82116, GSE71779, GSE89336, GSE71776, GSE61785, GSE10358, GSE37642, GSE106291, and GSE71014); FIMM AML scRNA data, https://www.synapse.org (https://doi.org/10.7303/syn21991014); DiseaseMeth, http://bio-bigdata.hrbmu.edu.cn/diseasemeth/analyze.html; TIMER 2.0, http://timer.comp-genomics.org/; and TIDE, http://tide.dfci.harvard.edu/.