Skip to main content

ORIGINAL RESEARCH article

Front. Immunol., 16 May 2022
Sec. Cancer Immunity and Immunotherapy

Development of a Phagocytosis-Dependent Gene Signature to Predict Prognosis and Response to Checkpoint Inhibition in Clear-Cell Renal Cell Carcinoma

Updated
Kunping Li&#x;Kunping Li1†Yuqing Li&#x;Yuqing Li1†Yinfeng Lyu&#x;Yinfeng Lyu1†Linyi TanLinyi Tan1Xinyi ZhengXinyi Zheng2Haowen Jiang*Haowen Jiang1*Hui Wen*Hui Wen1*Chenchen Feng,,*Chenchen Feng1,2,3*
  • 1Department of Urology, Huashan Hospital, Fudan University, Shanghai, China
  • 2Department of Pharmacology, Huashan Hospital, Fudan University, Shanghai, China
  • 3Beijing Advanced Innovation Center for Food Nutrition and Human Health, Beijing Technology and Business University (BTBU), Beijing, China

Aim: The action of immune checkpoint inhibition (ICI) largely depends on antibody-dependent cellular phagocytosis (ADCP). We thus aim to develop ADCP-based ccRCC risk stratification as both prognostic and therapeutic markers of ICI.

Method: Genomic data from multiple public datasets (TCGA, etc.) were integrated. A cancer-intrinsic ADCP gene set for ccRCC tailored from a recent report was constructed based on the association with prognosis, immune infiltrates, and response to ICI. Therapeutic potential was profiled using genome-drug sensitivity datasets.

Results: ADCP genes were selected from a recent CRISPR/Cas9 screen report. Following a four-module panel based on clinical traits, we generated a six-gene signature (ARPC3, PHF19, FKBP11, MS4A14, KDELR3, and CD1C), which showed a strong correlation with advanced grade and stage and worsened prognosis, with a nomogram showing predictive efficacies of 0.911, 0.845, and 0.867 (AUC) at 1, 3, and 5 years, respectively. Signatures were further dichotomized, and groups with a higher risk score showed a positive correlation with tumor mutation burden, higher expressions of inhibitory checkpoint molecules, and increased antitumor immune infiltrates and were enriched for antitumor immune pathways. The high risk-score group showed better response to ICI and could benefit from TKIs of axitinib, tivozanib, or sorafenib, preferentially in combination, whereas sunitinib and pazopanib would better fit the low risk-score group.

Conclusion: Here we showed a six-gene ADCP signature that correlated with prognosis and immune modulation in ccRCC. The signature-based risk stratification was associated with response to both ICI and tyrosine kinase inhibition in ccRCC.

Highlights

* Antibody-dependent cellular phagocytosis (ADCP) is critical for action of immunotherapy and has yet been reported in kidney cancer.

* We for the first time developed an ADCP-based gene signature encompassing ARPC3, PHF19, FKBP11, MS4A14, KDELR3, and CD1C.

* The signature showed a strong correlation with prognosis and response to immunotherapy and may predict response to tyrosine kinase inhibitors.

* The signature holds promise as both a therapeutic and prognostic marker in kidney cancer.

Introduction

Clear-cell renal cell carcinoma (ccRCC) is the most common subtype of renal cancer and is curable at the localized stage. Metastatic ccRCC consists of ~30% of cases and faces challenges in definitive treatment. Given the immunogenic nature of ccRCC, immunotherapy has re-immerged after initial attempts using interleukins and tumor-necrosis factor (TNF) (1, 2). The current frontline treatment for metastatic disease mandates the inclusion of immune checkpoint inhibition (ICI) that targets inhibitory immune checkpoints including PD-1, PD-L1, and CTLA-4 (3). Recent trials have demonstrated overwhelming advantages of the combination of ICI and tyrosine kinase inhibitors (TKIs), leaving monotherapy of TKI solely indicated for patients with favorable International Metastatic RCC Database Consortium (IMDC) risk (4).

However, the IMDC risk score is based on clinicopathological parameters that are dismally related to the pharmaceutical mechanisms of either TKI or ICI, let alone the synergy of combination (5). Besides the major anti-angiogenic effect and minor cancer-intrinsic inhibition effect of TKIs, the synergy is believed to result from increased immune infiltrates due to local inflammation by TKI and subsequent ICI antibodies kick in to facilitate the recognition of tumor cells by immune cells (6). Nonetheless, no consensus has been made on predictive marker(s) for ICI response even in a single type of cancer despite various attempts. PD-1/PD-L1 positivity, tumor mutation burden, and microsatellite instability (MSI) together with multiple experimental markers all showed efficacy in a certain context and can hardly be extrapolated across cancer types (7). Thus far, no genetic marker has been approved in the clinical setting for ccRCC.

Antibody-dependent cellular phagocytosis (ADCP) is the terminal step of ICI in which macrophages eliminate tumor cells upon previous activation of recognition (8). Albeit abundant in number, negative regulatory factors that dampen phagocytic activity remain the key obstacle hampering the full efficacy of ICI (9, 10). In ccRCC, the objective response rate of lenvatinib plus pembrolizumab reached 71%, among which 16.1% had complete response (11). While it is plausible to see a rapid change in the treatment paradigm in ccRCC, there are still ~30% of cases that fail to benefit from the combination. Whether ADCP plays a role in identifying potential beneficiaries from ICI and TKI in ccRCC remains unreported, and we thus aim to evaluate both prognostic and predictive merits with an in silico approach.

Materials and Methods

Data Acquirement and Processing

