Integrated Pan-Cancer Analysis Reveals Distinct Clinical, Genomic and Immunological Features of the LILRB Immune Checkpoint Family in Acute Myeloid Leukemia

Background Leukocyte immunoglobulin (Ig)-like receptor Bs (LILRBs), a family of type I transmembrane glycoproteins, are known to inhibit immune activation. Methods We comprehensively evaluated the transcriptional levels and prognostic signicances 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 AML patients with altered LILRBs expression. Results Here, we show that LILRBs were signicantly dysregulated in a number of cancers, especially in acute myeloid leukemia (AML). 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 involving in immune activation were deleted in LILRB1-high patients. Epigenetically, LILRB4 was signicantly hypomethylated and marked by MLL-associated histone modications in AML. Immunologically, LILRBs were positively associated with monocytic cells including M2 macrophages, but were negatively associated with tumor-suppressive CD8 T cells. suggest important immunological and clinical implications of LILRBs in AML, which warrants further clinical investigation with immunotherapy specically targeting AML with LILRBs dysregulations.


Background
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 [4], yet chemoresistance and relapse are commonly seen in in clinical practice. Some regimens under study, such as the combination of hypomethylating agents (HMAs) and Venetoclax, have shown promising results in certain subsets of AML patients [5,6]. However, there remains an urgent need to develop novel effective therapies for various subsets of AML.
Noteworthy, 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 [7,8]; however, the transfer of immunotherapy to AML has been less successful than to other cancers [9]. Indeed, the AML microenvironment is predominantly immunosuppressive. For example, we have previously demonstrated that M2 macrophages, a classical immunosuppressive component, was preferentially enriched in AML than other hematological malignancies and normal controls [10]. Also, a recent single-cell RNA-seq study has reported proportionally fewer T cells and CTLs in AML than normal controls, and the function of these T cells are profoundly impaired, probably mediated by CD14 + monocyte-like cells [11,12]. Moreover, Noviello et al. have reported that bone marrow T cells at AML relapse showed an exhausted phenotype, which was absent in patients maintaining long-term complete response [13]. These ndings suggest encouraging therapeutic opportunities by modulating the immune environment in AML.
The 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) [14]. This group of receptors contains 5 members (LILRB1-LILRB5) that are mainly expressed in hematopoietic lineage cells and also various types of tumors [14]. As these proteins negatively regulate immune activation [15][16][17], they are often considered as immunosuppressive component in the tumor microenvironment (TME). In AML, the TMEmodulating role of LILRBs have recently came into focus, especially for LILRB4. As demonstrated by the Gui group that LILRB4 facilitates tissue in ltration of AML cells by substantially suppressing T cell activities, and blocking LILRB4 activity e ciently inhibited AML development in vitro and in vivo [18,19]. In addition, LILRB1, which showed predominant expression in monocytic AML cells as LILRB4 did [20], was found to be up-regulated in dysfunctional CD8 + T cells from AML than T cells from healthy controls [21]. Interestingly, a non-immunological AML-promoting role was reported for LILRB2, which binds Angptl2 to maintain stemness of normal stem cells and support leukemia development by inhibiting differentiation of AML cells [22]. Despite the functional importance of LILRBs in cancers, there currently lacked a systematic study to explore the expression patterns and clinical implications of all LILRB members in pan-cancers, especially in AML. Therefore, in this study, drawing on rich multi-omics data in the public domain, we comprehensively evaluated the transcriptional levels and prognostic signi cances 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 AML patients with altered LILRBs expression.

Materials And Methods
Analysis of gene expression data Brie y, the mRNA expression data of LILRB family in normal tissues were obtained from the Genotype-Tissue Expression (GTEx) project (www.gtexportal.org/) [23]. To con rm the expression patterns of LILRBs in normal tissues, we then explored the HPA (Human protein atlas) and FANTOM5 dataset from the human protein atlas database (http://www.proteinatlas.org/) [24]. Expression data of LILRBs for over 1000 cancer cell lines from various organ sites were accessed through Cancer Cell Line Encyclopedia (CCLE) (https://www.broadinstitute.org/ccle) [25]. Furthermore, RNA-seq data of 64 cell lines from The Human Protein Atlas (HPA) ((https://www.proteinatlas.org/) [24] were used to validate expression patterns of LILRBs in cancer cell lines.
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 9465 tumor and 7831 normal samples based on RNA sequencing data from the TCGA and the GTEx projects. All these datasets were downloaded from the UCSC Xena project and were normalized between arrays using the limma package [26]. To validate the differential expression of LILRBs between AML and normal controls, we further retrieved two datasets containing both healthy and AML samples from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) (accession number GSE63270 and GSE30029). We also used a gene-expression dataset of normal hematopoietic cells (GSE42519) to assess the expression patterns of LILRB family genes at various stages of normal hematopoiesis.
We used the Hemap dataset curated by Dufva et al. [27], which includes datasets of AML, pre-B-ALL, DLBCL, and MM, to analyze the association between LILRBs expression and common molecular subtypes. The molecular subtypes were de ned as previously described [27]. Three datasets-BeatAML, TCGA, and GSE13159-which have detailed cytogenetic information, were used to determine whether LILRB4 expression was associated with MLL-rearranged AML. The counts per million (CPM)

