https://orcid.org/0000-0002-2650-8051
Figures
Abstract
Purpose
Merkel cell carcinoma (MCC) is a neuroendocrine carcinoma originating in the skin. Studies are needed to determine the mechanisms of immune escape in patients with MCC, and malignant cell conditions that promote immune evasion.
Methods
We used Single-cell RNA sequencing (scRNA-seq) to determine cellular features associated with MCC disease trajectory. A longitudinal multi-omics study was performed using scRNA-seq data of peripheral blood harvested from four-time points. Six major cell types and fifteen cell subgroups were identified and confirmed their presence by expression of characteristic markers. The expression patterns and specific changes of different cells at different time points were investigated. Subsequently, bulk RNA data was used to validate key findings.
Results
The dynamic characteristics of the cells were identified during the critical period between benign improvement and acquisition of resistance. Combined with the results of the validation cohort, the resistance program expressed in the relapse stage is mainly associated with T cell exhaustion and immune cell crosstalk disorder. Coinciding with immune escape, we also identified a decrease non-classical monocytes and an expansion of classical monocytes with features of high inflammation and immune deficiency.
Conclusion
Changes in cellular status, such as depletion of T cells and dysregulation of B cell proliferation and differentiation, may lead to drug resistance in MCC patients. Meanwhile, the widespread decreased antigen presentation ability and immune disorders caused by deletion of MHC class II gene expression should not be ignored.
Citation: Tao Q, Du J-x, Zhang S, Lin W, Luo Y, Liu Y, et al. (2023) Longitudinal multi-functional analysis identified responses of T cells, B cells, and monocytes as hallmarks of immunotherapy tolerance in patients with merkel cell carcinoma. PLoS ONE 18(11): e0293922. https://doi.org/10.1371/journal.pone.0293922
Editor: Yoon-Seok Chung, Korea Disease Control and Prevention Agency, REPUBLIC OF KOREA
Received: April 13, 2023; Accepted: October 21, 2023; Published: November 20, 2023
Copyright: © 2023 Tao et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript.
Funding: This study was funded by the 2021 Special Fund for Science and Technology Innovation Strategy of Guangdong Province [grant number pdjh2021b0120]; 2020 National Innovation and Entrepreneurship Program for Chinese College Students [grant number 202010572017]; and the 2020 Guangzhou University of Chinese Medicine College Student Innovation and Entrepreneurship Project [grant number 202010572201]. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Merkel cell carcinoma (MCC) is a rare and aggressive neuroendocrine carcinoma originating in the skin, usually caused by Merkel cell polyomavirus (MCPyV) [1–3]. Due to the characteristics of high misdiagnosis rate, poor prognosis, easy local recurrence, and metastasis [4], MCC has a 5-year survival rate of 40% and a fatality rate of 33%, which is even higher than melanoma [1]. Standard MCC therapies are mainly divided into the following two types: Primary and metastatic lymph node tumors are treated with surgery and radiotherapy as the mainstay. Unresectable and metastatic tumors were treated with cytotoxic chemotherapy and immune checkpoint inhibitors (ICIs) as the mainstay [5]. Nghiemet et al. reported a 56% response rate in patients with advanced MCC, receiving anti-PD1 antibody pembrolizumab to block the immune checkpoint [6]. The 6-month progression-free survival rate was 67% for immunotherapy, compared with 24% for chemotherapy [7,8]. These studies showed that immunotherapy had good curative effect and great potential for patients with advanced MCC. However, drug tolerance and immune escape associated with immunotherapy have become key issues in the clinical use of ICIs.
ICIs have changed the therapeutic landscape of several types of cancer [9], especially MCC [10]. Nonetheless, many patients with MCC acquired resistance after immunotherapy, which is often intrinsic [11–13]. The mechanism of ICIs resistance in MCC was unknown. ICIs target cell-cell interactions, so resistance may derive from different cells and their interactions in the tumor ecological system. Among several possible mechanisms for the acquired resistance, immune damage caused by loss of expression of major histocompatibility complex (MHC) molecules was the most prevalent [14–16]. T cell infiltration level was also thought to be associated with acquired resistance [17], but the complex mechanisms involved in drug resistance obviously more than that. Among many factors, the biological relationship between cell subsets and different temporal pathological stages may be important for us to elucidate the mechanisms of drug resistance and recurrence.
Compared with conventional bulk sequencing, single-cell RNA sequencing (scRNA-seq) is more outstanding in profiling cancer immune interactions and producing reliable response biomarkers [18]. In this study, different cell clusters in the peripheral blood mononuclear cells (PBMCs) of a patient with drug-resistant MCC were identified. We explored the proliferation potential and genetic relationships of T cells at the levels of differentiation and evolution. By studying the expression patterns of genes in PBMCs, we found a series of highly expressed genes associated with drug resistance and immune avoidance pathways. Meanwhile, these genes may serve as cellular markers. Moreover, we further revealed the evolutionary routes of various cellular from recovery to relapse in the same patient. Dissecting the immune checkpoint dysregulation and immune escape mechanisms occurring in the immunotherapy of this patient, is beneficial for understanding the reasons behind the immunotherapy failure and cancer recurrence in MCC. This, in turn, further promotes the development of targeted treatment approaches for MCC.
Materials and methods
Single-cell RNA-seq data acquisition and processing
ScRNA-seq data (GSE117988) has been described by Paulson KG et al [19]. Patients with advanced MCC were treated with autologous ex vivo expanded CD8+ T cells, and experienced an initial reduction in overall tumor pressure but eventually relapsed. Drug resistance occurred during the 22-month clinical response period, and PBMCs samples were collected at four time points for 10X single-cell sequencing: pre-treatment (Pre), early post treatment day + 27 (D+27), responding post treatment day + 376 (RespD376), relapse/acquired resistance post treatment day + 614 (ARD614).
We performed simple filtering of the expression matrix by quality control, dimension reduction, and feature selection with the Seurat package (R version 4.0.2) [20]. Cells expressing more than 500 genes and genes expressing more than 10% of the cells were selected for further analyses. Cells with detected gene numbers < 500 or > 6,000 and mitochondrial unique molecular identifiers (UMIs) > 10% were removed. Then the filtered expression matrix was normalized using NormalizeData function, and the top 2,000 genes with a high degree of intercellular variation were calculated using FindVariableFeatures function. RunPCA function was used to reduce the dimensionality of the data followed by nonlinear dimensional reduction.
Cell type annotation and cluster marker identification
Dimensionality reduction was performed for the gene-barcode matrix using PCA and UMAP analyses. The first 10 principal components (PCs) were retained and selected for nonlinear dimensional reduction. These PCs were projected into two-dimensional space using UMAP. To cluster the cells, share-nearest neighbor (SNN) graph method was implemented in Seurat [20]. Subsequently, the distinct clusters of cells were divided according to a neighborhood size of 40 and a resolution parameter of 0.6. To optimize feature selection, three different cell clusters were excluded because of likely blood contamination. Cluster-specific markers were identified via differential expression using FindMarkers function in Seurat. Clusters were then annotated and classified on the basis of the expression of typical markers of specific cell types.
Differential expression genes analysis
In order to further explore the mechanism of tumor recurrence, we conducted a longitudinal differential expression analysis of cells at RespD376 and ARD614 time points. Differentially expressed genes (DEGs) between the two above time points were identified by applying FindMarkers function in Seurat. DEGs in each cell type were identified with false discovery rate (FDR) adjusted p-value of less than 0.05.
Co-expression analysis
To reflect the expression regulatory relationship among DEGs, the WGCNA package (version 1.70.3) was used to identify the co-expressed modules of DEGs. All the DEGs identified from the above longitudinal analysis (981 genes) were used to obtain the gene co-expression modules. First, we calculated pairwise gene correlations based on the log-transformed normalized expression counts matrix across all cells. A soft threshold function was used to construct a signed adjacency matrix, which was the basis for constructing the topology overlap matrix (TOM) for the next step. The TOM was used to construct a gene tree by hierarchical clustering. DEGs were then split into different modules according to the gene tree by using CutreeDynamic function, while the minimum of the module size was set to 15. Then mergeCloseModules function with a parameter cut height set to 0.45 was used to merge modules that were closely associated. Moreover, the module eigengene values for each time point were calculated.
Cell-Cell communication analysis
The receptor-ligand interactions that exist in different cell types were confirmed based on the significant DEGs that have been identified in these cell clusters. The reference data set for these ligand-receptor pairs come from the R package iTALK [21]. If the average expression of a ligand or receptor in a certain cell type is above the threshold of 0.5 (log2 (average expression)>0.5) and the ligand or receptor is expressed more than 10% of the cells, the cell type is defined as "expressing" the ligand or receptor transcript. After the expression of the receptor or ligand is determined, the corresponding interaction will be defined as incoming or outgoing of that cell type, respectively.
Subclustering of major interested cell types
For the cell types identified as important or of interest by the previous analysis, cells were first extracted from the integrated data list. Next, these major cell types were divided into further subclustering by linear and nonlinear dimensionality reduction methods.
Exhaustion score calculation
In order to assess the level of exhaustion of cell subtypes and identify cell subtypes with exhaustion features, we utilize the AddModuleScore function in Seurat to perform a characteristic scoring of cells based on the expression of TIGIT, PDCD1, CTLA4, CD27, LAG3.
Bulk RNA data validation
A total of 124 samples were included in the validation dataset. These samples were collected from GSE39612 cohort of GEO database, which included tumor tissue from 60 MCC patients and 64 normal tissues as controls. Deconvolution algorithm Cibersort was used to quantify tumor-infiltrating immune cells from second-generation sequencing data. The expression matrix of top 50 genes specifically expressed in each cell type was extracted as the reference data set according to the result of cell cluster annotation.
Functional enrichment analysis
Gene set variation analysis (GSVA) was executed using R package GSVA (version 1.38.2) [22] to quantify the pathway enrichment score across time-points or cell types. Metabolic pathway activities were evaluated using a curated dataset that contained 85 metabolic pathways [23]. Limma R package (version 3.38.3) was used to calculate differential activities of pathways between time-points or cell types. Pathway with a corrected P ≤ 0.01 was identified significantly disturbed.
Results
MCC patients exhibit acquired resistance that depends upon changes in the distribution of cell clusters in different stages
We projected the MCC scRNA-seq data onto a UMAP figure and resolved the cell populations characterized by transcriptional signatures and stages of MCC patients respectively (Fig 1A and 1B). Six cell phenotypes were identified: T cell, B cell, NK cell, monocyte, erythroblast, and platelets. The expression of common markers in these cell types were shown in Fig 1C and 1D. Relative percentages of the six cell phenotypes were calculated relative to the stages of patients with MCC (Fig 1E and 1G). There was an increase in the percentage of lymphocytes, including T cell, B cell, and monocytes (Fig 1G). Cell-cell communication among six cell phenotypes were presented in the interaction plot (Fig 1H). The relationship between ligands and receptors was shown in the inner circle. T cells played an important role in cell-cell communication, followed by monocytes and B cells.
(A) Cellular populations identified. The UMAP projection of 12874 single cells from MCC samples shows the formation of 6 main clusters with label names. Each dot corresponds to a single cell, colored according to cell type. (B) Pseudotime UMAP representation of all merged cells, colored by time point. (C) Canonical cell markers were used to label clusters by cell identity as represented in the UMAP plot. (D) Dot plot for cell-type-specific signature genes. Color discriminates genes with increased (red) or decreased (blue) expression, and point size represents the number of cells per group expressing the corresponding gene. (E) UMAP plot of cluster breakdown by time point. (F) The average proportion of six main cells among different time points. (G) Percentages of the six types among four groups. Groups are shown in different colors. (H) Circos plot showing the interactions among different cell types in MCC.
The heatmap exhibited some of the genes characteristically expressed in different cell types at different times (Fig 2A). To capture the determinants that drive MCC from a benign response to resistance, DEGs in the stages of RespD376 and ARD614 were further analyzed (Fig 2B). In other words, we divided six kinds of cells into group RespD376 and group ARD614 for difference analysis. The number of DEGs obtained in each cell type were illustrated in Fig 2C (A total of 981 DEGs). Subsequently, 981 DEGs identified from RespD376 and ARD614 were used for weighted gene co-expression network analysis. Five modules were identified in total, which were defined as M1–M5, of co-expressed genes following specific expression patterns throughout the MCC disease phases (Fig 2D). For example, M4 module genes were mainly involved in monocyte differentiation, cytokine-mediated signaling pathway, and lymphocyte-mediated immune response (Fig 2E).
(A) Heatmap showing the expression of the top 18 differentially expressed genes for each cell cluster (q < 0.05). (B) Differential genes (DEGs) of six major cell types between groups RespD376 and ARD614. Colors discriminate DEGs with upregulation (red) or downregulation (blue) expression. The Y-axis represents the log2-fold change of DEGs. (C) Number of DEGs between RespD376 and ARD614 in each cell type. Colors discriminate upregulation (red) or downregulation (blue) expression. (D) Heatmap showed the average expression of module hub genes in different cell types and time points. (E) The dot plot showed the pathways enriched in specificity in the five modules.
A proliferative exhausted CD8 + T cell subpopulation emerges in the late relapse stage of MCC patient
Single-cell transcriptome data of T cells were extracted for re-clustering, and four subpopulations with transcriptional signatures were resolved (Fig 3A). First, the high expression of GAPDH, which indicated that glycolysis and exhaustion overexpressed in cell cluster 3 attracted our attention (Fig 3B). Further, investigation confirmed that the exhaustion-related markers LAG3, SRGN and CD27 were all upregulated in cluster 3 (Fig 3C). Besides, both CD8A and CD3E were highly expressed, so cluster 3 was identified as an exhausted CD8+ T cell. Cytotoxic and proliferative phenotypes were also identified. For instance, the cytotoxicity-associated markers GNLY, CSTW and CST7 were up-regulated in cluster 2, and cluster 3 expressed proliferation markers MKI67 and TYMS counterintuitively (Fig 3D). Transcript levels of exhaustion marker GAPDH showed a positive correlation with MKI67 for cluster 3 cells (Fig 3E). That may reflect the proliferative hierarchy for exhausted T cells maintenance in the context of chronic infections and cancer [24]. Quantitative analysis was performed to verify the existence of exhausted cells and elevated exhaustion score were observed for T cell cluster 3 (Fig 3F). Meanwhile, the developmental trajectory of the four clusters presented a binary branched structure (Fig 3G), cluster 2 cells as the root, and exhausted CD8+ T cells as the end state, which was consistent with previous studies [25,26]. In this process, T cell exhaustion-related immune checkpoints (LAG3, CTLA-4, PD-1 and TIGIT) tended to be up-regulated, while synergist checkpoints (KLRG1) tended to be downregulated (Fig 3H). In addition, expressions of LAG3, TIGIT and proliferation markers TYMS and MKI67 in MCC were also found to be higher than those in normal samples in the validation cohort (Fig 3I). Tumor immune invasion analysis also showed an increased estimated proportion of this cluster (C3) cells in MCC patients (Fig 3J). The cells in cluster 3 were exhausted CD8 + T cell subpopulation in a state of proliferation, and T cell phenotypic composition evolved with MCC severity.
(A) T cell compartment subtypes are represented as a UMAP. (B) Expression of GAPDH in four T cell subtypes. (C) Expression of the exhaustion-related markers. (D) Expression levels of various subtype markers in T cells. Red dots indicate high expression and blue dots indicate low expression. (E) Correlation between proliferation marker expression and GAPDH expression in C3 cluster cells. (F) Exhaustion score calculated based on expression levels of exhaustion markers. “**” indicates that the T-test p-value is less than 0.05. (G) Differentiation trajectory of CD8+ T cells in MCC, with each color-coded for clusters (top) and pseudotime (bottom). (H) Expression of the common immune checkpoints. Red indicates high expression and blue indicates low expression. (I) The expression of some genes in the validation cohort. (J) The predicted proportion of C3 cluster cells in the validation cohort. (K) Number of DEGs between RespD376 and ARD614 in CD8+ T cells. (L) Venn diagrams represent 82 DEGs that were up-regulated in RespD376 and down-regulated in ARD614. (M) Functional enrichment modules of 82 DEGs, colored by module type.
Subsequently, genes expression was compared among various time points. The difference analysis revealed 125 downregulated genes in ARD614 compared with pre-treatment, and 99 upregulated genes in RespD376 compared with ARD614 (Fig 3K). A total of 82 genes were strongly associated with disease progression after intersecting these genes (Fig 3L). The 82 genes mainly focused on cellular response to interferon−gamma, regulation of lymphocyte proliferation, and MHC protein complex binding (Fig 3M). Interferons (IFNs) consisted of different types such as class I IFN (α, β, ε, κ, ν, ω, τ, δ, ζ), class II IFN (γ), and class III IFN (λ). IFNs not only initiated anti-viral but also anti-tumor responses. In addition, class I and II IFNs caused an up-regulation of MHC molecules on the surface of the tumor cells [14,27].
B cell and plasmablast changes across the MCC disease trajectory
We first extracted 3434 cells identified as B-cell lineages from integrated scRNA-seq data. UMAP package identified three distinct large clusters, including 778 naive B cells (NB), 210 transitional B cells (Trans), and 1142 plasmablasts (PB, Fig 4A). Since there were only 58 B cells at the stage of EarlyD27 and the cells at this stage were not of great significance for us to investigate the relapse mechanism, these 58 cells were excluded in the subsequent analysis (Fig 4B). We observed a significantly decreased abundance of NB (up-regulation of MS4A1 and IGHD) and increased abundance of PB (up-regulation of CD38) in ARD614 compared with Pre and RespD376 (Fig 4C).
(A) The UMAP projection of 3434 B cells, colored according to B cell subtype (right) and time point (left). (B) The number of B cells at each time point. (C) The proportion changes of three B cell subtypes at each time point. (D-F) Dot plots for disease trajectory signature genes in naive B cells (NB), transitional B cells (Trans), and plasmablasts (PB). Genes were selected on the basis of the expression amount of the most characteristic genes. Color discriminates genes with upregulation (red) or downregulation (blue) expression, and point size represents the number of cells per group expressing the corresponding gene. (G) Metabolic pathways enriched in B cell subtypes. Top 20 active metabolic pathways in PBs are shown. Significant differences in metabolic activity were determined by using a Kruskal-Wallis test. (H) Proportion change of plasma cells. (I) Expression of HLA class II genes in the validation cohort.
We next explored the longitudinal gene expression patterns of NBs, PBs and Trans (Fig 4D–4F). The expression level of MHC class II genes (HLA-DQA1, HLA-DMB) in NB was associated with disease progression. The expression of MHC II genes was highest at the time point of Pre and significantly down-regulated at the ARD614 (Fig 4D), which suggested a disorder of immune crosstalk between adaptive immune cells and a decline in B cell antigen presentation. The chemokine receptors (CCL5, CCL4, and CCR7) were present in the PBs at the stage of RespD376, yet such genes were absent from PBs in patients with recurrence (ARD614) (Fig 4E). This could inhibit germinal center (GC) reactions and ultimately lead to the dysregulated persistent immune responses in the later stage of immunotherapy [28]. The expression of the activator protein-1 complex member (JUN) was upregulated in the Trans at RespD376, which may indicate that the tissue was primed for malignant proliferation before the appearance of recurrence features (Fig 4F). The expression of G0S2 which was a potential regulator of metabolic changes increased at the stage of ARD614 [29]. Thus, metabolic pathway activities were assessed using a dataset of 85 metabolic pathways described in previous studies [23] (Fig 4G). There was a significant increase in glycolysis and gluconeogenesis in RespD376. Although Resp376 is a phase of recovery for MCC patients, it may indicate a delayed response between activation of drug metabolism and acquisition of drug resistance. PBs from the RespD376 state displayed a high metabolic activity, which was reduced only upon the ARD614 stage. Many metabolic processes in PBs increased: including propanoate and pyruvate metabolism, arachidonic acid metabolism, steroid biosynthesis, glycolysis and gluconeogenesis. PBs can act as a nutrient sink that regulates the immune response [30]. Metabolic hyperactivity of PBs in ARD614 may result in nutrient deprivation of the GC reaction which will hinder the generation of long-lived plasma cell and memory B cell responses [30]. In bulk RNA data, the abundance of PBs was found to be higher in MCC samples compared to normal samples (Fig 4H). Similarly, the NB marker HLA-DQA1 was upregulated in normal tissues compared to MCC (Fig 4I). Altogether, the analysis identified the widespread activation of PBs and suggested that increased drug metabolism and inhibition of GC reaction by the rapid development of PBs may be one of the main mechanisms of drug resistance and relapse.
S100AhighHLA-DR/DPlow monocyte subpopulation reflects coordinated changes with both the immunotherapy progression and MCC severity
UMAP visualization of all 2489 monocytes showed significant separation of non-classical and classical monocytes (Fig 5A). Classical monocytes were characterized by the upregulated transcripts CD14, FCN1, VCAN, and GPX1, while non-classical monocytes were identified based on the up-regulation of FCGR3A, CDKN1C, POU2F2, and MS4A7 (Fig 5B). The classical monocyte cluster exhibited expression of inflammation-related transcripts such as S100A8, S100A9, S100A10, and S100A12, and decreased levels of MHC II transcripts such as HLA-DR/DP (Fig 5B and 5C). The classical monocyte (S100AhighHLA-DR/DPlow) fraction increased in Pre and ARD614, compared to the RespD376. The non-classical monocyte cluster (S100AlowHLA-DR/DPhigh) percentage increased for RespD376 stage (Fig 5D). Interestingly, the characteristic gene expression patterns of the two cell subtypes in the protein family (Type I IFN response genes IFITM3 and IFI27 in the classical cluster, IFI27 and IFI30 in the non-classical cluster, for instance) were both highly expressed in RespD376, yet these genes were absent from monocytes in ARD614. Meanwhile, general down-regulation of the pro-inflammatory factors such as TNF family (TNF, TNFAIP2, TNFAIP3 and TNFSF10) and chemokine family (CCL3, CCL4 and CX3CR1) were present in RespD376 until disease relapse (ARD614) (Fig 5E and 5F). Differential gene analysis of these two cell types revealed that the classical monocyte had 35 genes upregulated and 35 genes downregulated, while the non-classical monocyte had 51 genes upregulated and 40 genes downregulated in expression (Fig 5G). The proportion changes and transcript characteristics of the two monocyte subtypes at different periods seemed to suggest that the improvement of MCC benefited from a peripheral environment with low inflammation and high immune response, prompting us to explore this conjecture more.
(A) UMAP embedding of all monocytes colored by monocyte compartment subtype. (B) Violin plots showing the expression of non-classical and classical monocytes, S100A, HLA-DR/DP markers in monocyte subtypes. (C) Several cell phenotypic markers were used to identify monocyte cell subtypes. (D) The proportion changes of two monocyte subtypes at each time point. (E) Dot plot for disease trajectory signature genes in classical monocytes. Color discriminates genes with upregulation (red) or downregulation (blue) expression, and point size represents the number of cells per group expressing the corresponding gene. (F) Dot plot for disease trajectory signature genes in non-classical monocytes. (G) Number of DEGs between RespD376 and ARD614 in two monocyte subtypes. (H) The top 10 GSVA results of non-classical monocytes. Y-axis represents the expression level of the entire pathway. (I) Changes in levels of “antigen processing and presentation pathways” in two monocyte subtypes. (J) Correlation between pathway “primary immunodeficiency” and pathway “antigen processing and presentation” in monocytes.
The GSVA results of non-classical monocytes (S100AlowHLA-DR/DPhigh) showed that apoptosis, B cell receptor signaling pathway, glycolysis gluconeogenesis, natural killer cell mediated cytotoxicity and T cell receptor signaling pathway were more active during convalescence (pseudotime RespD376) (Fig 5H). The corresponding biological significance of HLA class II genes expression signatures were shown in Fig 5I. In both classical and non-classical monocytes, the ability of antigen processing and presentation in convalescence was stronger than those in the pre-treatment and relapse. Furthermore, antigen processing and presentation was negatively correlated with the primary immunodeficiency (Fig 5J), suggesting that normal functioning of monocyte after treatment may significantly reduce tumor pressure until recurrence.
Discussion
A comprehensive understanding of the various cellular changes in MCC patients is essential to determine the effectiveness of treatment, predict disease prognosis, and understand the heterogeneity of reported disease severity. In this report, we mainly discussed the following aspects:
The significant difference observed between the response phase and the relapse phase in peripheral immune cells is that exhausted CD8+ T cells (T cells that gradually lose their function) [31] were significantly elevated and exhibited abnormal phenotypes. Previous studies found that increased infiltration of CD8+ T cells effectively improve the prognosis of patients with MCC [32]. The increase of CD8+ T cells was also observed in the response period (RespD376), which seemed to explain the decrease in tumor pressure in MCC patients during the RespD376. We also identified a cluster of CD8+ T cells that amplified during the RespD376 and showed both exhaustion and proliferation phenotypes. This was consistent with the previous finding that infiltrating T cells were usually characterized by an exhausted phenotype [33,34]. The impact of exhausted phenotype T cell effectors on patients may be greater than the activity of their CD8+ T cells expanded in vitro, leading to the gradual depletion of these T cells used for immunotherapy and loss of their therapeutic efficacy. Meanwhile, the expression of immune checkpoints LAG3, CTLA-4, PD-1, and TIGIT in this cluster was higher than those in other clusters. Both exhaustion of cells and upregulation of these immune checkpoints inhibit the immune response and can be used by tumor cells as an immune escape mechanism, which further led to the failure of immunotherapy [35]. In addition, the results revealed that T cell infiltration level was not an important independent predictor of survival in MCC patients.
The results of our longitudinal study identified two other cellular features associated with MCC recurrence that previous studies had not recognized. First, the relative abundance of PBs with low expression of chemokine receptors (CCL5, CCL4, and CCR7) and high metabolic activity increased during the relapse phase (ARD614). Dynamic changes in chemokine receptor expression induce antigen-stimulated B cells to target specific lymphoid tissue sub-compartments, resulting in further activation and proliferation of B cells [36]. During the immune process, B cells upregulate the receptor for CCL19/21 (CCR7), produced in the paracortex, causing them to migrate to the T cell zone. Activated B cells also secrete T chemokines CCL3 and CCL4, thus increasing the efficiency of interaction with activated T cells that also migrate to the T:B boundary [37]. So, high expression of chemokines not only activate B cells to produce an effective immune response but also improve the synergistic efficiency of T cells and B cells. At the same time, both the dysregulation of these chemokines and the metabolic hyperactivity of PBs in the ARD614 will hinder the formation and sustainment of germinal centers, which will lead to the interruption of B cell-mediated long-lasting and broadly protective antiviral immune response, further compromising the effectiveness of immunotherapy [30,38].
A significant increase of high inflammatory classical monocytes was identified, which carried a weak MHC II signature. Chronic inflammatory diseases such as rheumatoid arthritis was associated with a higher incidence of MCC [39]. There was also a significant correlation between inflammatory infiltration and PD-L1 expression in tumor cells [40]. Hence, we speculated that the upregulation of pro-inflammatory factors and inflammatory transcripts at the ARD614 may create a local pro-inflammatory environment that drives the expression of tumor PD-L1 [41]. Meanwhile, extensive downregulation of MHC II was observed in T cells, B cells, and monocytes during the ARD614, which confirmed the importance of immunity as a prognostic determinant of MCC patients. These results clearly suggested that regulatory events in T cells, B cells, and monocytes might act as pivotal components of an unfavorable course of MCC, which mandates further prospective exploration.
There are several limitations of our study. First, the work contains no wet lab cross-check. Second, with the current scRNA-seq strategy, peripheral blood samples are not a complete representation of the true state of MCC. Therefore, future studies with new scRNA-seq data and validation in a larger clinical cohort may help to further investigate the mechanisms of immunotherapy tolerance in MCC.
Conclusions
In this study, computational methods were used to integrate single-cell characteristics at different time points to form an integrated and comprehensive view of MCC that relapsed after immunotherapy. A series of adverse factors that may lead to immune escape in MCC were determined, including T cell exhaustion and immune crosstalk disorders. The peripheral environment with high inflammation and low immune response may also lead to resistance to immunotherapy.
Acknowledgments
We sincerely thank our respective college and institutes for technical assistance and valuable support in the completion of this project. We are also immensely grateful to Dr.Jian-ming Zeng (University of Macau), and all the members of his bioinformatics team, for generously sharing their experience and codes.
References
- 1. Paulson KG, Park SY, Vandeven NA, Lachance K, Thomas H, Chapuis AG, et al. Merkel cell carcinoma: Current US incidence and projected increases based on changing demographics. Journal of the American Academy of Dermatology. 2018;78(3):457–63.e2. Epub 2017/11/06. pmid:29102486; PubMed Central PMCID: PMC5815902.
- 2. Becker JC, Stang A, DeCaprio JA, Cerroni L, Lebbé C, Veness M, et al. Merkel cell carcinoma. Nature reviews Disease primers. 2017;3:17077. Epub 2017/10/27. pmid:29072302; PubMed Central PMCID: PMC6054450.
- 3. Feng H, Shuda M, Chang Y, Moore PS. Clonal integration of a polyomavirus in human Merkel cell carcinoma. Science (New York, NY). 2008;319(5866):1096–100. Epub 2008/01/19. pmid:18202256; PubMed Central PMCID: PMC2740911.
- 4. Lemos BD, Storer BE, Iyer JG, Phillips JL, Bichakjian CK, Fang LC, et al. Pathologic nodal evaluation improves prognostic accuracy in Merkel cell carcinoma: analysis of 5823 cases as the basis of the first consensus staging system. Journal of the American Academy of Dermatology. 2010;63(5):751–61. Epub 2010/07/22. pmid:20646783; PubMed Central PMCID: PMC2956767.
- 5. Bichakjian CK, Olencki T, Aasi SZ, Alam M, Andersen JS, Blitzblau R, et al. Merkel Cell Carcinoma, Version 1.2018, NCCN Clinical Practice Guidelines in Oncology. Journal of the National Comprehensive Cancer Network: JNCCN. 2018;16(6):742–74. Epub 2018/06/13. pmid:29891526.
- 6. Nghiem PT, Bhatia S, Lipson EJ, Kudchadkar RR, Miller NJ, Annamalai L, et al. PD-1 Blockade with Pembrolizumab in Advanced Merkel-Cell Carcinoma. The New England journal of medicine. 2016;374(26):2542–52. Epub 2016/04/20. pmid:27093365; PubMed Central PMCID: PMC4927341.
- 7. Nghiem P, Kaufman HL, Bharmal M, Mahnke L, Phatak H, Becker JC. Systematic literature review of efficacy, safety and tolerability outcomes of chemotherapy regimens in patients with metastatic Merkel cell carcinoma. Future oncology (London, England). 2017;13(14):1263–79. Epub 2017/03/30. pmid:28350180; PubMed Central PMCID: PMC6040046.
- 8. Iyer JG, Blom A, Doumani R, Lewis C, Tarabadkar ES, Anderson A, et al. Response rates and durability of chemotherapy among 62 patients with metastatic Merkel cell carcinoma. Cancer medicine. 2016;5(9):2294–301. Epub 2016/07/20. pmid:27431483; PubMed Central PMCID: PMC5055152.
- 9. Sharma P, Allison JP. The future of immune checkpoint therapy. Science (New York, NY). 2015;348(6230):56–61. Epub 2015/04/04. pmid:25838373.
- 10. Knepper TC, Montesion M, Russell JS, Sokol ES, Frampton GM, Miller VA, et al. The Genomic Landscape of Merkel Cell Carcinoma and Clinicogenomic Biomarkers of Response to Immune Checkpoint Inhibitor Therapy. Clinical cancer research: an official journal of the American Association for Cancer Research. 2019;25(19):5961–71. Epub 2019/08/11. pmid:31399473; PubMed Central PMCID: PMC6774882.
- 11. Gide TN, Wilmott JS, Scolyer RA, Long GV. Primary and Acquired Resistance to Immune Checkpoint Inhibitors in Metastatic Melanoma. Clinical cancer research: an official journal of the American Association for Cancer Research. 2018;24(6):1260–70. Epub 2017/11/12. pmid:29127120.
- 12. Syn NL, Teng MWL, Mok TSK, Soo RA. De-novo and acquired resistance to immune checkpoint targeting. The Lancet Oncology. 2017;18(12):e731–e41. Epub 2017/12/07. pmid:29208439.
- 13. Kelderman S, Schumacher TN, Haanen JB. Acquired and intrinsic resistance in cancer immunotherapy. Molecular oncology. 2014;8(6):1132–9. Epub 2014/08/12. pmid:25106088; PubMed Central PMCID: PMC5528612.
- 14. Paulson KG, Tegeder A, Willmes C, Iyer JG, Afanasiev OK, Schrama D, et al. Downregulation of MHC-I expression is prevalent but reversible in Merkel cell carcinoma. Cancer immunology research. 2014;2(11):1071–9. Epub 2014/08/15. pmid:25116754; PubMed Central PMCID: PMC4221542.
- 15. Ritter C, Fan K, Paulson KG, Nghiem P, Schrama D, Becker JC. Reversal of epigenetic silencing of MHC class I chain-related protein A and B improves immune recognition of Merkel cell carcinoma. Scientific reports. 2016;6:21678. Epub 2016/02/24. pmid:26902929; PubMed Central PMCID: PMC4763224.
- 16. Willmes C, Adam C, Alb M, Völkert L, Houben R, Becker JC, et al. Type I and II IFNs inhibit Merkel cell carcinoma via modulation of the Merkel cell polyomavirus T antigens. Cancer research. 2012;72(8):2120–8. Epub 2012/03/06. pmid:22389452.
- 17. Le DT, Durham JN, Smith KN, Wang H, Bartlett BR, Aulakh LK, et al. Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Science (New York, NY). 2017;357(6349):409–13. Epub 2017/06/10. pmid:28596308; PubMed Central PMCID: PMC5576142.
- 18. Grün D, van Oudenaarden A. Design and Analysis of Single-Cell Sequencing Experiments. Cell. 2015;163(4):799–810. Epub 2015/11/07. pmid:26544934.
- 19. Paulson KG, Voillet V, McAfee MS, Hunter DS, Wagener FD, Perdicchio M, et al. Acquired cancer resistance to combination immunotherapy from transcriptional loss of class I HLA. Nature communications. 2018;9(1):3868. Epub 2018/09/27. pmid:30250229; PubMed Central PMCID: PMC6155241.
- 20. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nature biotechnology. 2018;36(5):411–20. Epub 2018/04/03. pmid:29608179; PubMed Central PMCID: PMC6700744.
- 21. Wang Y, Wang R, Zhang S, Song S, Jiang C, Han G, et al. iTALK: an R Package to Characterize and Illustrate Intercellular Communication. bioRxiv. 2019:507871.
- 22. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC bioinformatics. 2013;14:7. Epub 2013/01/18. pmid:23323831; PubMed Central PMCID: PMC3618321.
- 23. Xiao Z, Dai Z, Locasale JW. Metabolic landscape of the tumor microenvironment at single cell resolution. Nature communications. 2019;10(1):3763. Epub 2019/08/23. pmid:31434891; PubMed Central PMCID: PMC6704063.
- 24. Bengsch B, Ohtani T, Khan O, Setty M, Manne S, O’Brien S, et al. Epigenomic-Guided Mass Cytometry Profiling Reveals Disease-Specific Features of Exhausted CD8 T Cells. Immunity. 2018;48(5):1029–45.e5. Epub 2018/05/17. pmid:29768164; PubMed Central PMCID: PMC6010198.
- 25. Guo X, Zhang Y, Zheng L, Zheng C, Song J, Zhang Q, et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nature medicine. 2018;24(7):978–85. Epub 2018/06/27. pmid:29942094.
- 26. Qian J, Olbrecht S, Boeckx B, Vos H, Laoui D, Etlioglu E, et al. A pan-cancer blueprint of the heterogeneous tumor microenvironment revealed by single-cell profiling. Cell research. 2020;30(9):745–62. Epub 2020/06/21. pmid:32561858; PubMed Central PMCID: PMC7608385.
- 27. Gehrcken L, Sauerer T, Schaft N, Dörrie J. T-Cell Responses in Merkel Cell Carcinoma: Implications for Improved Immune Checkpoint Blockade and Other Therapeutic Options. International journal of molecular sciences. 2021;22(16). Epub 2021/08/28. pmid:34445385; PubMed Central PMCID: PMC8395396.
- 28. Reimer D, Lee AY, Bannan J, Fromm P, Kara EE, Comerford I, et al. Early CCR6 expression on B cells modulates germinal centre kinetics and efficient antibody responses. Immunology and cell biology. 2017;95(1):33–41. Epub 2016/07/29. pmid:27465674.
- 29. Heckmann BL, Zhang X, Xie X, Liu J. The G0/G1 switch gene 2 (G0S2): regulating metabolism and beyond. Biochimica et biophysica acta. 2013;1831(2):276–81. Epub 2012/10/04. pmid:23032787; PubMed Central PMCID: PMC3698047.
- 30. Vijay R, Guthmiller JJ, Sturtz AJ, Surette FA, Rogers KJ, Sompallae RR, et al. Infection-induced plasmablasts are a nutrient sink that impairs humoral immunity to malaria. Nature immunology. 2020;21(7):790–801. Epub 2020/05/20. pmid:32424361; PubMed Central PMCID: PMC7316608.
- 31. Ritter C, Fan K, Paschen A, Reker Hardrup S, Ferrone S, Nghiem P, et al. Epigenetic priming restores the HLA class-I antigen processing machinery expression in Merkel cell carcinoma. Scientific reports. 2017;7(1):2290. Epub 2017/05/26. pmid:28536458; PubMed Central PMCID: PMC5442125 interest. Author S.U. has received advisory board honorariums from Roche, BMS and Novartis. Author J.C.B. has received speaker honorariums from Amgen, MerckSerono, and Pfizer, advisory board honorariums from Amgen, CureVac, eTheRNA, Lytix, MerckSerono, Novartis, Rigontec, and Takeda as well as research funding from Boehringer Ingelheim, BMS and MerckSerono; none of these activities are related to the submitted report.
- 32. Paulson KG, Iyer JG, Tegeder AR, Thibodeau R, Schelter J, Koba S, et al. Transcriptome-wide studies of merkel cell carcinoma and validation of intratumoral CD8+ lymphocyte invasion as an independent predictor of survival. Journal of clinical oncology: official journal of the American Society of Clinical Oncology. 2011;29(12):1539–46. Epub 2011/03/23. pmid:21422430; PubMed Central PMCID: PMC3082974 found at the end of this article.
- 33. Walsh NM, Fleming KE, Hanly JG, Dakin Hache K, Doucette S, Ferrara G, et al. A morphological and immunophenotypic map of the immune response in Merkel cell carcinoma. Human pathology. 2016;52:190–6. Epub 2016/03/17. pmid:26980039.
- 34. Dowlatshahi M, Huang V, Gehad AE, Jiang Y, Calarese A, Teague JE, et al. Tumor-specific T cells in human Merkel cell carcinomas: a possible role for Tregs and T-cell exhaustion in reducing T-cell responses. The Journal of investigative dermatology. 2013;133(7):1879–89. Epub 2013/02/20. pmid:23419694; PubMed Central PMCID: PMC3691077.
- 35. Torphy RJ, Schulick RD, Zhu Y. Newly Emerging Immune Checkpoints: Promises for Future Cancer Therapy. International journal of molecular sciences. 2017;18(12). Epub 2017/12/07. pmid:29211042; PubMed Central PMCID: PMC5751245.
- 36. Lam JH, Smith FL, Baumgarth N. B Cell Activation and Response Regulation During Viral Infections. Viral immunology. 2020;33(4):294–306. Epub 2020/04/25. pmid:32326852; PubMed Central PMCID: PMC7247032.
- 37. Damdinsuren B, Zhang Y, Khalil A, Wood WH 3rd, Becker KG, Shlomchik MJ, et al. Single round of antigen receptor signaling programs naive B cells to receive T cell help. Immunity. 2010;32(3):355–66. Epub 2010/03/17. pmid:20226693; PubMed Central PMCID: PMC3607434.
- 38. Ryg-Cornejo V, Ioannidis LJ, Ly A, Chiu CY, Tellier J, Hill DL, et al. Severe Malaria Infections Impair Germinal Center Responses by Inhibiting T Follicular Helper Cell Differentiation. Cell reports. 2016;14(1):68–81. Epub 2016/01/05. pmid:26725120.
- 39. Sahi H, Sihto H, Artama M, Koljonen V, Böhling T, Pukkala E. History of chronic inflammatory disorders increases the risk of Merkel cell carcinoma, but does not correlate with Merkel cell polyomavirus infection. British journal of cancer. 2017;116(2):260–4. Epub 2016/12/16. pmid:27978533; PubMed Central PMCID: PMC5243985.
- 40. Lyford-Pike S, Peng S, Young GD, Taube JM, Westra WH, Akpeng B, et al. Evidence for a role of the PD-1:PD-L1 pathway in immune resistance of HPV-associated head and neck squamous cell carcinoma. Cancer research. 2013;73(6):1733–41. Epub 2013/01/05. pmid:23288508; PubMed Central PMCID: PMC3602406.
- 41. Lipson EJ, Vincent JG, Loyo M, Kagohara LT, Luber BS, Wang H, et al. PD-L1 expression in the Merkel cell carcinoma microenvironment: association with inflammation, Merkel cell polyomavirus and overall survival. Cancer immunology research. 2013;1(1):54–63. Epub 2014/01/15. pmid:24416729; PubMed Central PMCID: PMC3885978.