Article Text
Abstract
Background Recent studies have focused on the correlation between N6-methyladenosine (m6A) modification and specific tumor-infiltrating immune cells. However, the potential roles of m6A modification in the tumor immune landscape remain elusive.
Methods We comprehensively evaluated the m6A modification patterns and tumor immune landscape of 513 clear cell renal cell carcinoma (ccRCC) patients, and correlated the m6A modification patterns with the immune landscape. The m6Ascore was established using principal component analysis. Multivariate Cox regression analysis was performed to evaluate the prognostic value of the m6Ascore.
Results We identified three m6Aclusters—characterized by differences in Th17 signature, extent of intratumor heterogeneity, overall cell proliferation, aneuploidy, expression of immunomodulatory genes, overall somatic copy number alterations, and prognosis. The m6Ascore was established to quantify the m6A modification pattern of individual ccRCC patients. Further analyses revealed that the m6Ascore was an independent prognostic factor of ccRCC. Finally, we verified the prognostic value of the m6Ascore in the programmed cell death protein 1 (PD-1) blockade therapy of patients with advanced ccRCC.
Conclusions This study demonstrated the correlation between m6A modification and the tumor immune landscape in ccRCC. The comprehensive evaluation of m6A modification patterns in individual ccRCC patients enhances our understanding of the tumor immune landscape and provides a new approach toward new and improved immunotherapeutic strategies for ccRCC patients.
- immunotherapy
- kidney neoplasms
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
Background
Methylation of N6 adenosine to produce N6-methyladenosine (m6A) is the most common type of RNA modification1; it is thought to regulate multiple RNA-related processes, such as RNA stability,2 translation,3 alternative splicing4 5 and nuclear export.6 m6A modification is a dynamic and reversible process which is regulated by m6A methyltransferases (‘writers’), m6A demethylases (‘erasers’) and m6A-binding proteins (‘readers’).7 The m6A methyltransferases—consisting of METTL3, METTL14, WTAP, RBM15, RBM15B, ZC3H13, CBLL1 and VIRMA—catalyze m6A modification as m6A writers, while a set of m6A demethylases— including ALKBH5 and FTO—mediate the reversal of m6A modification of RNA as m6A erasers.8 9 Moreover, m6A-binding proteins—such as IGF2BP1/2/3, YTHDF1/2/3 and YTHDC1/2—recognize and bind to the m6A methylation sites in RNA as m6A readers.8 9 m6A is an essential RNA modification that regulates multiple key cellular processes including cellular differentiation, stem cell renewal and response to DNA damage.10 Evidently, aberrant expression of m6A regulators is associated with tumorigenesis, malignant tumor progression and immunomodulatory abnormality.10 11
Immune checkpoint therapy (ICT)—such as programmed cell death protein 1 (PD-1)/PD ligand 1 (PD-L1) blockade therapy—is transformative in the treatment of advanced clear cell renal cell carcinoma (ccRCC).12 13 However, there remains a considerable proportion of patients with no response or resistance to ICT.14 In solid malignant tumors, the PD-1 blocking response is associated with numerous tumor-intrinsic15 16 and tumor immune microenvironment (TIME) characteristics.17 18 A common paradigm in the immunology of solid tumors is that effective responses to anti-PD-1 therapy occur when the TIME is characterized by high infiltration of CD8+ T cells and that resistance to this therapy occurs when the TIME is characterized by the lack of such an infiltration.19 20 Understanding the biology of the TIME that drives the ICT response is crucial to the design of immunotherapeutic strategies.21 22
Several studies have recently focused on the special relationship between m6A regulators and immune cells. Wang et al23 reported that METTL3-mediated m6A modification increased the translation of certain immune transcripts and physiologically promoted the activation of dendritic cells (DCs) and DC-based T-cell responses. Li et al24 showed that deletion of METTL3 in T cells disrupted the homeostasis and differentiation of T cells. Han et al25 found that deletion of YTHDF1 elevated the antitumor response of antigen-specific CD8+ T cells and enhanced the efficacy of anti-PD-L1 therapy. However, limited by existing experimental technology, the above research is confined to a limited number of m6A regulators and cell types, while the development and progression of cancers depend on cross-talk among multiple m6A regulators of RNA methylation.9 Therefore, a comprehensive evaluation of the immune landscape mediated by a variety of m6A regulators will enhance our overall understanding of the immunomodulatory effect of m6A regulators on the TIME. Recently, the m6A modification patterns of gastric cancer were comprehensively evaluated based on multiple m6A regulators and systematically correlated with the tumor immune landscape, indicating the important role of m6A modification in TIME diversity in gastric cancer.26
In the present study, we integrated the molecular and clinical data of 513 ccRCC patients to comprehensively evaluate the m6A modification patterns and tumor immune landscape and correlated the m6A modification patterns with the immune landscape. We identified three distinct m6A modification patterns and were surprised to find that they had distinct immune landscapes and prognoses, indicating the crucial roles of m6A modification in the formation of individual tumor immune landscapes in ccRCC patients. We went on to quantify the m6A modification patterns of individual ccRCC patients by establishing the gene signature of m6A regulators.
Methods
Molecular and clinical data
The workflow of our study is shown in online supplemental figure S1. RNA sequencing data (count values) for gene expression analysis, genetic mutations (VarScan), and clinical data were downloaded from the Genomic Data Commons (https://portal.gdc.cancer.gov/).27 The count values were transformed into transcripts per kilobase million (TPM) values (the gene lengths used for the above transformation were measured as total non-overlapping exon length) and the Ensembl gene IDs of the RNA-seq data were converted to gene symbols by referring to the annotation file (https://www.gencodegenes.org/human/release_22.html). The copy number variation (CNV) data were downloaded from the Broad GDAC Firehose (https://gdac.broadinstitute.org/). The normalized data from another ccRCC cohort (91 cases) were downloaded from the International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/).
Supplemental material
Model-based clustering analysis for m6A regulators
Gene expression levels were quantified using the metric log2 (TPM +1), then used to identify m6A modification patterns based on the expression of 24 m6A regulators by model-based clustering analysis implemented in the R package ‘mclust’.28 The optimal number of clusters was determined based on the Bayesian information criterion.
Immune cellular fraction estimates
CIBERSORT—a deconvolution algorithm reported by Newman et al29 and verified by fluorescence-activated cell sorting—was used to quantify the 22 infiltrated immune cells according to normalized gene expression profiles. The 22 immune cells included memory B cells, naïve B cells, plasma cells, resting/activated DCs, resting/activated natural killer (NK) cells, resting/activated mast cells, eosinophils, neutrophils, monocytes, M0–M2 macrophages, and seven T-cell types (CD8+ T cells, regulatory T cells (Tregs), resting/activated memory CD4+ T cells, follicular helper T cells, naïve CD4+ T cells and gammadelta T cells (γδ T cells)) For each sample, the sum of all estimated values for the proportion of immune cells was equal to 1. CIBERSORT results were obtained from the following website: https://gdc.cancer.gov/about-data/publications/panimmune (online supplemental table S1).30 The relative abundance of Th1/Th2/Th17 cell infiltration in the ccRCC TIME was quantified by single-sample gene-set enrichment analysis. The gene sets for marking the Th1/Th2/Th17 cell types were obtained from a study published by Thorsson et al.30 The prognostic value of infiltrated immune cells was assessed by univariate Cox regression analysis.
Supplemental material
Evaluation of values of key immune characteristics and measures of DNA damage among m6Aclusters
Values of key immune characteristics (including leukocyte fraction, Th1/Th2/Th17 cells, single nucleotide variant neoantigens, indel neoantigens, proliferation, aneuploidy score, intratumor heterogeneity (ITH), B-cell receptor (BCR) evenness, T-cell receptor (TCR) evenness and cancer testis antigens (CTA) score) and measures of DNA damage (including CNV burden (number of segments and fraction of genome alterations, respectively), loss of heterozygosity (LOH; number of segments with LOH events, and fraction of bases with LOH events, respectively), homologous recombination deficiency, and mutation load (non-silent mutation)) were obtained from the following website: https://gdc.cancer.gov/about-data/publications/panimmune (online supplemental table S1).
Correlations between the expression characteristics of m6A regulators and immunomodulators
A list of 78 immunomodulators (IMs) was obtained (online supplemental table S2),30 three of which (HLA-DRB3, HLA-DRB4 and KIR2DL2) had no corresponding mRNA expression and were excluded from subsequent analysis. The median expression levels of the samples were used to represent the expression of each ccRCC subtype. In order to examine the differences in IM expression among different subtypes, we carried out the Kruskal-Wallis test on the gene expression levels for each ccRCC subtype. CNVs for each IM gene were obtained from the following website: https://gdc.cancer.gov/about-data/publications/panimmune.30 We calculated the difference between the observed and expected amplification frequencies (deletions) for each IM gene in each ccRCC subtype, where the expected frequency is the overall amplification frequency (deletions) of all ccRCC cases.
Gene set variation analysis
Gene set variation analysis (GSVA)—a non-parametric and unsupervised method commonly used for estimating pathway variations in the samples of expression datasets—was performed to explore the differences in biological processes among m6A modification patterns.31 The ‘c2.cp.kegg.v6.2.symbols’ gene sets for GSVA were downloaded from the Molecular Signatures Database (MSigDB). A p<0.05 was considered statistically significant.
Identification of differentially expressed genes among m6Aclusters
To identify genes related to m6A modification patterns, we classified patients into m6Aclusters based on the expression of 24 m6A regulators. Differentially expressed genes (DEGs) among these clusters were determined using the R package ‘DESeq2’, which was applied using the raw count values of RNA sequencing data. Genes with adjusted p<0.01 and at least two-fold changes in expression were identified as DEGs.
Construction of the m6A gene signature
We applied a methodology to quantify the m6A modification pattern (m6Ascore) of individual ccRCC patients. The m6Ascore was established as follows. First, we extracted the overlapping DEGs among m6Aclusters and classified the ccRCC patients into several groups using model-based clustering to analyze overlapping DEGs. Univariate Cox regression analysis was performed to evaluate the prognosis of each overlapping DEG. Genes with a significant prognosis (p<0.05) were extracted for further analysis. Next, principal component analysis (PCA) was performed to establish the m6A gene signature. We selected both principal components 1 and 2 as signature scores. Finally, the m6Ascore was defined using a method similar to Genomic Grade Index26 32 33:
where i is the expression of overlapping genes with a significant prognosis of DEGs among m6Aclusters.
Correlation between m6Ascore and other relevant biological processes
Spearman’s correlation analysis was performed to investigate the correlation between m6Ascore and other relevant biological processes using the gene sets reported by Mariathasan et al (online supplemental table S3),18 including (1) antigen processing machinery (APM), (2) effector CD8 T-cell signature, (3) immune checkpoint, (4) nucleotide excision repair, (5) mismatch repair, (6) DNA replication, (7) DNA damage repair, (8) epithelial-mesenchymal transition markers, (9) Wnt targets, (10) pan-fibroblast transforming growth factor-β response signature, and (11) angiogenesis signature.
Genomic and clinical data with anti-PD-1 therapy for ccRCC
A systematic search for the genomic and transcriptomic datasets of ccRCC patients treated with anti-PD-1 therapy was performed. We ultimately included one immunotherapeutic cohort—advanced ccRCC with treatment of PD-1 blockade and mammalian target of rapamycin (mTOR) inhibition—obtaining the genomic, transcriptomic and clinical data (online supplemental table S4) from the online supplemental data appended to the published paper.34
Supplemental material
Statistical analysis
Statistical significance for three or more groups was estimated using the Kruskal-Wallis test and association between categorical variables was explored using the χ2 test. The correlation coefficient was calculated via Spearman’s correlation analysis. Continuous variables were dichotomized for patient survival using optimal cut-off values determined by ‘survminer’ R package. The Kaplan-Meier method was used to generate survival curves and the log-rank test was used to determine the statistical significance of differences. The independent prognostic factors, determined by multivariate Cox regression analysis, were visualized by ‘forestplot’ R package. The ‘oncoplot’ function of R package ‘maftools’ was used to depict the mutation landscape of The Cancer Genome Atlas Kidney Renal Clear Cell Carcinoma (TCGA-KIRC) cohort and immunotherapeutic cohort. The protein–protein interaction (PPI) networks among m6A regulators were identified based on the STRING interaction database35 and visualized by Cytoscape.36 All tests were two sided, and p<0.05 was regarded as significant. All analyses were performed with R software V.3.62 (http://www.R-project.org).
Results
Molecular characteristics and clinical relevance of m6A regulators in ccRCC
On reviewing the literature, we identified 24 genes that mainly regulate RNA methylation including 8 writers (RBM15/RBM15B, METTL14, METTL3, WTAP, CBLL1, VIRMA and ZC3H13), 2 erasers (FTO and ALKBH5) and 14 readers (FMR1, ELAVL1, HNRNPC, HNRNPA2B1, YTHDF1/2/3, YTHDC1/2, RBMX, IGF2BP1/2/3 and LRPPPRC). Somatic mutations and CNVs were integrated to explore the prevalence of m6A regulator variations in ccRCC. The overall average mutation frequency of m6A regulators was low, with only 27 of 336 samples having m6A regulator mutations (figure 1A). We then studied the CNV alteration frequency of the m6A regulators and demonstrated that CNV alterations were prevalent (figure 1B). The mRNA expression levels of m6A regulators in ccRCC and adjacent tissues were also explored, revealing that 21 out of 24 m6A regulators were differentially expressed (figure 1C). The above analyses showed that the genetic and expressional variations in m6A regulators were highly heterogeneous between ccRCC and adjacent tissues, suggesting a crucial role for the imbalance of m6A regulator expression in the development and progression of ccRCC. Moreover, the function of genes is not isolated, in that it has been shown that collaboration among m6A regulators exists in the context of cancer.37 38 Thus, the correlation of mRNA expression among m6A regulators was explored. We identified that writers, erasers, and readers had a high expression correlation (figure 1D, (online supplemental table S5) and interacted with each other frequently in PPI networks (online supplemental figure S2). Taken together, these results indicate crucial cross-talk roles among m6A regulators of RNA methylation in the formation of distinct m6A modification patterns.
Supplemental material
Next, the clinical relevance of m6A regulators in ccRCC patients was explored. We found that many m6A regulators were related to prognosis in patients with ccRCC (figure 1D). Several m6A regulators (eg, IGF2BP1 and IGF2BP3) presented oncogenic characteristics, with higher expression levels of these genes related to poor prognosis in ccRCC patients. In contrast, we found that several m6A regulators (eg, LRPPRC and METTL14) presented characteristics of tumor suppressors, with higher expression levels of these genes correlated with favorable prognosis in ccRCC patients.
m6A modification patterns mediated by 24 m6A regulators
Model-based clustering was performed to classify ccRCC patients based on the expression of 24 m6A regulators. We ultimately uncovered three distinct methylation modification patterns (identified as m6Aclusters C1–C3), including 118 cases in m6Acluster-C1, 110 cases in m6Acluster-C2, and 285 cases in m6Acluster-C3 (figure 2A). The expression of m6A regulators with the greatest differences among subtypes (p<10−15) included two risk factors for overall survival (OS) (IGF2BP2 and IGF2BP3) and four favorable factors for OS (YTHDC1, RBMX, METTL14 and FTO) (online supplemental figure S3). m6Acluster-C3 was characterized by low expression levels of IGF2BP2 and IGF2BP3 and high expression levels of YTHDC1, RBMX, METTL14 and FTO (online supplemental figure S3). Therefore, it was not surprizing that m6Acluster-C3 had the most favorable prognosis (figure 2B).
Supplemental material
GSVA was performed to investigate the activity of biological processes among these distinct m6A modification patterns. As shown in figure 2C–D and online supplemental table S6, m6Aclusters-C1 and -C3 were markedly enriched in pathways related to immune activation, including the activation of cytokine–cytokine receptor interaction and the chemokine signaling pathway, TCR signaling pathway and BCR signaling pathway. Meanwhile, m6Acluster-C1 showed enrichment in stromal activation pathways such as cell adhesion and ECM receptor interaction. In contrast, m6Acluster-C2 was predominantly associated with the biological process of immunosuppression.
Immune characteristics in distinct m6A modification patterns
Thorsson et al30 explored the pan-cancer immune landscape and finally identified six immune subtypes (C1–C6) covering 30 cancer types that were assumed to define immune response patterns with implications for further exploration of immunotherapy. Immune subtype C3—characterized by elevated Th17, low to moderate tumor cell proliferation, and lower levels of overall CNVs and aneuploidy than the other immune subtypes—was enriched in most ccRCC patients. Strikingly, the three distinct methylation modification patterns had distinct proportions of the C3 immune subtype, with m6Acluster-C3 having the highest (96.14%), followed by m6Acluster-C1 (90.68%) and C2 (57.27%) (p<0.001) (figure 2A). We then explored the detailed immune characteristics in distinct m6A modification patterns. As shown in figure 3A–C, m6Acluster-C1 had a high proliferation rate, ITH, and lower levels of aneuploidy and overall CNVs. m6Acluster-C2 had the highest aneuploidy score and overall CNVs, as well as a high proliferation rate and ITH, and presented a more prominent macrophage signature dominated by M0 macrophages (online supplemental figure S4). m6Acluster-C3 was defined by elevated Th17, low tumor cell proliferation, ITH and lower levels of aneuploidy and overall CNVs.
Supplemental material
IMs are essential for cancer immunotherapy with multiple IM agonists and antagonists being investigated in clinical oncology.39 To advance this research, an understanding of their expression in different m6A modification patterns is needed. We explored IM gene expression and CNVs among the m6A subtypes (figure 3D). IM gene expression varied across m6A subtypes and genes with the greatest differences between subtypes (p<10−15) including ADORA2A, CX3CL1, EDNRB, ENTPD1, HMGB1, TNFRSF4, VEGFA and C10orf54 were most highly expressed in m6Acluster-C3 (online supplemental figure S5). CNVs affected numerous IMs and varied across m6A subtypes. m6Acluster-C1 and C2 exhibited frequent amplification and deletion of most IM genes.
Supplemental material
Immune landscape was significantly associated with the expression of m6A regulators
Spearman’s correlation analysis was performed to explore the specific correlation between each m6A regulator and immune cell infiltration. As shown in online supplemental figure S6, there was a widespread correlation between the expression of m6A regulators and immune cell infiltration. We focused on the regulator IGF2BP3—an m6A reader—demonstrating its association with poor survival in ccRCC patients (figure 4A). We revealed that ccRCC samples with high expression levels of IGF2BP3 demonstrated greater Th2-cell infiltration enrichment in both the TCGA-KIRC dataset and the immunotherapeutic cohort (figure 4B). Consistent with a previous study,40 Th2 cells were associated with negative outcomes in ccRCC patients (figure 4C). Furthermore, we explored whether the expression of IGF2BP3 and Th2-cell infiltration affected the efficacy of anti-PD-1 therapy. In the anti-PD-1 cohort, a trend in impaired survival was observed in patients with high expression levels of IGF2BP3 (figure 4D). As expected, high Th2-cell infiltration was also associated with poor survival with PD-1 blockade (figure 4E).
Supplemental material
Construction of the m6A gene signature
We applied a methodology (known as m6Ascore) to accurately evaluate the m6A modification pattern in individual ccRCC patients. 299 m6A subtype-related DEGs (figure 5A, (online supplemental table S6) were identified using the DESeq2 package of R software. Univariate Cox regression analysis was performed to evaluate the prognosis of each gene in the m6A subtype-related DEGs; 190 genes conferring significant prognoses were extracted for further PCA to establish the m6A gene signature. Changes in the attributes of individual ccRCC patients were visualized with an alluvial diagram which showed that m6Acluster-C2 had the lowest proportion of the C3 immune subtype and was linked to a high m6Ascore (figure 5B).
The correlation between m6Ascore and the known biological processes was analyzed to better demonstrate the features of the m6A gene signature (figure 5C, online supplemental table S8). It was shown that m6Ascore was negatively correlated with APM (r = –0.22, p<0.001), but positively correlated with mismatch repair-relevant signatures, including mismatch repair, DNA damage repair and DNA replication. Furthermore, the Kruskal-Wallis test showed a significant difference in m6Ascore among m6Aclusters (figure 5D). Next, the prognostic value of the m6Ascore in patients with ccRCC was explored. The patients were divided into high or low m6Ascore groups, with optimal cut-off values determined by the ‘survminer’ R package. Patients with high m6Ascores demonstrated significant survival impairment (figure 5E). We further performed multivariate Cox regression analysis (with factors related to patient sex, age, grade and TNM status included) to investigate the independent prognostic value of m6Ascore, revealing that m6Ascore serves as an independent prognostic biomarker for ccRCC patients (figure 5F). The prognostic value of m6Ascore in ccRCC was validated in another cohort (91 cases) from the ICGC database (figure 5G). The differences in the distribution of somatic mutations between the high and low m6Ascores in the TCGA-KIRC cohort were explored using the maftools package (figure 5H,I). We also found that high m6Ascore was relatively depleted for PBRM1 mutations (26% vs 44%) (p=0.009).
m6A modification patterns in the role of anti-PD-1 therapy
ICTs (eg, anti-PD-1/PD-L1 therapies) have emerged as a critical breakthrough in the field of tumor therapy. We explored the prognostic value of the m6A modification signature in the anti-PD-1 therapy of patients with ccRCC based on one immunotherapeutic cohort. In the anti-PD-1 cohort, the low m6Ascore group presented a markedly prolonged survival (figure 6A,B). In addition, we found that high m6Ascore was relatively depleted for PBRM1 mutations (29% vs 62%) (p<0.001) (figure 6C,D).
Discussion
Increasing evidence reveals that m6A modification plays critical roles in tumorigenesis, therapeutic resistance and immune response via collaboration among m6A regulators.41 Recently, the role of m6A modification in the tumor immune landscape was comprehensively explored in gastric cancer.26 In this study, we focused on the role of m6A modification in the immune landscape in ccRCC to enhance our understanding of the TIME antitumor immune response and provide more effective immunotherapeutic strategies for patients with ccRCC.
Many previous studies have identified ccRCC subtypes based on genomic profiling,42–44 improving the future application of precision-focused, personalized treatments for ccRCC. A 4-mRNA pattern with significant differences in patient survival was identified by unsupervised analyses based on mRNA expression data.42 In the present study—based on 24 m6A regulators—we identified three m6A modification patterns with significantly distinct immune landscapes, characterized by differences in Th17 signature, extent of ITH, overall cell proliferation, aneuploidy, overall somatic copy number alterations, expression of immunomodulatory genes and prognosis. m6Acluster-C3 was the lowest in both proliferation and ITH, suggesting low tumor growth rates in C3. In addition, C3 presented enrichment pathways related to full immune activation and demonstrated the most pronounced Th17 signature, consistent with a previous study demonstrating that Th17 expression is generally related to improved prognosis.45 m6Acluster-C1 also presented enrichment pathways related to full immune activation and relatively high infiltration of CD8+ T cells, while exhibiting high proliferation and ITH, suggesting high tumor growth rates in C1. Consequently, it was not surprizing that C1 exhibited activated immunity but poor survival. m6Acluster-C2 was predominantly associated with immune suppression of biological processes and relatively low infiltration of CD8+ T cells, exhibiting high proliferation and ITH.
In addition, the correlation between each m6A regulator and each TIME infiltration cell type was explored. Our results revealed that high expression levels of IGF2BP3 demonstrated significantly greater enrichment of Th2-cell infiltration. Strikingly, the high infiltration of Th2 cells and expression of IGF2BP3 were both associated with poor survival with PD-1 blockade. It has been reported that Th2-mediated immunosuppression reduced protective cellular immunity and was found to be related to tumor progression.46 Based on the results described above, we speculated that IGF2BP3-mediated m6A modification may promote the infiltration of Th2 cells, thus decreasing the intratumoral antitumor immune response.
Considering the individual heterogeneity of m6A modification, we applied a methodology to accurately evaluate the m6A methylation pattern of individual ccRCC patients known as m6Ascore. Integrated analyses revealed that m6Ascore is a robust and independent prognostic factor for ccRCC. Our study also found that m6Ascore was negatively correlated with APM and high m6Ascore was relatively depleted for PBRM1 mutations. It has been shown that APM was elevated in patients with a partial or complete response to anti-PD-1 therapy but decreased in those with progressive disease on anti-PD-1 therapy.40 PBRM1 mutations were found to be related to improved response, progression-free survival (PFS) and OS with anti-PD-1 therapy in patients with advanced ccRCC.34 Therefore, the above results indirectly revealed that m6A modification may be a critical factor mediating the clinical response to immunotherapy and indirectly confirmed the value of m6Ascore in predicting immunotherapeutic outcomes.
ICTs (such as anti-PD-1/PD-L1 therapies) have revolutionized the treatment of multiple advanced cancers, including ccRCC.12 13 Although significant clinical benefits are achievable in ccRCC patients receiving ICTs, the immunotherapeutic outcomes exhibited individual heterogeneity.14 Therefore, it is of clinical significance to search for markers to predict the outcomes of immunotherapy. The common paradigm in solid tumor immunology is that pre-existing CD8+ T cell infiltration, coupled with a high number of non-synonymous mutations, drives the response to anti-PD-1 therapy.15 16 19 47 48 However, in contrast to other cancer types, neoantigen load, tumor mutation burden and HLA zygosity were not related to anti-PD-1 therapy response in advanced ccRCC.34 Importantly, it was found that there were no statistical differences in response to or survival following anti-PD-1 therapy between immune infiltrated tumors and immune deserts/excluded tumors in advanced ccRCC.34 In this study, we verified the prognostic value of the m6Ascore in the anti-PD-1 therapy of patients with advanced ccRCC. Thus, the m6Ascore may serve as a predictive strategy for anti-PD-1 therapy.
Consequently, we herein provided a new perspective of immuno-oncology and individualized immunotherapy in ccRCC. However, several limitations should be addressed in our study. First, the infiltration of tumor immune cells was obtained based on algorithms owing to technical limitations. Our analyses were also limited by the lack of clinical cohorts to verify the correlation between m6A modification and tumor immune landscape and the prognostic value of m6Ascore in ccRCC. Therefore, further validation based on large-cohort prospective clinical trials are warranted in the future.
In conclusion, this study revealed the correlation between m6A modification and the tumor immune landscape in ccRCC. Our comprehensive evaluation of m6A modification patterns in individual ccRCC patients enhances our understanding of the tumor immune landscape and provides a new approach toward new and improved immunotherapeutic strategies for ccRCC patients.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
JZ and ZL contributed equally.
Contributors GZ and JZ contributed to the study conception; JZ and ZL conducted the data analysis and were responsible for writing the first draft of the paper. CC, XD and TD revised the paper; and all authors read and approved the final version of the manuscript.
Funding This work was financed by grants from the National Natural Science Foundation of China (No.81670643, No.81802821, No.81872437 and No.81870483), the Collaborative Innovation Project of Guangzhou Education Bureau (No.1201620011) the Guangzhou Science Technology and Innovation Commission (No.201704020193) and the Science and Technology Planning Project of Guangdong Province (No.2017B030314108).
Competing interests None declared.
Patient consent for publication Not required.
Provenance and peer review Not commissioned; externally peer reviewed.
Data availability statement All data relevant to the study are included in the article or uploaded as online supplemental information. All data used in this work can be acquired from the GDC portal (https://portal.gdc.cancer.gov/), Broad GDAC Firehose (https://gdac.broadinstitute.org/) and the website (https://gdc.cancer.gov/about-data/publications/panimmune).
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.