Processed and standardized RNA-seq data and clinical data of The Cancer Genome Atlas of Kidney Renal Clear Cell Carcinoma (TCGA-KIRC) cohort were downloaded from UCSC Xena (https://xena.ucsc.edu/). The E-MTAB-1980 dataset was acquired from the ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/). Patients without complete prognostic data were eliminated. Then, patients with overall survival (OS) for more than 30 days were included in the subsequent analysis. The gene hits, identified with the CRISPR/Cas9 screen according to the 5% false discovery rate (FDR) or 95% credible interval, were collected from the previous research (8). These genes were defined as antibody-dependent cellular phagocytosis (ADCP)-related genes. Details on study designs have been published in a previous research (8).

Co-Expression Module Construction

Weighted gene co-expression network analysis (WGCNA) is a system biology method used to describe the correlation patterns between different genes based on expression data. ADCP-related genes in TCGA-KIRC cohort were used to perform WGCNA analysis by the WGCNA R software package (12). GoodSamplesGenes function was used to eliminate outlier samples and genes. The appropriate soft power value was determined according to scale independence (more than 0.8). According to topological overlap matrix (TOM)-based dissimilarities, ADCP-related genes with similar expression profiles were classified into the same gene modules. The minimum number of genes was set as 30. The correlation between the module eigengene and the phenotype was evaluated by the Spearman correlation test. To functionally annotate different gene modules, the Metascape database (http://metascape.org) was utilized to annotate and visualize Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis.

Gene Signature Construction

The gene module with the highest correlation with clinical phenotypes was enrolled in a subsequent construction of gene signature. Univariate COX regression was performed on the training set, consisting of 97 ADCP-related genes and 513 ccRCC patients from TCGA-KIRC cohort, to select prognostic genes, and a cutoff of p-value was 0.01. Least absolute shrinkage and selection operator (LASSO) regression analysis is a common shrinkage method which can screen appropriate variables from multicollinear and high-dimensional data. In our study, LASSO regression was implemented on prognostic genes from the previous step using the “glmnet” R software package to screen the most valuable prognostic candidates. Next, stepwise regression was used to further shrink variables and select the best model based on the minimum Akaike information criterion (AIC) value principle by the “My.stepwise” R software package. Finally, the gene signature was constructed by multivariate COX regression. Based on the gene signature, the risk score of each patient was calculated by the formula below:

Riskscore =iCoefficient (i)Expression of gene(i)

Coefficient (i) is the regression coefficient of gene (i) in the multivariate COX regression model; Expression of gene (i) is the expression value of (i). Patients were dichotomized into high- and low-risk groups based on risk scores. The survival analysis of different risk groups and genes in signature was performed by “survival” R software packages. Additionally, the “survminer” R software package was used to calculate the optimal cutoff value of gene expression. The “survivalROC” R software package was applied to plot receiver operating characteristic (ROC) curves to estimate the predictive accuracy of the gene signature. The E-MTAB-1980 cohort was utilized to validate the above results.

Nomogram Construction

A nomogram was constructed based on risk score and clinical phenotypes by multivariate COX regression through “rms” and “survminer” R software packages. All-subset regression was performed to find the optimal model by the “leaps” R software package. The consistency between predicted and actual survival outcomes was evaluated according to calibration curves. Time-dependent ROC curves were plotted to estimate the predictive accuracy of the nomogram.

Somatic Variant Analysis

Genetic mutation information was downloaded from TCGA database. The “maftools” package was used to calculate tumor mutation burden (TMB), visualize, and compare the frequency of genetic mutation between high- and low-risk groups. Copy number alteration (CNA) data of TCGA-KIRC were downloaded from the cBioPortal database (https://www.cbioportal.org/). The CNA burden was defined as the total number of genes with copy number gain or loss (13, 14).

Identification of Differentially Expressed Genes and Functional Enrichment Analysis

The differentially expressed genes (DEGs) between high- and low-risk groups were identified by the limma R package (15). GO and KEGG enrichment analyses of upregulated genes in the high-risk group were implemented and visualized by the “clusterProfiler” R package (16). Gene set enrichment analysis (GSEA) was performed to analyze significantly enriched pathways in the high-risk group by “GSEA” software. In GSEA, significantly enriched pathways met the following criteria: normalized enrichment score >1; nominal p-value < 0.05; and FDR q-value < 0.25. The gene set of “c5.go.bp.v7.4.symbols.gmt” was selected to perform GSEA.

Immune Cell Infiltrate Analysis

ESTIMATE, CIBERSORT, and single-sample gene set enrichment analysis (ssGSEA) algorithms were applied to calculate the abundance of immune cell infiltrations between high- and low-risk groups (1719).

Correlation of the Gene Signature With Immunotherapy

Immunophenoscore (IPS) has been verified as a signature to predict the response to immunotherapy (20). The IPS score was calculated based on the expression of the representative genes or gene sets: immunomodulators, MHC molecules, effector cells, and suppressor cells. The IPS score of TCGA-KIRC was downloaded from The Cancer Immunome Atlas (TCIA) database (https://tcia.at/home). RNA-seq profiles and clinical data of anti-PD-1 antibody nivolumab in ccRCC were derived from three clinical trials: CheckMate 009 (CM-009; NCT01358721), CheckMate 010 (CM-010; NCT01354431), and CheckMate 025 (CM-025; NCT01668784) (2123). The normalized and integrated data were downloaded from the supplementary data of the previous research (24). Risk scores for the above patients were calculated based on the gene signature.

Prediction of Chemotherapeutic Effect and Exploring Potential Chemotherapy Drugs in ccRCC

For personalized treatment, the sensitivity of chemotherapy was predicted based on the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) database. The “pRRophetic” R package (25) was applied to estimate the half-maximal inhibitory concentration (IC50) of some common chemotherapy drugs in ccRCC. The DEGs between two risk groups were uploaded into the connectivity map (cMAP, https://portals.broadinstitute.org/cmap/) database (26) to explore potential chemotherapy drugs. The enrichment score varied from -1 to 1. The cutoff of the p-value was 0.05.

Sample Collection and Immunohistochemistry

A standard immunohistochemistry (IHC) protocol was followed. In brief, sections were sliced consecutively at 4 μm and were then deparaffinized followed by gradient rehydration in ethanol and subsequently rinsed. Samples were then mounted on polylysine-coated glass. A pilot staining was conducted using multiple antibodies. Multiple primary antibodies were tested, and antibodies against FKBP11 (Sigma-Aldrich, St. Louis, MO, USA, Cat. HPA041709; dilution 1:500), ARPC3 (Sigma-Aldrich, Cat. HPA006550; dilution 1:500), KDELR3 (Sigma-Aldrich, Cat. HPA043477; dilution 1:500), and MS4A14 (Abcam, Cat. ab182151) were adopted for the final evaluation. Assessment of positivity was referenced from The Human Protein Atlas (27), as none of the factors were reported before in ccRCC using IHC staining. According to the intensity and extensity of staining, we designated a three-tiered scoring system for FKBP11 and ARPC3 and a dichotomized scoring system for KDELR3 and MS4A14. However, CD1C and PHF19 were not IHC-detectable in ccRCC despite multiple trials with different antibodies, which was also supported in The Human Protein Atlas (27). The sections treated with PBS in lieu of the primary antibody were chosen as negative controls, and tumor-infiltrating lymphocytes (TILs) were stained and evaluated using the H-score, as previously reported (28). All patients signed an informed consent, and the study was approved by the Huashan Institutional Review Board (HIRB-2022-204).

Statistical Analysis

The Kaplan–Meier log-rank test was applied to perform survival analysis. The Wilcoxon test was used to evaluate continuous-variable data. The chi-square test and Fisher’s exact test were performed to analyze differences between categorical variables. The Spearman correlation analysis was used to calculate the correlation coefficient. The IHC data were regarded as non-parametric, and comparisons for medians were analyzed using Mann–Whitney’s U-test or Kruskal–Wallis test with post-hoc Dunn’s test. The p value of <.05 was accepted as statistically significant. The visualization of results was performed by “ggplot2,” “ggpubr,” and “ggstatsplot” R software packages. Statistical analysis and visualization were implemented by R software (version 4.1.1) and GraphPad software.

Results

Data Acquirement

The training set consisted of 513 ccRCC patients and 72 normal samples from TCGA-KIRC cohort. The testing set consisted of 100 ccRCC patients from the E-MTAB-1980 cohort. The detailed clinical information of training and testing cohorts is shown in Table 1. A total of 543 ADCP-related genes were extracted from a previous research (8), of which 511 were found in TCGA cohort (Figure 1A). In the CRISPR activation (CRISPRa) screen, genes with the combo casTLE effect less than 0 were considered as anti-phagocytic genes, whereas those with the combo casTLE effect more than 0 were pro-phagocytic genes. In the CRISPR knockout (CRISPRko) screen, the definition was the opposite. casTLE, which stands for cas9 High Throughput maximum Likelihood Estimator, was utilized to evaluate the gene effect sizes comparing negative controls in the CRISPR/Cas9 screen (29). Details on study designs and process of genes screening could refer to the previous research (8).

TABLE 1
www.frontiersin.org

Table 1 Clinical and pathologic features of eligible ccRCC patients from TCGA-KIRC and E-MTAB-1980.

FIGURE 1
www.frontiersin.org

Figure 1 (A) The casTLE score of ADCP-related genes in the CRISPR/Cas9 screen. (B) Clustering dendrogram and clinical traits heatmap of ccRCC patients based on ADCP-related genes. The shade of color in the heatmap increased with age, Fuhrman grade, pathology T stage, and pathology stage. Red, white, and gray are on behalf of metastasis, no metastasis, and data missing in pathology M and N stage, respectively. (C) The relationship between the soft threshold and the index of scale-free topologies and mean connectivity. (D) The heatmap of the correlation of gene modules and clinical traits. In each grid, the top number was the coefficient between gene modules and clinical traits, and the bottom was the p-value. (E–G) The functional clustering network of turquoise, blue, and brown modules. casTLE, cas9 high-throughput maximum likelihood estimator; ADCP, antibody-dependent cellular phagocytosis; ccRCC, clear cell renal cell carcinoma.

Construction of Co-Expression Modules and Identification of the Key Module

To find the gene modules with the highest correlation with clinical phenotypes, we implemented WGCNA analysis on TCGA-KIRC cohort. After deleting outliers and unqualified samples, 507 patients and 511 ADCP-related genes were enrolled into the subsequent analysis. The sample dendrogram and heatmap of clinical phenotypes are presented in Figure 1B. In order to construct the gene co-expression network, the optimal soft-thresholding value was set as 4 (scale-free R2 = 0.85, mean connectivity = 15.99, Figure 1C) and cut height was set as 0.5, after which the four-gene modules were identified (turquoise 369 genes, blue 97 genes, gray 6 genes, and brown 39 genes). Next, the correlations between gene modules and clinical traits were evaluated to find the key module, and the result showed that the blue module was most correlated with clinical features (Figure 1D). Then, we performed functional enrichment analysis to evaluate the biological function of each gene module based on the Metascape database. For the turquoise module, the top enriched terms were adaptive immune system, carbohydrate derivative biosynthetic process, and leukocyte differentiation (Figure 1E). For the blue module, lymphocyte activation, immune effector process, and positive regulation of immune response were significantly enriched (Figure 1F), which indicated that the blue module played important roles in immune regulation. Finally, the brown module was enriched in NADH dehydrogenase complex assembly, translation, and response to oxidative (Figure 1G).

Construction of the Gene Signature

To explore the prognostic value of ADCP-related genes, we constructed a gene signature based on TCGA-KIRC cohort. The workflow of the constructing gene signature is shown in Figure 2A. The genes in the blue module, which was most correlated with clinical phenotypes, were selected as candidate genes of the constructing gene signature. Taking advantage of univariate COX regression, we firstly screened 38 prognostic genes (p-value < 0.01, Supplementary Data Table S1) from the blue module. To avoid multicollinearity and overfitting, we then applied LASSO regression on these prognostic genes and recognized nine genes more related to prognosis in ccRCC (Figures 2B, C). Next, to further screen variates, stepwise regression was performed to identify the optimal model based on the principle of minimum AIC. Finally, FKBP11 (FKBP Prolyl Isomerase 11), CD1C (CD1c Molecule), PHF19 (PHD Finger Protein 19), ARPC3 (Actin Related Protein 2/3 Complex Subunit 3), MS4A14 (Membrane Spanning 4-Domains A14), and KDELR3 (KDEL Endoplasmic Reticulum Protein Retention Receptor 3) were incorporated into the multivariate COX regression model and risk score = (0.261 * FKBP11 expression) + (-0.627 * CD1C expression) + (0.667 * PHF19 expression) + (0.418 * ARPC3 expression) + (0.298 * MS4A14 expression) + (0.155 * KDELR3 expression) (Figure 2D). The risk score of each patient was calculated by this model, and each patient was divided into high- and low-risk groups according to the median of risk score.

FIGURE 2
www.frontiersin.org

Figure 2 Construction and validation of the gene signature. (A) The workflow of the constructing gene signature. (B) The coefficients of LASSO regression in the 10-fold cross validation based on 38 candidate genes. (C) Ten-fold cross-validation for tuning parameter selection. The dotted vertical lines are drawn at the optimal values by the minimum criteria and one-standard error criterion. (D) The forest plot of six hub genes based on univariate COX regression of overall survival in TCGA-KIRC cohort. (E, F) Kaplan–Meier survival curves of overall survival for patients in high- and low-risk groups in TCGA-KIRC cohort and E-MTAB-1980 cohort. (G, J) The distributions of the risk score and vital status in TCGA-KIRC cohort and E-MTAB-1980 cohort. (H, I) The time-dependent ROC curves at 1, 3, and 5 years based on the gene signature in TCGA-KIRC cohort and E-MTAB-1980 cohort. The proportion of high- and low-risk groups in different Fuhrman grades and pathology T stages in TCGA-KIRC cohort (K) and E-MTAB-1980 cohort (L) * 0.01 < p < 0.05; ** 0.001 < p < 0.01; ****p < 0.0001.

High-Risk Group Correlated With Worse Prognosis, Higher Fuhrman Grade, and Advanced Pathologic Stage

We then investigated the association between the gene signature and clinical phenotypes. In the training cohort, most of the dead cases were concentrated in the high-risk group (Figure 2G). Except for CD1C, other genes were upregulated in the high-risk group according to the heatmap (Figure 2G). The result of survival analysis demonstrated that the high-risk group had a worse prognosis than the low-risk group (p-value < 0.0001, Figure 2E). Additionally, the AUC (area under the curve) values of the gene signature at 1, 3, and 5 years were 0.784, 0.717, and 0.718, respectively (Figure 2H). Furthermore, with the elevation of grades and stages, the proportions of the high-risk group were significantly increased accordingly (Figure 2K). In the testing cohort, the above conclusions were consistent with the training cohort. More dead cases distributed in the high-risk group (Figure 2J). CD1C was overexpressed in the low-risk group, whereas the expressions of the others were upregulated in the high-risk group. The prognosis in the high-risk group was worse likewise (p-value = 0.013, Figure 2F). A higher predictive accuracy was presented in the E-MTAB-1980 cohort, and the AUC values at 1, 3, and 5 years were 0.798, 0.767, and 0.788 (Figure 2I). The distribution of the high- and low-risk groups in different stages and grades was consistent with the training cohort (Figure 2L). Finally, we investigated the relationship between six key genes and prognosis and clinical phenotypes. FKBP11, PHF19, ARPC3, MS4A14, and KDELR3 were significantly upregulated in M1 cases (Figure 3A). FKBP11, PHF19, ARPC3, and MS4A14 were overexpressed in tumor tissues (Figure 3B). Moreover, based on the optimal cutoff value of gene expression, higher expressions of FKBP11, PHF19, ARPC3, MS4A14, and KDELR3 had a worse prognosis in the training and testing cohorts, while patients with higher-level CD1C had a more prolonged OS in both cohorts (Figures 3C, D). The above conclusions indicate that the gene signature has a good performance in predicting prognosis and a robust correlation with clinical phenotypes. In addition, the gene signature is stable and verifiable.

FIGURE 3
www.frontiersin.org

Figure 3 The relationship between expression of six key genes and tumor, metastasis, and prognosis. The comparison expression of six hub genes in patients with and without distant metastasis (A), and tumor tissues and normal tissues (B). Kaplan–Meier survival curves of the overall survival of six hub genes in TCGA-KIRC cohort (C) and E-MTAB-1980 cohort (D) according to the optimal cutoff value of gene expression. * 0.01 < p < 0.05; ** 0.001 < p < 0.01; ****p < 0.0001; ns no significance.

The Gene Signature for Nomogram Construction

To quantitatively predict the prognosis of ccRCC, we integrated the gene signature and clinical phenotypes of patients from TCGA-KIRC cohort to construct a predictive nomogram. According to the median of risk score, patients were divided into high- and low-risk groups. Patients then were further subdivided into different subgroups: distant metastasis yes and no, lymph node metastasis yes and no, T1/T2, and T3/T4. The all-subset regression method was performed, and the model with maximal adjust R square was selected. Finally, risk score, age, distant metastasis, lymph node metastasis, pathology T stage, and shortest, intermediate, and longest-dimension of primary tumor were incorporated into multivariate COX regression analysis to construct the nomogram (Figure 4A). The overall survival at 1, 3, and 5 years for each patient could be predicted based on the total points from the nomogram. Furthermore, the time-dependent ROC curves presented that the nomogram had the highest AUC values at 1, 3, and 5 years of 0.911, 0.845, and 0.867, respectively, compared to other clinical features of ccRCC (Figures 4B–D). There was a good consistency between the actual observed OS and predictive OS by nomogram according to the calibration curves at 1, 3, and 5 years (Figures 4E–G). Taken together, the nomogram had a good performance on predictive accuracy in TCGA-KIRC cohort; however, there was no independent external cohort to validate it due to lack of tumor size data.

FIGURE 4
www.frontiersin.org

Figure 4 Construction of nomogram. (A) Nomogram to predict overall survival at 1, 3, and 5 years based on the gene signature, age, pathology T, N, and M stages, and tumor size. (B–D) Time-dependent ROC curves at 1, 3, and 5 years of the nomogram and clinical features of ccRCC. (E–G) Calibration curves of nomogram for predicting overall survival at 1, 3, and 5 years.

Somatic Variant Analysis Between High- and Low-Risk Groups

As somatic gene mutation could act as trunk events promoting tumorigenesis, especially in renal cell carcinoma, we figured out top 15 mutated genes in high- and low-risk groups separately (Figure 5A) and summarized variant classification displaying the number of variants in each sample as a stacked bar plot and variant types as a boxplot (Figure 5B). We identified the top five differentially mutated genes (Figure 5C) including VHL, PBRM1, TTN, BAP1, and SETD2 in high- and low-risk groups. The high-risk group with worse prognosis possessed a higher frequency of BAP1 mutation, and a previous study has proved that BAP1 loss is associated with high tumor grade (30). What is more, the frequency of PRBM1 mutation was higher in the low-risk group. A previous study has shown that PRBM1 is associated with angiogenesis and patients with PRBM1 mutation are more likely to benefit from antiangiogenics than from immunotherapy (31). We also compared tumor mutation burden (TMB) (Figure 5D) and copy number alteration (CNA) burden (Figure 5E), which have been considered as response biomarkers for immune checkpoint blockade in solid tumors such as melanoma, non-small cell lung cancer, and urothelial carcinoma. In our study, it was found that patients in the high-risk group had significantly higher TMB and CNA burden. According to the evidence that tumors with higher TMB implied better ICI efficacy (32), we speculated that the high-risk group with higher genomic alteration burden possessed a better response to anti-PD-1 or anti-PD-L1 antibody treatment in ccRCC.

FIGURE 5
www.frontiersin.org

Figure 5 Somatic variant analysis between high- and low-risk groups in TCGA. (A) Top 15 mutated genes in high- and low-risk groups. (B) Landscape of somatic mutations in TCGA-KIRC cohort. (C) Comparison of top five mutated genes in high- and low-risk groups. Differences in tumor mutation burden (D) and copy number alteration burden (E) between the high- and low-risk groups ** 0.001 < p < 0.01; ****p < 0.0001.

Functional Enrichment Analysis in High- and Low-Risk Groups

To further estimate the biological difference in distinct risk groups, we screened out differentially expressed genes (DEGs) between the high- and low-risk groups. A total of 393 genes were significantly upregulated in the high-risk group, while 214 genes were downregulated (Figure 6A). KEGG and GO analyses were performed based on 393 overexpressed genes in the high-risk group (Figures 6B, C). We found that the cytokine–cytokine receptor interaction and T-cell receptor signaling pathway in KEGG, and acute inflammatory response, regulation of T-cell activation, and T-cell activation in GO were enriched in the high-risk group. Similarly, GSEA showed that the high-risk group was significantly associated with several antitumor immune pathways (Figures 6D–I), which contained positive regulation of interleukin 17 production, positive regulation of interleukin 4 production, regulation of T helper I type immune response, response to interferon β, B cell-mediated immunity, and positive regulation of T-cell proliferation. To sum up, we revealed the biological differences between high- and low-risk groups mainly associated with immune-related pathways, and the high-risk group was significantly associated with antitumor immune pathways.

FIGURE 6
www.frontiersin.org

Figure 6 Functional enrichment between high- and low-risk groups. (A) Volcano plot of differentially expressed genes (DEGs) between high- and low-risk groups. Log2(fold change) > 1, adjusted p-value < 0.05. KEGG (B) and GO (C) enrichment analysis base on upregulated DEGs. Top 15 enriched pathways were presented. (D–I) GSEA between high- and low-risk groups. KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; GSEA, gene set enrichment analysis; NES, normalized enrichment score; NPval, nominal p-value; FDR qval, FDR q-value.

Immune Landscape of High- and Low-Risk Groups

We then evaluated tumor immune cell infiltration in the tumor microenvironment (TME) between high- and low-risk groups. The landscape of immune infiltration and clinical parameters displayed higher grade, advanced TNM stages, and more decreased events in the high-risk group (Figure 7A). CIBERSORT analysis was also used to assess the proportion of 22 lymphocytes in the TME of ccRCC (Figure 7B), and the result showed that the TME of the high-risk group was composed of more T-cell CD4+ memory activated/resting, T-cell CD4+ naïve, T-cell CD8+, and less macrophage M2. Specifically, the abundance of antitumor immune cells, including activated CD8 T cell, activated dendritic cell, central memory CD8 T cell, effector memory CD8 T cell, natural killer cell, and natural killer T cell, was higher in the high-risk group in TCGA and E-MTAB-1980 cohorts (Figures 7C, D). Besides, by using the ESTIMATE algorithm on TCGA-KIRC cohort, patients in the high-risk group had higher immune score, stromal score, and ESTIMATE score (Figures 7E–G).

FIGURE 7
www.frontiersin.org

Figure 7 Immune infiltrate landscape of high- and low-risk groups. (A) The landscape of immune infiltrates and clinical phenotypes between high- and low-risk groups. (B) The relative proportion of 22 immune cells in high- and low-risk groups by CIBERSORT analysis. Antitumor-related immune cells between high- and low-risk groups in TCGA-KIRC cohort (C) and E-MTAB-1980 cohort (D). (E–G) ESTIMATE analysis between high- and low-risk groups in TCGA-KIRC cohort. * 0.01 < p < 0.05; ** 0.001 < p < 0.01; *** 0.0001 < p < 0.001; ****p < 0.0001; ns no significance.

High-Risk Group Was Associated With a Better Response to Immunotherapy

To investigate whether the gene signature was related to the response to immunotherapy, we firstly performed IPS analysis on TCGA-KIRC cohort. In our study, the IPS-PD1/PDL1/PDL2 blocker score (PD-1/PD-L1/PD-L2 positive and CTLA-4 negative), IPS–CTLA4 blocker score (PD-1/PD-L1/PD-L2 negative and CTLA-4 positive), and IPS–CTLA4 and PD1/PDL1/PDL2 blocker score (PD-1/PD-L1/PD-L2 positive and CTLA-4 positive) were significantly higher in the high-risk group (Figure 8A), which indicated that patients with a high risk score might benefit from immunotherapy. Next, we explored expressions of the common inhibitory immune checkpoints (ICPs) in the high- and low-risk groups. We found that there were more inhibitory immune checkpoints (CTLA4, LAG3, and PD-1) significantly upregulated in the high-risk group (Figure 8B). In addition, it has been proved that the loss of the HLA gene family is associated with resistance to immunotherapy (33, 34). In our study, most of HLA genes were overexpressed in the high-risk group (Figure 8C), which further proved that patients in the high-risk group might be suitable for immunotherapy than the low-risk group. Moreover, in order to further elaborate the relationship between gene signature and response to immunotherapy, we directly validated it in a ccRCC cohort treated with the anti-PD-1 antibody. The ccRCC cohort consisted of 39 CR/PR (complete response/partial response) patients and 132 SD/PD (stable disease/progressive disease) patients. As is shown in Figure 8D, in the CR/PR group, patients with a higher risk score occupied more proportion (64.1% vs. 35.9%, p-value = 0.03). Taken together, patients with a higher risk score might be the optimal beneficiaries for immunotherapy and we provided a novel and potential gene signature to predict response to immunotherapy.

FIGURE 8
www.frontiersin.org

Figure 8 Correlation between the gene signature and response to immunotherapy. (A) The relationship between immunophenoscore and the gene signature. (B) Expression of inhibitory immune checkpoints in high- and low-risk groups. (C) Expression of the HLA gene family in high- and low-risk groups. (D) The association between risk groups and response to anti-PD-1 antibody treatment in the ccRCC cohort. CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease. * 0.01 < p < 0.05; ** 0.001 < p < 0.01; *** 0.0001 < p < 0.001; ****p < 0.0001; ns no significance.

Functional Analysis of Six Key Genes

To investigate the function of six key genes, we constructed the protein–protein interaction (PPI) networks evaluated by the STRING database. Within the PPI network, each node represented all the proteins produced by a single, protein-coding gene locus and each edge represented protein–protein associations. FKBP11 was mainly enriched in the biological process of the GO term associated with myoblast fusion involved in skeletal muscle regeneration, but the PPI network was of no statistical significance (Figure 9A). CD1C was enriched in several immune-related biological processes, such as positive regulation of interleukin-2 production (Figure 9B). PHF19 was mainly correlated with histone methylation (Figure 9C). ARPC3 was primarily involved in biological processes associated with meiosis (Figure 9D). However, MS4A14 was not enriched in any biological process of the GO term (Figure 9E). KDELR3 was associated with protein intracellular transport (Figure 9F). The detailed biological processes of GO terms about these genes are shown in Supplementary Data Table S2. In addition, we investigated the genetic dependencies of six key genes in cancer cells evaluated by the Cancer Dependency Map (DepMap) database. The Chronos dependency score, based on data from a cell depletion assay, was utilized to describe the importance of the gene of interest on a given cell line. A lower Chronos score indicated more essential to tumor cell proliferation. A score of 0 indicated a gene not essential; correspondingly, -1 means comparable to the median of all pan-essential genes. The results showed that only the medians of the Chronos score of ARPC3 were less than -1 in most tumor cells, including kidney cancer cells (Figure 9G), which suggested that ARPC3 was more essential to tumor cell proliferation.

FIGURE 9
www.frontiersin.org

Figure 9 Functional analysis of six key genes. (A–F) Protein–protein interaction networks of six key genes. (G) The Chronos dependency scores of ARPC3 gene in 31 tumor types.

Prediction for Response to Chemotherapy Drugs and Investigation of Potential Drugs Based on the Gene Signature

To select chemotherapy drugs suitable for patients in the high- and low-risk groups respectively, we predicted the IC50 values of common chemotherapy drugs in TCGA-KIRC based on the GDSC database. The results demonstrated that the IC50 values of pazopanib and sunitinib were lower in the low-risk group, whereas sorafenib, gefitinib, erlotinib, tivozanib, axitinib, and temsirolimus were lower in the high-risk group (Figure 10A). Based on the results, personalized treatments could be provided for each patient according to the risk score. Meanwhile, we utilized the cMAP database to explore small-molecule drugs which might be effective for ccRCC. Twelve small-molecule drugs were identified (p-value < 0.05, Supplementary Data Table S3) based on the 249 upregulated genes and 110 downregulated genes in the high-risk group. The 2D structures of the top 10 small-molecule drugs are presented on Figure 10B. Taken together, our gene signature may contribute to personalized treatment and the development of new drugs in ccRCC.

FIGURE 10
www.frontiersin.org

Figure 10 Exploring potential drugs. (A) Predicting half-maximal inhibitory concentration values of eight common chemotherapy drugs based on the GDSC database. (B) Identifying potential small-molecule drugs based on DEGs by the cMAP database. The 2D structures of top 10 small-molecule drugs were presented. IC50, half-maximal inhibitory concentration; GDSC database, Genomics of Drug Sensitivity in Cancer database; cMAP database, the connectivity map database. * 0.01 < p < 0.05; ** 0.001 < p < 0.01; *** 0.0001 < p < 0.001; ****p < 0.0001.

Validated Prognostic Markers and Correlation With TILs

Due to the availability and performance of IHC antibodies in the pilot assays, only FKBP11, ARPC3, KDELR3, and MS4A14 were validated for prognostics using our in-house FFPE samples with IHC staining (Figures 11A–D). A retrospective consecutive cohort comprising 72 samples from patients who underwent radical or partial nephrectomy from January 2009 to December 2011 was assembled for validation. We chose the patients from the last decade simply due to the relative latent disease course of ccRCC, and survival events were not mature in a more recent cohort. Demographic characteristics are shown in Table 2. Expressions of all four factors were associated with worsened overall survival and advanced tumor grade and stage, respectively (Figures 11A–D, Table 2). Of note, primary tumors in patients with metastasis (M1) also showed higher expressions of the four (Table 2). There was no gender disparity or age preference (Table 2 and Figure 11E). All four factors showed a moderate to strong correlation with TIL abundancy (Figure 11E).

FIGURE 11
www.frontiersin.org

Figure 11 Validation of ADCP-related genes in 72 in-house ccRCC samples. Shown were representative IHC images (“+” for score >=1 and “-” for score 0) and impact on overall survival of (A) FKBP11, (B) ARPC3, (C) KDELR3, and (D) MS4A14. (E) Representative H-score 1 and 0 for TIL abundance and heatmap for its correlations with ADCP-related genes and age (Spearman’s r). (Scale bar = 200 µm).

TABLE 2
www.frontiersin.org

Table 2 Association of with clinicopathological parameters of ccRCC patients in the IHC validation cohort.

Discussion

Thus far, there is only one study reporting ADCP in ccRCC and the process was used to evaluate the effect of the monoclonal antibody against CAIX (35). Our report is, to date, the first study comprehensively reporting the cancer-intrinsic ADCP-regulatory gene signature in ccRCC. Although the intercellular process of phagocytosis has been elucidated, the cancer-intrinsic genetic alteration that drives or expels ADCP remains largely undecided, until a recent study by Kamber et al. (8). They employed an unbiased CRISPR/Cas9 screen and identified a reservoir of novel cancer-intrinsic genes involved in ADCP. We took advantage of the gene set and further generated a six-gene signature specific in ccRCC which showed an association not only with clinicopathological parameters but also with response to ICI and TKI.

Our study identified six genes that are associated both with ADCP and with prognosis in ccRCC. ARPC3 showed the highest HR and has not been reported in ccRCC. The gene encodes one of seven subunits of the human Arp2/3 protein complex (36). The Arp2/3 protein complex has been conserved through evolution and is implicated in the control of actin polymerization in cells (37). The process of phagocytosis is highly complex and involves major rearrangements of the cytoskeleton, a process in which the Arp2/3 protein complex plays a role (38). Padilla et al. reported that mir-124-5p regulates phagocytosis of human macrophages by targeting the actin cytoskeleton via the Arp2/3 complex (39). Silencing of the ARP2/3 complex has also been shown to disturb pancreatic cancer cell migration (40). Taken together, ARPC3 could be of both prognostic and therapeutical value in ccRCC, which warrants further study.

In contrast, PHF19 is not reported in phagocytosis but has vastly been reported in cancer research, most prominently in malignant melanoma (41). PHF19 is the subunit of polycomb repressive complex 2 (PRC2) and regulates the expression of key genes involved in cell growth and differentiation by modulating both PRC2/EZH2 catalytic activity and recruitment (42). Of note, its role has not been reported in ccRCC either. FKBP11 belongs to the FK506-binding protein (FKBP) that can bind specifically to the immunosuppressant FK506 and rapamycin (43). Whereas FKBP11 has not been reported in phagocytosis, immunosuppression of FK506 is vastly applied clinically in organ transplantation (44). Notably, an in-silico study using datasets overlapping ours also showed a similar prognostic role of FKFB11 in ccRCC (45). Neither MS4A14 nor KDELR3 was reported related to phagocytosis, and the genes are functionally distant from ADCP in previous studies. Identification of such genes shed light on further understanding of phagocytosis in ccRCC. Interestingly, CD1C is the only protective gene shown in our signature. Prior studies of CD1C in ccRCC focused on its expression in TILs and in peripheral lymphocytes showing increased CD1C-positive cells upon treatment of sunitinib and bevacizumab (46, 47). In our cohort, TIMER (https://cistrome.shinyapps.io/timer/) was used to determine the distribution of CD1C in ccRCC and found that it was predominantly expressed in tumor cells in ccRCC (48, 49). Our finding that CD1C expression was higher in the low-risk group that corresponded to lower IC50 of sunitinib echoes the previous studies, supporting the signature as a therapeutic marker.

Our findings could have translational implications. Besides prognosis, our signature was also associated with drug sensitivity. The sensitivity profile to TKIs and ICI corresponds to the current guidelines. We speculate that the low-risk group may benefit more from sunitinib and pazopanib whereas the high-risk group may benefit from ICI. Of note, TKIs that fall into later lines of treatment in ccRCC show inversed sensitivity between high- and low-risk groups from sunitinib and pazopanib, further echoing the fact that axitinib, tivozanib, sorafenib, and temsirolimus cannot reach first-line therapy, possibly due to the inherent ADCP status of patients. Interestingly, the high-risk group that showed lower IC50 in axitinib and sensitivity to ICI further corroborates the recent trial supporting the front-line use of axitinib-based ICI+TKI combination (50, 51). This also implies that the combination of sunitinib or pazopanib with ICI may not be as synergistic. Last but not least, we have listed several agents that could modify the ADCP status based on our signature, which may direct further investigation.

To find the optimal biomarker to predict prognosis and drug response in renal cell carcinoma, many gene signatures have been constructed based on the same methods in recent years. Wu et al. (52), Yin et al. (53), and Zhang Z et al. (54) constructed gene signatures according to differential expression genes between tumor and normal tissues; Zhang et al. (55) and Lv et al. (56) established gene signatures based on tumor microenvironment-related genes and glycolysis-related genes, respectively. They investigated the relationship between gene signatures and clinical features, immune infiltration, and prognosis, and these gene signatures had a good performance on predicting prognosis. Notably, the gene signature from Yin et al. demonstrated a significant correlation with response to immune checkpoint inhibitors (ICIs); the gene signature from Lv et al. was significantly associated with response to ICIs and tyrosine kinase inhibitors (TKIs), which might be the potential biomarker to predict drug response. Compared with previous similar studies, our research focused on ADCP-related genes and emphasized the important effect of ADCP-related genes on prognosis, tumor microenvironment, response to immunotherapy, and targeted therapy in ccRCC. Our monogram, integrating our gene signature and clinical features, possessed higher accuracy of prediction, whose AUC at 1 year reached 0.911. Additionally, we not only investigated the relationship between the gene signature and response to ICI but also predicted response to several common TKIs, which might contribute to developing a personalized treatment plan for patients with ccRCC. However, the patients enrolled in the validation cohorts were generally at the early stage and TKI usage was not enough to undergo statistical analyses.

Our study has limitations, First, this is an in-silico analysis without external validation. A complete elaboration of ADCP in ccRCC still requires profound in vitro and in vivo studies. Second, some of the drug-sensitivity analyses are profiled using simulation. A more rigorous drug screen is now in progress by our group. Finally, just like Kamber’s study, the ADCP-related genes were selected upon a CRISPR screen and many of its functionality in phagocytosis remains elusive. How the cancer-intrinsic signaling format ADCP by those genes should be thoroughly studied.

Conclusion

Using an in-silico analysis by integrating multiple genomic datasets, we showed a six-gene ADCP signature that correlated with prognosis and immune modulation in ccRCC. The signature-based risk stratification was associated with response to both ICI and tyrosine kinase inhibition in ccRCC.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The studies involving human participants were reviewed and approved by the Huashan Institutional Review Board (HIRB). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

CF, YqL, XZ, and YfL carried out the in silico analysis. KL and CF participated in the study design. CF, HJ, and HW retrieved the data. CF and KL drafted the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was sponsored in part by the National Natural Science Foundation of China (Grant No. 81874123 and No. 81772709) and by the Beijing Advanced Innovation Center for Food Nutrition and Human Health, Beijing Technology and Business University (BTBU).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer FW declared a shared affiliation with the authors to the handling editor at the time of review.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.853088/full#supplementary-material

References

1. Harrison ML, Obermueller E, Maisey NR, Hoare S, Edmonds K, Li NF, et al. Tumor Necrosis Factor Alpha as a New Target for Renal Cell Carcinoma: Two Sequential Phase II Trials of Infliximab at Standard and High Dose. J Clin Oncol (2007) 25(29):4542–9. doi: 10.1200/JCO.2007.11.2136

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Pili R, Quinn DI, Hammers HJ, Monk P, George S, Dorff TB, et al. Immunomodulation by Entinostat in Renal Cell Carcinoma Patients Receiving High-Dose Interleukin 2: A Multicenter, Single-Arm, Phase I/II Trial (NCI-CTEP#7870). Clin Cancer Res (2017) 23(23):7199–208. doi: 10.1158/1078-0432.CCR-17-1178

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Andrews LP, Yano H, Vignali DAA. Inhibitory Receptors and Ligands Beyond PD-1, PD-L1 and CTLA-4: Breakthroughs or Backups. Nat Immunol (2019) 20(11):1425–34. doi: 10.1038/s41590-019-0512-0

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Bedke J, Albiges L, Capitanio U, Giles RH, Hora M, Lam TB, et al. The 2021 Updated European Association of Urology Guidelines on Renal Cell Carcinoma: Immune Checkpoint Inhibitor-Based Combination Therapies for Treatment-Naive Metastatic Clear-Cell Renal Cell Carcinoma Are Standard of Care. Eur Urol (2021) 80(4):393–7. doi: 10.1016/j.eururo.2021.04.042

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Ko JJ, Xie W, Kroeger N, Lee JL, Rini BI, Knox JJ, et al. The International Metastatic Renal Cell Carcinoma Database Consortium Model as a Prognostic Tool in Patients With Metastatic Renal Cell Carcinoma Previously Treated With First-Line Targeted Therapy: A Population-Based Study. Lancet Oncol (2015) 16(3):293–300. doi: 10.1016/S1470-2045(14)71222-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Yi M, Jiao D, Qin S, Chu Q, Wu K, Li A. Synergistic Effect of Immune Checkpoint Blockade and Anti-Angiogenesis in Cancer Treatment. Mol Cancer (2019) 18(1):60. doi: 10.1186/s12943-019-0974-6

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Yuan J, Khilnani A, Brody J, Andtbacka RHI, Hu-Lieskovan S, Luke JJ, et al. Current Strategies for Intratumoural Immunotherapy - Beyond Immune Checkpoint Inhibition. Eur J Cancer (2021) 157:493–510. doi: 10.1016/j.ejca.2021.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Kamber RA, Nishiga Y, Morton B, Banuelos AM, Barkal AA, Vences-Catalán F, et al. Inter-Cellular CRISPR Screens Reveal Regulators of Cancer Cell Phagocytosis. Nature (2021) 597(7877):549–54. doi: 10.1038/s41586-021-03879-4

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Cunha LD, Yang M, Carter R, Guy C, Harris L, Crawford JC, et al. LC3-Associated Phagocytosis in Myeloid Cells Promotes Tumor Immune Tolerance. Cell (2018) 175(2):429–41e16. doi: 10.1016/j.cell.2018.08.061

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Su S, Zhao J, Xing Y, Zhang X, Liu J, Ouyang Q, et al. Immune Checkpoint Inhibition Overcomes ADCP-Induced Immunosuppression by Macrophages. Cell (2018) 175(2):442–57e23. doi: 10.1016/j.cell.2018.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Motzer R, Alekseev B, Rha SY, Porta C, Eto M, Powles T, et al. Lenvatinib Plus Pembrolizumab or Everolimus for Advanced Renal Cell Carcinoma. N Engl J Med (2021) 384(14):1289–300. doi: 10.1056/NEJMoa2035716

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Langfelder P, Horvath S. WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

13. Budczies J, Seidel A, Christopoulos P, Endris V, Kloor M, Győrffy B, et al. Integrated Analysis of the Immunological and Genetic Status in and Across Cancer Types: Impact of Mutational Signatures Beyond Tumor Mutational Burden. Oncoimmunology (2018) 7(12):e1526613. doi: 10.1080/2162402X.2018.1526613

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Lu Z, Chen H, Li S, Gong J, Li J, Zou J, et al. Tumor Copy-Number Alterations Predict Response to Immune-Checkpoint-Blockade in Gastrointestinal Cancer. J Immunother Cancer (2020) 8(2):e000374. doi: 10.1136/jitc-2019-000374

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-Cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Choueiri TK, Fishman MN, Escudier B, McDermott DF, Drake CG, Kluger H, et al. Immunomodulatory Activity of Nivolumab in Metastatic Renal Cell Carcinoma. Clin Cancer Res (2016) 22(22):5461–71. doi: 10.1158/1078-0432.CCR-15-2839

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Motzer RJ, Rini BI, McDermott DF, Redman BG, Kuzel TM, Harrison MR, et al. Nivolumab for Metastatic Renal Cell Carcinoma: Results of a Randomized Phase II Trial. J Clin Oncol (2015) 33(13):1430–7. doi: 10.1200/JCO.2014.59.0703

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Motzer RJ, Tannir NM, McDermott DF, Arén Frontera O, Melichar B, Choueiri TK, et al. Nivolumab Plus Ipilimumab Versus Sunitinib in Advanced Renal-Cell Carcinoma. N Engl J Med (2018) 378(14):1277–90. doi: 10.1056/NEJMoa1712126

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant' Angelo M, Forman J, et al. Interplay of Somatic Alterations and Immune Infiltration Modulates Response to PD-1 Blockade in Advanced Clear Cell Renal Cell Carcinoma. Nat Med (2020) 26(6):909–18. doi: 10.1038/s41591-020-0839-y

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Geeleher P, Cox N, Huang RS. Prrophetic: An R Package for Prediction of Clinical Chemotherapeutic Response From Tumor Gene Expression Levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, et al. The Connectivity Map: Using Gene-Expression Signatures to Connect Small Molecules, Genes, and Disease. Science (2006) 313(5795):1929–35. doi: 10.1126/science.1132939

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Tissue-Based Map of the Human Proteome. Science (2015) 347(6220):1260419. doi: 10.1126/science.1260419

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Feng C, Ding G, Ding Q, Wen H. Overexpression of Low Density Lipoprotein Receptor-Related Protein 1 (LRP1) Is Associated With Worsened Prognosis and Decreased Cancer Immunity in Clear-Cell Renal Cell Carcinoma. Biochem Biophys Res Commun (2018) 503(3):1537–43. doi: 10.1016/j.bbrc.2018.07.076

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Morgens DW, Deans RM, Li A, Bassik MC. Systematic Comparison of CRISPR/Cas9 and RNAi Screens for Essential Genes. Nat Biotechnol (2016) 34(6):634–6. doi: 10.1038/nbt.3567

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Pena-Llopis S, Vega-Rubin-de-Celis S, Liao A, Leng N, Pavía-Jiménez A, Wang S, et al. BAP1 Loss Defines a New Class of Renal Cell Carcinoma. Nat Genet (2012) 44(7):751–9. doi: 10.1038/ng.2323

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Carril-Ajuria L, Santos M, Roldan-Romero JM, Rodriguez-Antona C, de Velasco G. Prognostic and Predictive Value of PBRM1 in Clear Cell Renal Cell Carcinoma. Cancers (Basel) (2019) 12(1):16. doi: 10.3390/cancers12010016

CrossRef Full Text | Google Scholar

32. Yarchoan M, Hopkins A, Jaffee EM. Tumor Mutational Burden and Response Rate to PD-1 Inhibition. N Engl J Med (2017) 377(25):2500–1. doi: 10.1056/NEJMc1713444

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Flores-Martin JF, Perea F, Exposito-Ruiz M, Carretero FJ, Rodriguez T, Villamediana M, et al. A Combination of Positive Tumor HLA-I and Negative PD-L1 Expression Provides an Immune Rejection Mechanism in Bladder Cancer. Ann Surg Oncol (2019) 26(8):2631–9. doi: 10.1245/s10434-019-07371-2

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Garrido F. HLA Class-I Expression and Cancer Immunotherapy. Adv Exp Med Biol (2019) 1151:79–90. doi: 10.1007/978-3-030-17864-2_3

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Chang DK, Moniz RJ, Xu Z, Sun J, Signoretti S, Zhu Q, et al. Human Anti-CAIX Antibodies Mediate Immune Cell Inhibition of Renal Cell Carcinoma In Vitro and in a Humanized Mouse Model In Vivo. Mol Cancer (2015) 14:119. doi: 10.1186/s12943-015-0384-3

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Machesky LM, Reeves E, Wientjes F, Mattheyse FJ, Grogan A, Totty NF, et al. Mammalian Actin-Related Protein 2/3 Complex Localizes to Regions of Lamellipodial Protrusion and Is Composed of Evolutionarily Conserved Proteins. Biochem J (1997) 328( Pt 1):105–12. doi: 10.1042/bj3280105

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Welch MD, DePace AH, Verma S, Iwamatsu A, Mitchison TJ. The Human Arp2/3 Complex is Composed of Evolutionarily Conserved Subunits and is Localized to Cellular Regions of Dynamic Actin Filament Assembly. J Cell Biol (1997) 138(2):375–84. doi: 10.1083/jcb.138.2.375

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Rotty JD, Brighton HE, Craig SL, Asokan SB, Cheng N, Ting JP, et al. Arp2/3 Complex Is Required for Macrophage Integrin Functions But Is Dispensable for FcR Phagocytosis and In Vivo Motility. Dev Cell (2017) 42(5):498–513 e6. doi: 10.1016/j.devcel.2017.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Herdoiza Padilla E, Crauwels P, Bergner T, Wiederspohn N, Förstner S, Rinas R, et al. Mir-124-5p Regulates Phagocytosis of Human Macrophages by Targeting the Actin Cytoskeleton via the ARP2/3 Complex. Front Immunol (2019) 10:2210. doi: 10.3389/fimmu.2019.02210

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Rauhala HE, Teppo S, Niemelä S, Kallioniemi A. Silencing of the ARP2/3 Complex Disturbs Pancreatic Cancer Cell Migration. Anticancer Res (2013) 33(1):45–52.

PubMed Abstract | Google Scholar

41. Ghislin S, Deshayes F, Middendorp S, Boggetto N, Alcaide-Loridan C. PHF19 and Akt Control the Switch Between Proliferative and Invasive States in Melanoma. Cell Cycle (2012) 11(8):1634–45. doi: 10.4161/cc.20095

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ji Y, Fioravanti J, Zhu W, Wang H, Wu T, Hu J, et al. miR-155 Harnesses Phf19 to Potentiate Cancer Immunotherapy Through Epigenetic Reprogramming of CD8+ T Cell Fate. Nat Commun (2019) 10(1):2157. doi: 10.1038/s41467-019-09882-8

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Pinto D, Duarte M, Soares S, Tropschug M, Videira A. Identification of All FK506-Binding Proteins From Neurospora Crassa. Fungal Genet Biol (2008) 45(12):1600–7. doi: 10.1016/j.fgb.2008.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Vaughn CJ. Drugs and Lactation Database: LactMed. J Electron Resou Med Libraries (2012) 9:272–7. doi: 10.1080/15424065.2012.735134

CrossRef Full Text | Google Scholar

45. Sun Z, Qin X, Fang J, Tang Y, Fan Y. Multi-Omics Analysis of the Expression and Prognosis for FKBP Gene Family in Renal Cancer. Front Oncol (2021) 11:697534. doi: 10.3389/fonc.2021.697534

PubMed Abstract | CrossRef Full Text | Google Scholar

46. van Cruijsen H, van der Veldt AA, Vroling L, Oosterhoff D, Broxterman HJ, Scheper RJ, et al. Sunitinib-Induced Myeloid Lineage Redistribution in Renal Cell Cancer Patients: CD1c+ Dendritic Cell Frequency Predicts Progression-Free Survival. Clin Cancer Res (2008) 14(18):5884–92. doi: 10.1158/1078-0432.CCR-08-0656

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Garcia JA, Mekhail T, Elson P, Triozzi P, Nemec C, Dreicer R, et al. Clinical and Immunomodulatory Effects of Bevacizumab and Low-Dose Interleukin-2 in Patients With Metastatic Renal Cell Carcinoma: Results From a Phase II Trial. BJU Int (2011) 107(4):562–70. doi: 10.1111/j.1464-410X.2010.09573.x

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res (2017) 77(21):e108–e10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, et al. Comprehensive Analyses of Tumor Immunity: Implications for Cancer Immunotherapy. Genome Biol (2016) 17(1):174. doi: 10.1186/s13059-016-1028-7

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Rini BI, Plimack ER, Stus V, Gafanov R, Hawkins R, Nosov D, et al. Pembrolizumab Plus Axitinib Versus Sunitinib for Advanced Renal-Cell Carcinoma. N Engl J Med (2019) 380(12):1116–27. doi: 10.1056/NEJMoa1816714

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Motzer RJ, Penkov K, Haanen J, Rini B, Albiges L, Campbell MT, et al. Avelumab Plus Axitinib Versus Sunitinib for Advanced Renal-Cell Carcinoma. N Engl J Med (2019) 380(12):1103–15. doi: 10.1056/NEJMoa1816047

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Wu X, Zhao Z, Khan A, Cai C, Lv D, Gu D, et al. Identification of a Novel Signature and Construction of a Nomogram Predicting Overall Survival in Clear Cell Renal Cell Carcinoma. Front Genet (2020) 11:1017. doi: 10.3389/fgene.2020.01017

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Yin X, Wang Z, Wang J, Xu Y, Kong W, Zhang J. Development of a Novel Gene Signature to Predict Prognosis and Response to PD-1 Blockade in Clear Cell Renal Cell Carcinoma. Oncoimmunology (2021) 10(1):1933332. doi: 10.1080/2162402X.2021.1933332

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Zhang Z, Lin E, Zhuang H, Xie L, Feng X, Liu J, et al. Construction of a Novel Gene-Based Model for Prognosis Prediction of Clear Cell Renal Cell Carcinoma. Cancer Cell Int (2020) 20:27. doi: 10.1186/s12935-020-1113-6

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Zhang L, Li J, Zhang M, Wang L, Yang T, Shao Q, et al. Identification of a Six-Gene Prognostic Signature Characterized by Tumor Microenvironment Immune Profiles in Clear Cell Renal Cell Carcinoma. Front Genet (2021) 12:722421. doi: 10.3389/fgene.2021.722421

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Lv Z, Qi L, Hu X, Mo M, Jiang H, Li Y. Identification of a Novel Glycolysis-Related Gene Signature Correlates With the Prognosis and Therapeutic Responses in Patients With Clear Cell Renal Cell Carcinoma. Front Oncol (2021) 11:633950. doi: 10.3389/fonc.2021.633950

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: clear-cell renal cell carcinoma, antibody-dependent cellular phagocytosis, immune checkpoint inhibition, biomarker, bioinformatics

Citation: Li K, Li Y, Lyu Y, Tan L, Zheng X, Jiang H, Wen H and Feng C (2022) Development of a Phagocytosis-Dependent Gene Signature to Predict Prognosis and Response to Checkpoint Inhibition in Clear-Cell Renal Cell Carcinoma. Front. Immunol. 13:853088. doi: 10.3389/fimmu.2022.853088

Received: 12 January 2022; Accepted: 12 April 2022;
Published: 16 May 2022.

Edited by:

Alejandro López-Soto, University of Oviedo, Spain

Reviewed by:

Consuelo Amantini, University of Camerino, Italy
Fangning Wan, Fudan University, China
Xiaming Liu, Huazhong University of Science and Technology, China

Copyright © 2022 Li, Li, Lyu, Tan, Zheng, Jiang, Wen and Feng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Hui Wen, drwenhui@163.com; Haowen Jiang, drjianghaowen@163.com; Chenchen Feng, drfengchenchen@163.com

These authors have contributed equally to this work

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.