What this is
- Esophageal squamous cell carcinoma (ESCC) is characterized by a complex immune microenvironment.
- This research provides a detailed immune cell atlas of ESCC using single-cell RNA sequencing.
- Key findings include the presence of exhausted T and NK cells, , and macrophages contributing to immune suppression.
- The study identifies potential immunosuppressive mechanisms that could be targeted for therapeutic interventions.
Essence
- ESCC features a suppressive immune microenvironment dominated by exhausted T and NK cells, , and macrophages, indicating challenges for effective immunotherapy.
Key takeaways
- Exhausted CD8 T cells progress from a pre-exhausted state to an exhausted state, indicating a dynamic immune response within the .
- The study reveals significant inter-tumoral heterogeneity in immune cell infiltration, with some tumors showing increased T cell and macrophage presence compared to adjacent tissues.
- Crosstalk between macrophages and may enhance immune suppression, suggesting that targeting these interactions could improve anti-tumor immunity.
Caveats
- The study's findings are based on a limited sample size of seven patients, which may affect the generalizability of the results.
- The analysis primarily focuses on immune cell populations and does not address the broader tumor biology or microenvironmental factors.
Definitions
- exhausted T cells: T cells that have lost their ability to effectively respond to antigens due to persistent stimulation in the tumor microenvironment.
- regulatory T cells (Tregs): A subset of T cells that modulate the immune response, often promoting tolerance and preventing autoimmunity.
- tumor microenvironment (TME): The environment surrounding a tumor, including immune cells, stromal cells, and signaling molecules that influence tumor behavior.
AI simplified
Introduction
Esophageal cancer is one of the most commonly diagnosed and deadly cancer types, especially in East Asia1. However, esophageal cancer is significantly understudied compared with other common tumor types, and in the recent decades there has been limited progress in therapeutics. Histologically, esophageal cancer can be classified into two subtypes: adenocarcinoma (EAC) and squamous cell carcinoma (ESCC)2. ESCC is the dominant subtype and accounts for ~90% of esophageal cancer cases worldwide3. Esophageal cancer is also among the tumor types with the highest median mutation burden, and is ranked higher than kidney, head and neck, and colorectal cancer4. Recently, PD-1 antibodies, pembrolizumab and nivolumab have been used in clinical trials for subsets of patients with advanced ESCC for whom first-line chemotherapy failed. However, they show only moderate improvement in the overall survival compared with chemotherapy5,6. A systematic interrogation of infiltrating immune cells in ESCC will help to profile the immune status of ESCC, evaluate the application of current checkpoint blockades, and, most importantly, lead to innovative immunotherapies.
Single-cell transcriptome analysis of immune cells in tumors provides a way to comprehensively study these cells in a highly complex tumor microenvironment (TME). Recently, single-cell RNA sequencing (scRNA-seq) has been applied to tumor-infiltrating immune cells isolated from limited types of cancers, including cutaneous melanoma7,8, non-small cell lung cancer9,10, hepatocellular carcinoma11,12, basal cell carcinoma13, colorectal cancer14,15, and breast cancer16. These studies uncover significant inter-tumoral and intra-tumoral heterogeneity in tumor immune profiles, diverse immunosuppressive populations, less-defined immune cell subsets, and signal transduction networks in related cancer types. However, such analyses have not been applied to ESCC.
To describe the immune landscape of ESCC, we used high-dimensional scRNA-seq to total immune cells isolated from seven surgically removed ESCC tumors and their matched adjacent tissues. T cell receptor (TCR) sequencing was also conducted to retrieve information on T cell clonality. Our analyses revealed inter-tumoral heterogeneity among individual ESCC patients. A subgroup of ESCC tumors presented significantly increased infiltration and clonal expansion of T cells, compared with their matched adjacent tissues. However, we identified exhausted T cells, exhausted NK cells, regulatory T (Treg) cells, alternatively activated macrophages (M2), and tolerogenic dendritic cells (tDCs) in these tumors, indicating an inflamed but immune-suppressed TME in ESCC. Transcriptional profile coupled TCR-sequencing revealed lineage connections in CD4 and CD8 T cells populations. We discovered exhausted CD8 T cells showing continuous progression from a pre-exhausted state to an exhausted state. Exhausted CD4, CD8 T, and NK cells were major proliferative cell components in the TME. Additionally, we identified crosstalk among macrophages and Tregs through ligand–receptor interactions that may contribute to the immune suppressive state and disease progression. Furthermore, we identified a gene signature that was significantly associated with the survival of patients with ESCC. Our results comprehensively characterized tumor-infiltrating immune cells, revealed the landscape of the suppressive immune state, and set the baseline for applying and developing immunotherapies for ESCC.
Results
scRNA-seq of immune cells isolated from ESCC
To enable a systematic analysis of immune cell populations, we normalized and pooled single-cell data from all samples and conducted unsupervised clustering to identify distinguishable populations. The whole procedure was performed using Seurat v3.0 with default parameters17. We annotated these populations using their canonical markers and successfully identified the major types of tumor-infiltrating immune cells as shown in other cancers, including T cells, NK cells, monocytes/macrophages, dendritic cells (DCs), B cells, plasma cells, and mast cells, as well as a very small fraction (1.31%) of other non-immune cells that were mixed in with the sorted cells (Fig. 1b). The expression of classic markers of these cell types was consistent with the annotation (Fig. 1c, d). We then analyzed “other” cluster form tumors, and found that most cells had copy number variations (CNVs), including both amplifications and deletions, suggesting that this cluster included tumor cells (Supplementary Fig. 2)
By comparing the percentages of each cell type in CD45+ cells between tumor and adjacent tissues, we found an increase of T cells and monocytes/macrophages in tumors. In contrast, the percentages of B and NK cells were decreased (Fig. 1e and Supplementary Fig. 3a). In agreement with recent studies18, we found a large degree of variation in the immune composition among tumors (Fig. 1f, g, and Supplementary Fig. 3b). T lineage cells were the most abundant immune cell type in most tumors, making up 30–71% of the total CD45+ cells (Fig. 1g). However, considering the ratios of each immune cell type to all cells analyzed by flow cytometry during CD45+ cell isolation, there was high variation between matched tumor and adjacent tissues, as well as among individuals (Supplementary Data 1). Seven pairs of samples were roughly divided into two groups. There were only minor differences between the matched adjacent and tumor tissues in three tumor-adjacent tissue pairs (S133, S134, and S150). T cells made up to fewer than the 2% of total cells in these tumors. In contrast, the immune profiles of four other tumor-adjacent pairs (S135, S149, S158, S159) presented a significant shift in a PCA, in which 6–12% of total cells were T cells in tumors (Fig. 1h, i). These tumors also showed increased numbers of monocytes/macrophages, compared with other tumors and adjacent tissues (Supplementary Fig. 3c). In addition, we found inter-patient variation in biologic signatures, including hypoxia, inflammation response, and TNFA-via NFKB pathways in lymphocytes. Interestingly, S135 and S158 showed similar gene signatures enrichment, and S133 and S134 showed similar gene signatures enrichment in these pathways (Supplementary Fig. 3d–f).
Next, we further validated our results for the major immune cell types with additional samples by flow cytometry and immunohistochemistry (IHC). We found an increase in T cells and macrophages and a decrease in NK and B cells in tumors, compared to adjacent tissues, which is consistent with the scRNA-seq data (Supplementary Figs. 4 and 5). Notably, neutrophils were not identified in scRNA-seq as a population like others reported12,18–20, but they were detected in low abundance by flow cytometry and IHC. The failure to detect neutrophils in scRNA-seq may be caused by the combination of the low abundance of neutrophils in ESCC and the limitation of the current 10× scRNA-seq technique. Neutrophils’ low RNA content and abundance of RNases may lead to increased sensitivity to prolonged processes of scRNA-seq, which could potentially result in fewer transcripts being detected, resulting in these cells not passing quality control.
Next, we compared the major compartments of infiltrating immune cells in ESCC to other cancer types with available data7,12,18–22. We found that ESCC was among the tumor types with a higher number of infiltrating T cells and monocytes/macrophages and a lower number of infiltrated B cells (Fig. 1j). This is consistent with our observation that ESCC had increased T cells and monocytes/macrophages and decreased B cell ratios, compared to their adjacent tissues (Fig. 1e). Notably, recent studies suggested positive effects of tumor-infiltrating B cells, especially those in tertiary lymphoid structures, on increasing response to immunotherapies23. Whether it is responsible for the response of checkpoint blockade in ESCC needs an additional investigation.
Profile of immune infiltrates in human ESCC with scRNA-seq and scTCR-seq. Schematic diagram of the experimental design and analysis.UMAP plot of 80,787 high-quality immune cells to visualize cell-type clusters based on the expression of known marker genes.Heatmap of the relative expression level of genes across cells, sorted by cell type. The expression was measured as the-score normalized log2 (count+1).Expression levels of relative marker genes across 80,787 cells illustrated as UMAP plots. The expression was measured as the log2 (count+1).Pie charts of cell-type fractions for ESCC tumor and adjacent tissues’ infiltrating immune cells, colored by cell type.UMAP plot of complete immune systems from seven ESCC tumor and adjacent tissues, representative of S133 and S159; other samples were shown in Supplementary Fig.. Cells were colored by clusters and labeled with the inferred cell types.The proportion of cells that contributed to each cluster by each sample, colored by cell type. A, Adjacent; T, Tumor.PCA was performed on the absolute ratio of cell types of 14 samples from seven patients to see the heterogeneity between samples, colored by sample.The absolute T cell ratio in total cells for the 14 samples.Comparison of the fractions of T cells, B cells, NK, monocytes/macrophages in tumors from patients with ESCC (= 7), BRAC (= 8), LC1 (= 7), LC2 (= 11), SKCM (= 16), HCC (= 3), COAD (= 10),n and KIRC (= 3). Skin cutaneous melanoma (SKCM); breast invasive carcinoma (BRCA); lung cancer (LC); hepatocellular carcinoma (HCC); colon cancer. (COAD); kidney cancer (KIRC). Each box represented the interquartile range (IQR, the range between the 25th and 75th percentiles) with the mid-point of the data, and whiskers indicate the upper and lower values within 1.5 times the IQR.value was calculated by two-tailed Wilcoxon sum rank test. a b c d e f g h i j z n n n n n n n n P 3b
Clustering and subtype analyses of T and NK cells
To investigate gene networks in cytotoxic and exhausted CD8 T cells and Treg cells, we used the public naïve, Treg, exhuastion, and cytotoxic signatures10 (Supplementary Data 3) and applied these signatures to CD8 and CD4 clusters and computed a transcriptional score. In CD4 T cells, CD4-C6-FOXP3 had the strongest Treg signature, while CD4-C5-STMN1 had the strongest exhaustion signature; CD4-C4-IFIT3 was enriched in cytotoxic signature (Fig. 2c). Consistent with previous analyses, CD8-C1-NKG7 was the most active cytotoxic CD8 T cells, whereas CD8-C5-CCL5 and CD8-C6-STMN1 had lower exhaustion scores than CD8-C7-TIGIT (Fig. 2d). The naïve score was very low in CD8 T cell clusters, which is consistent with Fig. 2b; no naïve clusters of CD8 T cells were identified. This result suggested that most of the tumor-infiltrating CD8 T cells were in the active, memory, or exhausted states in ESCC (Fig. 2d). We further analyzed the genes whose expression was highly correlated with the expression of FGFBP2, LAG3, and FOXP37. The top 50 genes were then used as signatures for cytotoxicity, exhaustion, and Treg with the top 30 genes shown in Supplementary Fig. 6c. We then use these signatures to analyze T cells clusters and found that the enrichment scores were consistent with the published signatures (Supplementary Fig. 6d, e).
There were lineage connections within CD4 populations. Some of the genes activated in Tregs (CD4-C6-FOXP3) overlapped with genes characteristic of the exhaustion program in CD4 T cells (CD4-C5-STMN1). We compared the gene expression of both clusters with naïve-like cell population (CD4-C1-CCR7). Genes enriched in both exhausted and Treg cells included regulatory molecules and many co-inhibitory and co-stimulatory receptors, such as TNFRSF9, CSF1, and TIGIT. In contrast, CD4-C6-FOXP3 expressed much higher levels of FOXP3, IL2RA, and CTLA4 than CD4-C5-STMN1, while CD4-C5-STMN1 expressed higher levels of CCL5, CCL4, IFI6, TOX, PDCD1, CXCL13, IFNG, and ID2 than CD4-C6-FOXP3 (Fig. 2e and Supplementary Fig. 6f). Visualization of the exhaustion and Treg scores confirmed the overlap between these two clusters (Fig. 2c). Cytotoxic and exhausted CD8 T cells both expressed many effector molecules such as GNLY and GZMH, while exhausted CD8 T cells expressed a higher level of IFNG than cytotoxic cells (Fig. 2f), which suggested that exhausted T cells still expressed high levels of some effector molecules and tried to respond to tumor cells. Different stages of CD8 T cell exhaustion have been reported in other cancer types26; pre-exhausted T cells expressed intermediate levels of PD1, EOMES, TOX, and NFAT2, while exhausted T cells expressed high levels of these molecules. We found that CD8-C5-CCL5 expressed lower levels of PDCD1 and EOMES than CD8-C7-TIGIT, while CD8-C6-STMN1 was in the middle (Fig. 2g). We also checked the expression levels of TOX and NFATC2; both genes participate in establishing epigenetic programs to install permanent exhaustion status in CD8 T cells27,28. The data showed that CD8-C7-TIGIT expressed high levels of TOX and NFATC2 (Fig. 2g). This suggests that CD8-C5-CCL5 were at an early stage of exhaustion, while CD8-C7-TIGIT was in the exhaustion stage and CD8-C6-STMN1 was likely a transition stage between CD8-C5-CCL5 and CD8-C7-TIGIT.
Dr. Simoni and colleagues reported that large fractions of tumor-infiltrating CD8 T cells are bystanders that recognize cancer unrelated epitopes. These cells lack CD39 and are phenotypically distinct from tumor antigen-specific CD8 T cells29. We then analyzed CD39 expression in CD8 T cells in our samples. The results showed that CD39 expression was significantly higher in pre-exhausted and exhausted CD8 T cells (C5, C6, C7), compared to other clusters. Additionally, most CD8 T cells expressed higher CD39 in tumors than adjacent tissues except S158 (Fig. 2h). CD39-CD8 T cells may be bystander cells and recognize non-tumoral antigens, such as Epstein-Barr virus, human cytomegalovirus, or influenza virus, which are commonly found in the esophagus.
Detailed characterization of T cells in ESCC. UMAP plot of 44,634 single T and NK cells from 14 samples, showing the formation of 17 main clusters, including six for CD4 T cells, seven for CD8 T cells, one for CD4 and CD8 negative T cells cluster, and three for NK cells. Each dot corresponds to a single cell, colored according to the cell cluster.Heatmap of-score normalized log2 (count+1) expression of selected T cell function-associated genes in each cell cluster.Dot plot of representative cytotoxic, exhaustive, naïve, and Treg signatures in CD4 T cell clusters,-score normalized log2 (count+1).Dot plot of representative cytotoxic, exhaustive, and naïve signatures in CD8 T cell clusters,-score normalized log2 (count+1).,Violin plot showing the,,,,,,,, andin CD4-C6-FOXP3 and CD4-C5-STMN1 cells,, andin CD8-C7-TIGIT and CD8-C1-NKG7 cells. The expression was measured as the log2 (count+1).Violin plot showing,,, andin the CD8-C5-CCL5, and CD8-C6-STMN1, CD8-C7-TIGIT clusters. The expression was measured as the log2 (count+1).Violin plot showingexpression in CD8+ T cell clusters from 7 paired ESCC and adjacent samples. The expression was measured as the log2 (count+1). Data are presented as mean ± SD;value was calculated by two-tailed Student’s-test. a b c d e f e f g h Z Z Z CXCL13 CCL4 CCL5 IFI6 ID2 IFNG TNFRSF9 FOXP3 TIGIT IFNG GNLY GZMH PDCD1 NFATC2 EOMES TOX ENTPD1 P t
Altered status of T and NK cells in tumors
We also observed a substantial decrease in NK cells in tumor tissues, compared with matched adjacent tissues (Fig. 1e). Additionally, the major cluster of NK cells switched from NK-C1-NCR3 in adjacent tissues to NK-C3-KLRC1 in tumors, and NK-C2-STMN1 also increased dramatically in tumors (Fig. 3g). NK-C1-NCR3 expressed high levels of NCR3, CD266, NKG7, and LAMP1 (Fig. 3h). In contrast, NK-C3-KLRC1 and NK-C2-STMN1 clusters expressed KLRC1 and ITGA1 inhibitory receptors at high levels (Fig. 3h). Flow cytometry assays verified increased NKG2A (KLRC1) expression in NK cells in ESCC compared with adjacent tissues (Fig. 3i). Indeed, NK-C3-KLRC1 and NK-C2-STMN1 had extremely low-cytotoxic scores; in contrast, the exhaustion scores were elevated (Fig. 3j), which indicated that NK cells were insufficient and function impaired in ESCC.
We further analyzed the cell cycle in T and NK cell clusters to determine the proliferating ability of cells. We generated a proliferation score of cell cycle genes that were previously shown to denote G1/S or G2/M phases30,31 and used it to infer the proliferation status of T and NK cell clusters. Interestingly, we found that CD4-C5-STMN1, CD8-C6-STMN1, and NK-C2-STMN1, which were enriched in exhaustion genes, were all highly proliferative (Fig. 3k–m), which is consistent with recent research suggesting that exhausted T cells are the major intra-tumoral proliferating immune cell compartment7.
Altered status of T and NK cells in tumors. Pie charts of CD4 T cell cluster fractions for ESCC tumors and adjacent tissues, colored by cluster type.Proportions of clusters of CD4 T cells in the tumors and adjacent tissues (= 7). Each box represents the interquartile range (IQR, the range between the 25th and 75th percentile) with the mid-point of the data, whiskers denote 1.5 times the IQR. Two-tailed paired Student’s-test.Proportions of CD4+CD25hiFOXP3+ Tregs in ESCC and adjacent tissues (= 9), Data are presented as mean ± SEM;value was calculated by two-tailed paired Student’s-test.Pie charts of CD8 T cell cluster fractions for ESCC tumors and adjacent tissues, colored by cluster type.Proportions of clusters of CD8 T cells in tumor and adjacent tissues (= 7). Each box represents the interquartile range (IQR, the range between the 25th and 75th percentiles) with the mid-point of the data; whiskers indicate the upper and lower values within 1.5 times the IQR. Two- tailed paired Student’s-test.Flow cytometry measured CD8 T cell PD1 expression in ESCC and adjacent tissues (= 7). Data are presented as mean ± SEM;value was calculated by two-tailed paired Student’s-test.Pie charts of NK cell cluster fractions for ESCC tumor and adjacent tissues, colored by cluster type.Violin plots comparing the indicated gene expression in NK cell clusters. The expression was measured as the log2 (count + 1).Flow cytometry measured NK cells NKG2A expression in ESCC and adjacent tissues (= 7). Data are presented as mean ± SEM;value was calculated by two-tailed paired Student’s-test.Dot plot of cytotoxic and exhaustion signatures in NK clusters.-score normalized log2 (count + 1).Dot plot of proliferation signatures in CD4, CD8, and NKclusters.-score normalized log2 (count + 1). a b c d e f g h i j k–m k l m n t n P t n t n P t n P t Z Z
Clonality of CD4 and CD8 T cells
Clonality of CD4 and CD8 T cells. The association between the number of T cell clonotypes and the number of cells per clonotype. The dashed line separates non-clonal and clonal cells; with the latter being identified by repeated usage of TCRs. Solid line LOESS fitting shows the correlation between the two axes.The TCR distribution of T cells across different samples. Unique (= 1), and clonal (= 2,≥ 3) TCRs are labeled with different colors. A, Adjacent; T, Tumor.Representative examples of dominant clonotypes (top 10 in color) from each tumor (gray) identified by TCR sequencing; other samples are shown in Supplementary Fig..Bar plot showing the average fractions of clone T cells among all 14 samples in each cluster, data is presented as mean ± SEM. Dot colors represent different samples. A, Adjacent; T, Tumor.,Percentage of clonal cells in CD8-C1-NKG7and CD4-C6-FOXP3from tumor and adjacent tissues (= 7). Data are presented as mean ± SEM,value was calculated by two-tailed paired Student’s-test.,Share TCR clone types between different clusters in CD4 T cellsand CD8 T cells. The lines connecting the dots indicate the clusters sharing TCR clone types, the bar plot shows the shared TCR clone type number, and dot only suggests the unique clone types.The number (left panel) and percentage (right panel) of shared TCRs between CD8-C7-TIGIT and other CD8 clusters in tumors.The number (left panel) and percentage (right panel) of shared TCRs between CD4-C6-FOXP3 and other CD4 clusters in tumors. a b c d e f e f g h g h i j n n n n P t 6a
Distinct functional composition of myeloid cells in ESCC
We further applied single-cell regulatory network inference and clustering (SCENIC) method38 to explore the transcription factors that may regulate monocyte, M1, and M2, development. MITF, BHLHE40, ATF3, and USF2 were upregulated in M2, whereas IRF transcription factors, including IRF1, IRF7, IRF2, IRF5, and PRDM1 were upregulated in M1. RARA, FOSB, and NFKB2 were greatly increased in the monocyte clusters (Fig. 5h and Supplementary Fig. 9h). Strikingly, SCENIC also revealed a dichotomy between most tumor and adjacent tissue pairs, except S158 (Fig. 5i). The transcriptional factor BHLHE40 was specifically expressed in M2. We conducted the network analyses to identify the BHLHE40 downstream genes and analyzed the functions through Metascape (Fig. 5j). Data showed that BHLHE40 downstream genes were associated with myeloid cell differentiation, and negatively regulated the cellular response (Supplementary Fig. 9i). BHLHE40 has recently been reported to mediate tissue-specific control of macrophage self-renewal, and proliferation39 and we found that BHELHE40 was associated with poor prognosis of ESCC (Supplementary Fig. 9j). BHLHE40 may play a critical role in inducing TAMs toward the M2 phenotype, and further studies are needed to explore the detailed mechanisms. Our data suggested that an enrichment of suppressive TAMs in the ESCC microenvironment may contributed to the progression of disease, and they elucidated some compelling TF candidates associated with prognosis.
Detailed characterization of monocyte/macrophage cells. UMAP plot of 17,256 myeloid cells from 14 samples, showing the formation of the 14 main clusters, including nine for monocytes/macrophages (mono/Macro) and five for dendritic cells (DCs). Each dot corresponds to one single-cell, colored according to the cell cluster.Dot plot of representative M1, M2, monocyte, and MDSC signatures in the monocytes/macrophage clusters.-score normalized log2 (count + 1).The trajectory of Macro-C3-CSF1, Mono-C1-VCAN, Macro-C2-IL1RN, Macro-C1-IL6 state transition in a two-dimensional state-space inferred by Monocle. Each dot corresponds to one single-cell, colored according to its cluster label. Arrows show the increasing directions of certain cell properties.A gene co-expression network of single cell from monocytes/macrophages was constructed by weighted correlation network analysis. The heatmap showed the topological overlap matrix among all genes used in the analysis. The darker color represented a higher overlap. The hierarchical clustering and module assignment of genes were shown along the left side and top.Module–trait relationships in all monocytes/macrophages sub-cluster. Within each heatmap, red indicated a positive correlation and green indicated a negative correlation. The numbers in brackets were correlationvalues. The red box indicated the Turquoise module was highlighted.The Turquoise module membership and gene significance for Mono-C1-VCAN. The most correlated genes were labeled.Kaplan–Meier overall survival curves of TCGA ESCC patients with the top 50 most correlated genes generated in.-values were calculated using the two-sided log-rank test.Heatmap of thevalues of AUC scores of expression regulation by transcription factors of the indicated clusters, as estimated using SCENIC.Heatmap of thevalues of AUC scores of expression regulation by transcription factors of the different samples, as estimated using SCENIC. A, Adjacent; T, Tumor.Network analysis of the transcription factor BHLHE40 and its target genes through iGraph. The size of the circle represents the co-expression relationship score. a b c d e f g f h i j Z P P t t
Detailed characterization of DC. Heatmap of five DC cluster marker genes. The top bar indicated the origins, and the middle bar indicated the clusters, and the lower bar indicated the different patients, the color intensity indicated the expression level. The expression was measured as-score normalized log2 (count + 1).Box plots of each dendritic cluster between the tumor and the adjacent tissues (= 7). Each box represented the interquartile range (IQR, the range between the 25th and 75th percentiles) with the mid-point of the data; whiskers indicated the upper and lower values within 1.5 times the IQR.value was calculated by two-tailed paired Student’s-test.Dot plot representative of the activation, migration, and tolerogenic signatures in dendritic clusters. The expression was measured as-score normalized log2 (count + 1).Volcano plot showing the differentially expressed genes (LAMP3DC versus other DC clusters), each red dot denoted an individual gene with an adjustedvalue < 0.01 and fold change > =2, two-sided Wilcoxon test.FACS analysis of LAMP3 and PDL1 expression in DC cells from adjacent and ESCC tumor tissues (= 5). Data were presented as mean ± SEM;value was calculated by paired two-tailed-test.CCR7 and CD83 expressions in LAMP3+ and LAMP3− DC cells from adjacent and ESCC tumor tissues were measured by FACS (= 4). Data were presented as the mean ± SEM.value was calculated by paired two-tailed-test.Detection of CD11CLAMP3PDL1IDODCs in ESCC tissue by multi-color IHC staining. Representative data from three patients were shown.DC cells were stimulated with or without LPS and IFNγ for 24 h, and the expression of PDL1, LAMP3, and IDO was analyzed. One of the three similar experiments was represented; data were presented as the mean ± SEM,value was calculated by two-tailed Student’s-test.DC cells were stimulated with or without LPS and IFNγ or not for 24 h and cultured with CD4CD45RAnaïve T cells for 4 days, and FOXP3 expression were measured. One of the six similar experiments was represented, data were presented as the mean ± SEM, andvalue was calculated by two-tailed Student’s-test.Heat-map of thevalues of AUC scores of expression regulation by transcription factors of the indicated clusters, as estimated using SCENIC. a b c d e f g h i j Z n P t Z P n P t n P t P t P t t + + + + + + +
Cell–cell interaction between immune cells in ESCC
Cell–cell communications, mainly through ligand–receptor interactions, play key roles in determining the TME and responding to therapeutics. Thus, we performed systematic analyses on the potential cell–cell interactions based on co-expression of known ligand–receptor pairs in any two types of tumor-infiltrating immune cells43 and compared them between tumor and adjacent tissues. Because of the high infiltration levels and the important immune regulatory roles of macrophage and Tregs in tumors, we applied scTHI (https://github.com/miccec/scTHI↗), another widely used method to analyze macrophages and Tregs interactions and found that both of them had comparable results (Supplementary Data 7).
We also predicted an interaction between MHC in Tregs and LILRB1 in macrophages (Fig. 7g). The MHC receptor LILRB1 is a negative regulator of myeloid cell activation and a promotor of the M2 suppressive state. The MHC–LILRB1 interaction suppresses macrophages and is a target of cancer immunotherapy46. We first analyzed LILRB1 expression in macrophages by scRNA-seq and further validated it by FACS. We found that the expression of LILRB1 in macrophages increased in ESCC, compared to adjacent tissues (Fig. 7h, i and Supplementary Fig. 11d). Tregs expressed high levels of HLA-A, HLA-B, and HLA-C (Supplementary Fig. 11e, f). Multi-color IHC staining also showed potential physical interaction (co-localization) between LILRB1 expression macrophages and Tregs (Fig. 7j). When Tregs co-cultured with macrophages, we found that Tregs promoted macrophages expressing M2 markers, including CD163 and PDL1, and decreased TNFα expression. However, LILRB1 antibody blockade suppressed these affects (Fig. 7k, l and Supplementary Fig. 11g, h). These data suggested that Tregs may modulate macrophage function through HLA and LILRB1 interactions, and blocking this pathway may promote antitumor immunity in ESCC.
Ligand–receptor-based interaction between immune cells in ESCC. Circos plot showing the predicted average strength in seven patients of the indicated interaction between macrophage-provided ligands and receptors in other cell types within tumor or adjacent tissues.Violin plot showing the average expression value of a given gene in CD4-C6-FOXP3 from the tumor and adjacent tissues. The expression was measured as the log2 (count + 1), two-tailed Wilcoxon test.Expression levels of relative marker genes across NK-T cells illustrated in UMAP plots, the expression was measured as the log2 (count + 1).Detection of CD4, Foxp3, and IL1RII in ESCC tissue by multi-color IHC staining. Representative data from three patients was shown.,CD4+ CD25− (Teff) cells were stimulated with IL1β, Treg cells were co-cultured with Teff cells in the presence of anti-IL1R2 or IgG antibodies, the percentage of Teff dividing cells (= 6)and the IFNγ expression in Teff cells (= 4)was measured by FACS after 72 h culture.Circos plot showing the predicted interaction strength of the indicated L–R pairs between tumor and adjacent tissues as Treg-provided ligands and interactions with the indicated cell types averaged across seven patients.Violin plot showing the average LILRB1 expression in macrophages from tumor and adjacent tissues. The expression was measured as the log2 (count + 1), two-tailed Wilcoxon test.FACS analysis of LILRB1 expression in CD68+ macrophage cells from adjacent and ESCC tumor tissues (= 3).Detection of CD4, Foxp3, CD68, and LILRB1 in ESCC tissue by multi-color IHC staining. Representative data from three patients was shown.,Macrophages were co-cultured with or without Treg cells, and flow cytometry measured the expression of CD163, CD206, and PDL1. Representative data from four patients was shownand TNFα expression. Representative data from three patients was shown (= 3)in CD68 macrophage cells in the presence of anti-LILRB1 or control IgG antibody. a b c d e f e f g h i j k l k l n n n n
Discussion
Here, we present a comprehensive characterization of immune cells in seven pairs ESCC tumors and matched adjacent tissues. Immune signature profiling and TCR β-chain repertoire analysis have been studied in ESCC using mRNA microarray and bulk RNA-seq, respectively47,48. Recently, a mouse model mimicking human ESCC development and construction of a single-cell ESCC developmental atlas have been reported49. Here, we combined deep sc-RNA-seq and TCR-seq, and illustrated the whole immune landscape, including the innate and adaptive immune cell atlas in ESCC and adjacent tissue. Our work will lay the foundation for developing and applying immune-targeted strategies for ESCC diagnosis and treatment.
Our study showed that ESCC was enriched in immune-suppressive cell populations, including Tregs, exhausted CD8 T, CD4 T and NK cells, M2 macrophages, and tDCs. All these immune-inhibitory cells may contribute to immune escape and promote tumor progression. Interestingly, we demonstrated that exhausted CD4, CD8 T cells, and NK cells were the major intra-tumoral proliferating immune cell compartments, although these cells were enriched in exhaustion genes. Pre-exhausted clusters CD8-C5-CCL5 and CD8-C6-STMN1 may serve as better targets for immunotherapies compared to the exhausted cluster (CD8-C7-TIGIT), as the latter are in a permanent and less reversible exhausted stage, making them more resistant to checkpoint inhibition due to their epigenetic changes50,51. Furthermore, we found that tumor-infiltrating NK cells were not only commonly reduced in ESCC, but also expressed high levels of checkpoint molecules, including NKG2A and CD49d, suggesting an exhausted state. Anti-NKG2A and anti-CD49d have been reported to be checkpoint inhibitors that promote anti-tumor immunity52,53. Our results suggested alternative pathways to re-activate anti-tumor immunity in ESCC, such as the blockade of NKG2A and CD49d alone or in combination with anti-PD1/PD-L1, which may improve the immunotherapeutic response.
Clonal amplification of cells carrying identical TCR-seq across clusters was a strong evidence to support the connection among these clusters and indicated the transition of cell status. While most cells contained unique TCRs, clonal amplification was observed to varying degrees in different clusters. Indeed, we found sharing of TCR sequences in almost all CD4 and CD8 clusters, suggesting a broad differentiation after T cell priming. Notably, exhausted CD8 T cells (CD8-C7-TIGIT) carried a much higher percentage of shared clones with other CD8 clusters, especially with the pre-exhaustion clusters CD8-C5-CCL5 and CD8-C6-STMN1, which was consistent with the related status of these clusters and the multi-step exhaustion hypothesis. In contrast, CD8-C1-NKG7, the most cytotoxic CD8 T cell cluster, had significantly fewer clonal T cells in tumors. This result suggested the inhibition of cellular proliferation of cytotoxic CD8 T cells by TEM. This may contribute to the decreased numbers of tumor-infiltrating effector CD8 T cells and to the immune-suppressive microenvironment. A similar finding was recently reported in bladder cancer, in which cytotoxic CD8 T cells (FGFBP2+ cluster in their analysis) were more clonal in normal tissues than in tumors54. On the other hand, it is possible that the clonal amplification of esophageal cytotoxic cells was due to occasional exposure to non-tumor antigens in the esophagus, and that most CD8-C1-NKG7 cells in adjacent tissues express low levels of CD39, indicating that they are bystander CD8 T cells. Interestingly, we did not find a significant population of naïve CD8 T cells in our samples, and similar studies in other cancer types, e.g. melanoma and liver carcinoma, identified a limited number of naïve CD8 T cells11,55. It is possible that naïve CD8 T cells less frequently infiltrate into the esophagus or are activated by the local environment, making them less likely to be detected as an independent population. Additional studies are needed to validate this phenomenon and identify the mechanism.
The critical role of monocytes/macrophages in tumors has been described in liver, breast, and lung cancers using scRNA-seq. Typically, macrophage activation is classified into either a pro-inflammatory M1 state or an M2 state associated with the resolution of inflammation56. Our analyses revealed that monocytes/macrophages reside along a spectrum of monocyte, M1, and M2 states; both M1-associated and M2-associated genes were frequently co-expressed in the same cells. Zhang, et al. reported six macrophage clusters identified in hepatocellular carcinoma, and found that a macrophage cluster (M4-C1-THBS1 in their analysis) was enriched for signatures of MDSC, similar to ours MDSC-C1-C1QC and MDSC-C2-APOE; they also found another cluster M4-C2-C1QA, co-existed with M1 and M2 signatures20. Lambrechts, et al. reported that macrophages in lung TME show rheostatic phenotypes and become M2 polarized19. Indeed, we found the M2 cluster Macro-C3-CSF1 was enriched in ESCC compared with adjacent tissue. The co-existence of M1 and M2 signatures indicated that TAMs were more complex than the classical M1/M2 model, and this phenomenon was also found in breast cancer and liver carcinoma12,18. When using WGCNA to analyze gene correlation from monocytes/macrophages, we found a gene set that was positively related to monocytes, patients who expressed higher of this signature had significantly better prognosis in ESCC and other tumor types with similar pathology, suggesting that this signature may serve as a prognostic biomarker in ESCC and other squamous cell carcinomas.
cDC2 is specifically involved in MHC class II-mediated antigen presentation and the activation and expansion of CD4 T cells. cDC1 is also necessary for anti-tumor immunity. Recently, some scRNA-seq papers have reported DCs in the TME12,40,57. We found that most of them were consistent with our data. Zilionise et al. reported that human DCs contain four distinct subsets, including hDC1, hDC2, hDC3 and pDC, hDC1, hDC2, and hDC3 were consistent with our cDC1 (DC-C1-CLEC9A), cDC2 (DC-C2-CLEC10A), and Lamp3+DC (DC-C3-LAMP3) subsets, respectively. However, the monocyte-derived DC (DC-C4-FCER1A) was absent in that paper. These differences may reflect a variation in DC states between different tumor tissues and/or the setting of analysis. Interestingly, we found that LAMP3 + DC had multiple functions, such as activation activity, migration activity, and tolerogenic ability, and this subset was also found and reported in lung and liver cancers, suggesting the conserved myeloid cells exist in many tumors.
Our work further demonstrated the interaction between Tregs and macrophages which may contribute to the immunosuppressive microenvironment in ESCC. Our data suggested that IL1R2 expressed on Tregs may enhance Treg function by blocking IL1β-dependent effector T cell activation. IL1R2, which serves as the IL1β decoy receptor and binds to IL1β can block follicular cell activation58. It has been reported that IL1R2 expressed on activated tumor Tregs and is correlated with poor prognosis in lung adenocarcinoma59,60. Our study also demonstrated that Tregs may modulate macrophage function through HLA-A, B, C, and LILRB1 interaction. LILRB1 expression was upregulated on TAMs, disruption of either MHC class I or LILRB1 potentiated phagocytosis of tumor cells both in vitro and in vivo. Investigating the mechanisms underlying this immunosuppressive MHC class I-LILRB1 signaling axis in TAMs will be useful in developing therapies to restore macrophage function46,61, and blocking this pathway may promote antitumor immunity in ESCC.
In summary, our transcriptional map of immune cells from ESCC and adjacent tissue provided a framework for understanding the immune status and revealed the dynamic nature of immune cells in the ESCC setting. In addition, we illustrated the immune-suppressive state of ESCC from many aspects, all of which are potential novel targets for developing immunotherapies in ESCC and other cancers.
Methods
Sample collection and preparation
Seven patients who were pathologically diagnosed with ESCC were enrolled in this study for single cell RNA-seq analysis. None of the patients had been treated with chemotherapy, radiation, or any other anti-tumor medicines prior to tumor resection. The clinical information of these patients was summarized in Supplementary Table. Paired, freshly excised ESCC tumors and adjacent esophageal tissues were obtained immediately after surgical removal. The adjacent tissues were at least 5 cm from the tumor tissues. Clinical samples were collected from the Xinhua and Ruijin Hospital, Shanghai Jiaotong University School of Medicine. Prior to participation, written informed consent was obtained from all subjects. All studies were performed in accordance with the Declaration of Helsinki. The study was approved by the Research Ethics Board of the Xinhua and Ruijin Hospitals, Shanghai Jiao Tong University School of Medicine. 1
Tumor and adjacent tissues were isolated by mincing the freshly obtained surgical specimens into 1-mm cubic pieces, followed by enzymatic digestion using 0.1% collagenase IV, 0.002% DNAse I, and 0.01% hyaluronidase, and were incubated on a rocker for 20–40 min at 37 °C. The digested tissues were then passed through a 40 µm cell strainer and washed twice with PBS prior to surface staining. Immune cells were stained at 1 × 106 cells per ml with antibodies (CD45-APC, CD235-FITC) for 30 min at 4 °C, and then washed and suspended in 200 µL FACS buffer. PI was added 5 min before flow cytometry sorting for dead cell discrimination. CD45+ CD235−PI−cells were sorted with FACS Aria II Cell Sorter (BD Biosciences); post-sort purity was routinely >95% for the sorted populations. The CD45+CD235− PI− cell gating strategy is shown in Supplementary Fig. 12a.
RNA-Seq library preparation for 10× Genomics single-cell 5′ and VDJ sequencing
FACS-sorted CD45+ cells were suspended at a 1 × 106 cells/ml concentration in FACS buffer and the viability was higher than 85%. Single-cell library preparation was carried out according to the manual of the Chromium Single-cell V (D) J Reagent Kits. Briefly, cell suspensions were first loaded on the Single-cell A Chip for gel bead-in-emulsion (GEM) generation and barcoding, aiming for a recovery of 10,000 cells per sample. Next, reverse transcription, RT cleanup cDNA amplification and quantification were performed to build a 5′ gene library. In addition, a human T Cell V (D) J Enrichment Kit was used to isolate and enrich for the V (D) J sequence before TCR library construction. Finally, libraries were sequenced on the NovaSeq 6000 system with NovaSeq 6000 S4 Reagent Kit (300 cycles).
scRNA-seq data analysis
The 10× Genomics Cell Ranger (3.0.1 version) pipeline was used to demultiplex raw files into FASTQ files, extract barcodes and UMI, filter, and map reads to the GRCh38 reference genome, and generate a matrix containing normalized gene counts versus cells per sample. This output was then imported into the Seurat (v3) R toolkit for quality control and downstream analysis. All functions were run with default parameters, unless otherwise specified. Low-quality cells (<400 genes/cell and >10% mitochondrial genes) were excluded. As a result, 80,787 cells with a median of 1170 detected genes per cell were included in downstream analyses. To remove the batch effect, the datasets collected from different samples were integrated using Seurat v3 with default parameters.
Dimensionality reduction, clustering, and annotation
We then identified a subset of genes that exhibit high cell-to-cell variation in the dataset, which helped to represent the biological signal in downstream analyses. The Seurat function ‘FindVariableFeatures’ was applied to identify the highly variable genes (HVGs). The top 2000 HVGs were used for data integration. The data were scaled using ‘ScaleData’ and the first 20 principle components were adopted for auto-clustering analyses using ‘FindNeighbors’ and ‘FindClusters’ functions. For all 80,787 cells, we identified clusters setting the resolution parameter as 1.5, and the clustering results were visualized with the UMAP scatter plot. The marker genes of each cell cluster were identified using the ROC analysis function provided by the Seurat ‘FindAllMarkers’ function for the top genes with the largest AUC (area under curve). We used GSEA (http://software.broadinstitute.org/gsea/msigdb/annotate.jsp↗) to perform biological process enrichment analysis on the differentially expressed genes in each cluster or subset. The whole dataset was then categorized into NK cells, T cells, monocytes/macrophage cells, B cells, mDCs, pDCs, plasma cells, mast cells, and other cells (including fibroblast cells and basal cells) according to the known markers: KLRC1, KLRD1 (NK cells), CD3G, CD3D, CD3E, CD2 (T cells), CD14, VCAN, FCGR2A, CSF1R (monocyte/macrophages), CD19, CD79A (B cells), CD1C, FCER1A (mDC), CLEC4C (pDC), SLAMF7, IGKC (plasma cells), TPSB2, CPA3 (mast cells), KRT19, IGFBP4, CTSB (basal cells/fibroblasts). Clusters were also confirmed by identifying significantly highly expressed marker genes in each cluster and then comparing them with the known cell-type-specific marker genes. The 44,634 NK-T cells were further categorized into 17 clusters with a resolution parameter of 1. The 17,256 myeloid cells were further categorized into 14 clusters with a resolution parameter of 0.4.
Definition of exhaustion, naivete, cytotoxicity, and other scores
We considered published signature gene lists for cytotoxicity, exhaustion, Treg, M1, M2, monocyte, and MDSC (Supplementary Data 3), and the DC cell activation, migration and tolerogenic scores have been previously described62 and are provided in Supplementary Data 5. We used the Z-scores of the average expression of four well-defined naïve markers (CCR7, TCF7, LEF1, and SELL) to define the naïve gene score for both CD8 and CD4 T cells.
Development of trajectory inference
We applied Monocle (version 2) and Slingshot using default parameters on the identified monocyte and macrophage cells to determine the potential lineage differentiation. First, the top 1000 HVGs were selected using the function ‘FindVariableFeatures’ in Seurat v3. Then the ‘estimateSizeFactors’, ‘estimateDispersions’, and ‘dispersionTable’ functions with default parameters were used to build statistical models to characterize the data. The genes with mean expression >0.1 were retained for subsequent analyses. The dimensionality reduction was performed by the ‘reduceDimension’ function using the DDRTree method.
SCENIC analysis
The SCENIC analysis was run as described38, using the pyscenic (version 0.9.19) and hg19-500bp-upstream-10species databases for RcisTarget, GRNboost, and AUCell. The input matrix was the normalized expression matrix that was from Seurat.
Cell–cell interaction
Cell–cell interaction analysis was analyzed by two algorithms, one has been described previously43, and the other was scTHI (https://github.com/miccec/scTHI↗). We determined the existence of the potential interactions between two cell types mediated by evaluating the expression levels of annotated ligands and receptors between the two cell types. A ligand–receptor score for each sample was calculated by multiplying the average expression of a ligand from one cell type and the average expression of a receptor from the other cell type. The interactions between the two cell types via the ligand–receptor pair were identified when the scores of the samples from the seven patients were significantly larger than zero (Wilcox test P < 0.05).
TCR data analysis
A 10× Genomics Cell Ranger pipeline was used to identify clonotypes by alignment and annotation with the default settings, according to the manufacturer’s recommendations. TCR reads were aligned to the GRCh38 reference genome, and consensus TCR annotation was performed using Cell Ranger VDJ. Only in-frame TCR alpha–beta pairs were considered to define the dominant TCR of a single-cell. Each unique dominant alpha–beta pair was defined as a clonotype. For two cells to be assigned to the same clonotype, both alpha and beta sequences had to be shared. If one clonotype was present in at least two cells, cells harboring this clonotype would be considered clonal, and the number of cells with such dominant alpha–beta pair indicated the degree of clonality of the clonotype.
Survival analysis
An online website (http://kmplot.com/analysis/index.php?p=service↗)63 was used to analyze the relationship between BHLHE40, and the 50-gene signature and ESCC survival.
Flow cytometric analysis
For cell surface marker staining, antibody cocktails were used (CD45, PI, CD3,CD4, CD8, CD19, CD138, CD25, CD127, CCR7, CD83, PD1, CD68, CD163, CD206, CD11C, CD14, CD56, NKG2A, CD66B, CD15, CD11B, CD117, HLA-DR, PDL1, CD11C, CD121b, LILRB1) For intracellular cytokine staining, cells were stimulated with a cell stimulation cocktail plus protein transport inhibitors (eBioscience) for 5 h. Then, the cells were fixed and permeabilized with Cytofix/Cytoperm buffer, and intracellular cytokines were stained with antibodies against FOXP3, IFNγ, TNFα, IDO1, LAMP3, and isotype control according to the manufacturer’s instructions. The antibody dilutions are in accordance with the instructions. Flow cytometric analysis was performed with a FACS Canto II instrument (BD Bioscience) and FlowJo software (TreeStar). Figure 3c, f, i gating strategies were shown in Supplementary Fig. 12b–d. Figure 6e gating strategies were shown in Supplementary Fig. 12e. Figure 6f gating strategies were shown in Supplementary Fig. 12f. Figure 6h gating strategy was shown in Supplementary Fig. 12g. Figure 6i gating strategy was shown in Supplementary Fig. 12h. Figure 7e, f gating strategies were shown in Supplementary Fig. 12i. Supplementary Figs. 7i and 9e gating strategies were shown in Supplementary Fig. 12j. Supplementary Fig. 11a, f gating strategies were shown in Supplementary Fig. 12k. Figure 7k, l gating strategies were shown in Supplementary Fig. 12l.
In vitro stimulation of DC subsets
Fresh blood was obtained from 30 healthy volunteers, (age range from 34 to 62, median age is 46). Peripheral blood mononuclear cells (PBMCs) were isolated from heparinized fresh blood by standard density gradient centrifugation with Ficoll-Paque Plus (GE Healthcare). Subsequently, Pan-DCs were obtained by negative selection using a Human Pan-DC pre-enrichment kit (Stem Cell Technologies). DCs were cultured in CellGenix GMP DC medium (CellGenix), and stimulated with LPS-EB Ultrapure (Sigma-Aldrich) and INFγ (Peprtotech); both at 100 ng/ml or not for 24 h, and then co-cultured with CD4+CD45RA+ naïve T cells in the anti-CD3 (eBioscience) 1 μg/ml for 4 days, and FOXP3 expression was detected by flow cytometry.
Macrophage differentiation
To obtain macrophage cells from healthy PBMC, the cells of each subject were cultured by using RPMI 1640 (Gibco) medium that included 10% fetal bovine serum (FBS) (Gibco), 100 U/ml penicillin, 100 μg/ml streptomycin (Gibco) and 1% 2 mM l-glutamate (Gibco) at 37 °C in a %5 CO2 incubator for 4 h. After the culture period, non-adherent cells were washed with phosphate-buffered saline (PBS) and cultured in different culture plates by using RPMI 1640 complete medium and 10 ng/ml human IL-2 (Peprotech) until to co-culture experiments. Adherent cells were detached by trypsin-EDTA (Gibco) and divided into 1 × 105 cells per well in two 12-well plates. Cells were cultured using RPMI 1640 complete medium and human 25 ng/ml macrophage colony-stimulating factor (M-CSF) (Peprotech) at 37 °C in a 5% CO2 incubator for 4 days.
IL1R2 antibody blocking
FACS sorted CD4+CD25− Teff cells were stained with CFSE (Sigma-Aldrich) at 2 μM for 8 min at room temperature. A total of 3 × 104 Teff cells were stimulated with or without IL-1β at 20 ng/ml, then co-cultured with or without FACS-sorted CD4+CD25hiCD127lo/- Treg cells at 3 × 104 in the presence of anti-CD3 1 μg/ml, anti-CD28 at 0.5 μg/ml (eBioscience), stimulated with anti-IL1R2 or control IgG (R&D) at 20 μg/ml, and then cultured for an addition 72 h. The expression of IFNγ and the division of Teff cells were analyzed by flow cytometry.
LILRB1 antibody blocking
Macrophages (3 × 104) were co-cultureed with FACS sorted CD4+CD25hiCD127lo/- Treg cells (3 × 104) in the presence of anti-CD3 1 μg/ml, stimulated with anti-LILRB1 or IgG control (R&D) at 10 μg/ml, and cultured for an additional 48 h. The expression of CD163, CD206, PDL1, and TNFα in macrophage cells was analyzed by flow cytometry.
Multi-color IHC
Human tissue specimens were provided by Xinhua and Ruijin Hospital, Shanghai Jiaotong University School of Medicine under an approved Institutional Review Board protocol. The specimens were collected within 30 min after tumor resection and fixed in formalin for 48 h. Dehydration and embedding in paraffin were performed following routine methods. These paraffin blocks were cut into 4 mm slides and fixed to glass slides. Then the paraffin sections were placed in the 70 °C paraffin oven for 1 h before deparaffinization in xylene and then rehydrated in 100%, 90%, 70% alcohol successively. Antigen was retrieved by immersion in boiling EDTA buffer (pH 9.0) for 15 min. Following a preincubation with Antibody Diluent/Block (PerkinElmer) to block nonspecific sites for 15 min, the sections were incubated with primary antibodies in a humidified chamber at 4 °C overnight. After the sections were washed with PBS twice for 5 min, they were incubated with Opal Polymer HRP Ms+Rb (PerkinElmer) for 10 min at 37 °C. The antigenic-binding sites were visualized using the Opal Seven-Color IHC Kit (PerkinElmer) according to the manufacturer’s protocol. Images were captured and analysis was performed with PerkinElmer Vectra 3.0.5 (PerkinElmer) and informs (PerkinElmer), respectively.
The primary antibodies, IHC metrics and fluorophores used in the validation of tolerogenic LAMP3+DC, rabbit anti-human CD11C (Abcam, 1:200) (Oapl690), rabbit anti-human LAMP3 (Abcam, 1:100) (Opal620), rabbit anti-human IDO1 (Abcam, 1:100) (Opal520), and rabbit anti-human PDL1 (Abcam, 1:100) (Oapl570). The validation of the potential physical interaction (co-localization) between macrophages and Treg cells were: rabbit anti-human CD4 (Abcam, 1:500) (Oapl690), mouse anti-human FOXP3 (Abcam, 1:100) (Oapl650), mouse anti-human CD68 (Abcam, 1:100) (Opal620), and rabbit anti-human LILRB1 (Abcam, 1:500) (Opal520). IHC was performed using a fully automated IHC instrument BOND MAX (LEICA). The primary antibodies used were anti-human CD3 (LEICA, 1:200), anti-human CD20 (LEICA, 1:300), anti-human MPO (Long island, 1:100), anti-human CD56 (Long island, 1:400), and anti-human CD68 (Long island, 1:100). All the antibodies information is included in Supplementary Table. 2
Weighted gene co-expression network analysis (WGCNA)
To construct the gene co-expression network, the R package WGCNA was implemented and a signed network was constructed using genes that were selected as noise robust by OGFSC using default parameters64. To group genes with coherent expression profiles into modules, average linkage hierarchical clustering was applied, using the topological overlap measure as the dissimilarity. The modular gene centrality, defined as the sum of within-cluster connectivity measures, was used to rank modular genes for hubness within each gene module.
Statistical analysis
The statistical methods used for each analysis are described in the above “Methods” sections and in the figure legends.
Reporting summary
Further information on research design is available in thelinked to this article. Nature Research Reporting Summary