Analysis of genetic alteration data
The genetic alterations of LILRBs from TCGA PanCancer Atlas studies (10967 patients), including somatic mutations, ampli cation, and deep deletion were assessed through the cbioportal for Cancer Genomics (http://www.cbioportal.org). To determine the association between LILRBs expression and common gene mutations, patients from the TCGA LAML cohort were rst strati ed into two groups by the median expression value of respective member genes, then mutation status for 24 most frequently mutated genes were identi ed using TCGA mutational data. The relationships between mutation status and LILRBs expression were analyzed by two-sided Fisher exact tests. The mutational pro les of patients with high or low LILRB1/LILRB5 expression were displayed as co-bar plots using the "Maftools" package [29]. To detect copy number alterations (deletions and amplications) in high-and low-LILRB1 expressers, we analyzed ltered segmented copy number data in each group (Affymetrix SNP 6.0 platform) using the GISTIC 2.0 algorithm [30].

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 ve genes were retrieved by DiseaseMeth version 2.0 web portal (http://bio-bigdata.hrbmu.edu.cn/diseasemeth/analyze.html). Correlation between LILRBs methylation with expression and survival data across cancers were analyzed through the GSCALite platform [31]. For AML, the DiseaseMeth dataset contains methylation data of 271 AML (from TCGA, GSE62303, and GSE64934) and 10 normal samples (from GSE58477). Another methylation dataset-GSE63409-contains DNA methylation pro les of 20 leukemia stem cells, 24 blast cells and 30 normal hematopoietic stem and progenitor cells. Beta values of Illumina 450k probes nearest the transcription start site of the genes were selected to represent methylation level of the gene promoter area. Heatmaps displaying methylation levels of LILRBs across samples were generated using 'pheatmap' package.

Analysis of Chromatin immunoprecipitation-sequencing (ChIP-seq) data
ChIP-seq data for MLL-fusion proteins and three histone marks (H3K79me2, H3K27ac, and H3K4me3) from MV4-11 and THP-1 cells were obtained from GSE79899 [32]. For validation purpose, the three epigenetic marks from ve other ChIP-seq datasets (H3K79me2 from GSE82116 and GSE71779; H3K27ac from GSE89336 and GSE71776; H3K4me3 from GSE61785 and GSE82116) were also included in our analyses. The gene tracks were generated by uploading the wiggle les as custom tracks onto the UCSC Genome Browser, assembly hg19. Wiggle les that have failed to be loaded as custom tracks were rst converted to bigwig using the UCSC wigToBigWig tool.

Survival Analysis
We used the "Gene Outcome" module of TIMER2.0 (http://timer.cistrome.org/) [33] to 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 overall survival (OS) across cancers were assessed by univariate Cox regression. To con rm the prognostic value of LILRBs in AML, we further obtained ve independent GEO datasets (GSE10358, n = 304; GSE37642 [U133A], n = 422; GSE37642 [U133plus2], n = 140; GSE106291, n = 250; GSE71014, n = 104) with available survival information. AML patients from these datasets and the TCGA dataset were divided into those with high and low gene expression, according to the optimal cut-off determined by the X-tile method [34]. We then performed Kaplan-Meier analysis (log-rank test) to compare the survival differences of two groups regarding overall survival (six datasets) and event-free survival (EFS) (only in TCGA dataset).

Immune response analysis
The relative abundances of 22 immune cell populations in AML patients were estimated using the CIBERSORT algorithm as previously described [10]. As CIBERSORT may not be suitable for the use of the RNA-seq data [35], this algorithm was exclusively applied to the TCGA LAML microarray dataset. For validation purpose, the relative fractions of immune cells were also estimated in two relatively large GEO datasets: 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 uni ed interface by Sturm et al. [36] 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 describe by De Simone et al. [37]. Spearman correlation analysis was used to test the association between LILRBs expression and these parameter estimates.

Differential gene expression analysis and functional enrichment analysis
Differential gene expression analysis for RNA sequencing data was performed using the raw read counts with the R/Bioconductor package "DESeq2", controlled for the false discovery rate (FDR) by the Benjamini-Hochberg procedure. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of LILRB1-coexpressed genes were performed using the STRING database (http://www.string-db.org/). GO and KEGG terms with false discovery rate (FDR)-corrected p values less than 0.05 were considered as signi cantly enriched. For displaying purposes, the top 10 GO terms of each three GO categories-biological process (BP), cellular component (CC), and molecular function (MF), and the top 10 KEGG pathway terms were visualized as bar plots.

Protein-protein interaction (PPI) network analysis
We applied STRING (http://string.embl.de/) to construct a protein-protein interaction (PPI) network of the differentially expressed genes (DEGs). We chose a con dence score > 0.9 as the judgment criterion.
Cytoscape visualization software (version 3.6.1) was used to present the LILRB1-related sub-network.

Gene Set Enrichment Analysis (GSEA)
Gene set enrichment analysis (GSEA) was performed on the TCGA dataset using GSEA v4.1.0 software (http://www.broad.mit.edu/gsea). Statistical signi cance of GSEA results was determined by 1,000 gene set permutations, with signal-to-noise gene ranking. All the gene sets used in this study were obtained from GSEA MSigDB website (http://software.broadinstitute.org/gsea/msigdb/index.jsp), and the gene sets were considered to be signi cantly enriched at a false discovery rate < 0.25 and normalized P-value < 0.05. Three categories of gene sets were used in this study: C2, curated gene sets containing genes coregulated in response to speci c perturbations; C7, immunologic signature gene sets that represent cell states and perturbations within the immune system; and H, hallmark gene sets which represent wellde ned biological states or processes.

Statistical analysis and visualization
Wilcoxon rank sum tests were used to compare differences between two groups. Speci cally, for Fig. 6a, we determined differential expression of LILRBs between each molecular subtype and the remaining samples for each disease subtype. The average fold changes (FCs) and Bonferroni-adjusted p-values (false discovery rate [FDR]) were computed using Wilcoxon rank sum tests. The FDR values were then combined for each LILRB members using Stouffer's method and categorized into ve groups based on signi cance cutoffs for visualization (0.05, 0.01, 0.001, 1e-5, 1e-16). All statistical analyses and visualizations were performed using either indicated web servers or R version 4.0.4. Speci cally, the expression map of LILRBs was generated using the "gganatogram" package [38], the box, bar, scatter and bubble plots were produced with the R package "ggplot2", "ggpubr" and "ggsci", the volcano plots were generated using the "EnhancedVolcano" package, and survival curves were made using the "survival" package. The correlation matrix was calculated and visualized using the "corrplot" or "ggcor" library.
Finally, "circlize" was used to create the chord diagrams. All statistical tests were two-sided with p-values less than 0.05 considered signi cant.

Expression patterns of LILRBs in normal tissues and cancer cell lines
In our rst 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 pro les 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 con rmed 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 ndings indicated a cellular-, tissue-, and disease-speci city 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 signi cantly 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 signi cantly 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 signi cantly up-regulated in AML and dynamically expressed during normal hematopoiesis Our observations, together with previous ndings [18], re ect an AML-speci c 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 upregulated 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, ampli cations, and deletions) frequencies of LILRBs across pan-cancers. The average alteration frequencies of ve 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 ve 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).
Ampli cations 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 rst 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 signi cantly 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, signi cantly 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 signi cant 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 signi cance and direction of the prognostic signi cances 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 ve 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 ve 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 speci c 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 ve 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 con rm 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 con rm the relevance of LILRB4 expression in MLL-rearranged AML, we collected four MLLrearrangement-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 LILRB4expressing samples with those in low LILRB4-expressing samples. We found gene-sets down-regulated in MLL-rearranged AML (MULLIGHAN_MLL_SIGNATURE_1_DN) showed signi cantly lower ssGSEA scores in LILRB4-high patients than in LILRB4-low patients; whereas for gene-sets up-regulated in MLLrearranged 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 signi cantly 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 signi cant 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 ve 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) classi cation 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 pro les characterized for AML, we identi ed signi cantly mutated genes occurred in patients with high and low LILRB1-LILRB4 expression (as strati ed 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 signi cantly deleted regions and four signi cantly ampli ed region (FDR = 0.25) (Fig. 7d). Interestingly, the majority of genes deleted in LILRB1-high patients were involved in in ammatory 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 signi cantly deleted. In contrast, LILRB-high AML patients had recurrent ampli cation at loci essential in AML pathogenesis, including KMT2A and ERG [47,48] (Fig. 7d).

Correlations Between LILRBs and Tumor Immune In ltrating 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 in ltration 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 nding that LILRBs were preferentially expressed in monocytic AML [19,49]. This monocytic preference was also con rmed 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 in ltrating 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 ndings further con rmed 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 signi cant 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 signi cance of LILRBs expression in AML
We then sought to investigate the biological features associated with LILRBs in AML. Since the expression of ve LILRB members were highly correlated, a comparison of gene expression pro les 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 LILRB1 high versus LILRB1 low 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 signi cantly up-regulated in high LILRB1 expressers (Fig. 9a), in line with our previous ndings. Next, we used STRING database to construct a protein-protein interaction (PPI) network of the differentially expressed genes (DEGs), with a con dence score > 0.90. Genes interacted with LILRB1 and their subnetworks 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 signi cant terms of BP, MF and CC enrichment analysis were shown (Fig. 9c). Notably, in terms of BP, immune response-related processes were signi cantly enriched, such as in ammatory 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 LILRB1 high and LILRB1 low 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 LILRB1 high 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 LILRB1 high 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 LILRB1 high group for HALLMARK gene sets (Fig. 10c).

Discussion
The LILRB family members-LILRB1-5-are a group of proteins containing the immune-inhibitory ITIM motifs which negatively regulate immune cell activation [14]. 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 [50], hepatocellular carcinoma (HCC) [51], and certain types of subtypes of adenocarcinoma [41]. In this study, based on combined datasets from TCGA and GTEx, we comprehensively analyzed LILRBs expression between tumor and adjacent normal tissue across 28 cancer types (9465 tumor and 7831 normal samples). Our data showed that LILRBs were signi cantly dysregulated in the majority of tumor types. For LILRB1-4, the most striking difference was seen between AML and its normal counterparts. Although previous study has reported increased expression of LILRB4 in AML [19,49], this was not seen in two independent validation sets of AML patients, probably due to technical/biological heterogeneity across studies. However, we do note a strong enrichment for LILRBs in the monocytic lineage; this observation was con rmed in HPA pan-cancer cell lines dataset, single-cell transcriptomics of immune cells, immune cell abundances estimated using bulk TCGA samples, and GSEA analysis of monocyte-related gene sets, in agreement with previous reports [19,20,22,49].
Despite being positively correlated with monocytes, LILRB1-4 were negatively correlated with the density of CD8 + T and NK cells, which are considered essential for effective anti-tumor immunity [27]. It has been shown that activated LILRB4 on monocytic AML cells recruits SHP-2 and upregulates NFκB, leading to increased ARG1 and uPAR accompanied by a concomitant suppression of T cell activity [18,19]. This might provide a potential mechanistic explanation to our observations. It should be noted that bone marrow (BM)-T cells in AML are often functionally impaired [11][12][13]52], possibly mediated by malignant monocyte-like cells from AML [11,19,21,53]. Further research aimed to unravel the underlying molecular mechanisms is clearly warranted, as this may provide opportunities for identi cation of new drugs 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 [54][55][56]. In a preliminary analysis, Deng et al. studied the prognostic relevance of several co-stimulating and co-inhibitory receptors in TCGA AML dataset, including LILRB1-LILRB4 [19], but there is limited data and require validation in large datasets. Here, we independently validated the prognostic signi cances of LILRB members in ve independent datasets. Strikingly, we showed that LILRB1-4 adversely impacted survival in almost all analyzed datasets. Of interest, we also noticed that LILRB4 was signi cantly 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 macrophages fractions were selectively upregulated in AML than other four hematological malignancies and normal controls [10]. Importantly, we also demonstrated a superior predictive performance of the M2 marker CD206 (MRC1) than classical prognosticators in AML. Interestingly, in this study we found that CD206 was signi cantly up-regulated in high LILRB1 expressers. As CD206 + and/or LILRB4 + monocytes could suppress T-cell proliferation and create an immunosuppressive microenvironment in AML [19,53], it could be hypothesized that at least part of the prognostic value of LILRBs could be attributed to the immunesuppressive 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 have found that LILRB2 could distinguish lymphoma and leukemia subtypes with high immune in ltration from those harboring lower cytolytic score [27]. We consistently found multiple genes involving 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 [57]. Preliminary evidence from the Gui group also suggests that blocking LILRB4 activation effectively reversed T-cell suppression and inhibited AML cell in ltration [18]. 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. Future cancer immunotherapy clinical trials will be critical to further validate these ndings.
In this study, we provided a comprehensive analysis of the expression patterns and clinical signi cances of LILRBs across pan-cancers, focusing its role in AML. We also analyzed the association of LILRBs 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 outcome. In addition, LILRBs were associated with an immune-suppressive TME in AML. Overall, these ndings suggest important immunological and clinical implications of LILRBs in AML, which warrants further clinical investigation with immunotherapy speci cally targeting AML with LILRBs dysregulations.  The color depicts the log2-transformed fold change (Log2FC) between tumor and normal tissues. *P < 0.05; **P < 0.01; ***P < 0.001. (e) Bar plot showing genes signi cantly upregulated and downregulated (P < 0.05) across different cancer types. Red, up-regulated expression; blue, down-regulated expression. (f and g) Box plots showing expression levels of LILRB1-4 in normal controls and AML in the GSE63270 (f) and GSE30029 (g) datasets.