m6A related exosome gene siganture in colon cancer identifies distinct microenvironment characterization and immunotherapeutical responses

doi:10.21203/rs.3.rs-895661/v1

Download PDF

Primary research

m⁶A related exosome gene siganture in colon cancer identifies distinct microenvironment characterization and immunotherapeutical responses

https://doi.org/10.21203/rs.3.rs-895661/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background

Gradually, more researches pay attention to the potential relationship between N6-methyladenosine (m⁶A) modification of RNA and immunity in tumor microenvironment (TME). Therefore, they developed unique scoring modes through detecting level of m⁶A in various tumors to learn immune phenotypes. Nevertheless, direct detection always is harmful even fatal to patients. Recent studies demonstrated that exosomes in body fluid would deliver RNA with m⁶A modification from tumors.

Methods

Through analyzing datasets from 1066 colon cancer samples, m⁶A-related exosome genes modification patterns which based on 59 m⁶A-related exosome genes was classified. Then we used MREGS to calculate modification patterns of each tumor with PCA algorithm.

Results

In accordance with 59 m⁶A-related exosome gene, we indicated that there were distinguishing characteristics of TME cell infiltration in three different m⁶A-related exosome gene modification patterns. Three immune cell infiltration profiles (immune-inflamed, immune-desert and immune-excluded) were pronouncedly relevant with the identified patterns in aspects of immune activation, EMT- and TGFβ-related pathways and activity of stroma in TME. We also discovered that tumor molecular subtypes, genetic aberrant expression and patients’ outcomes could be evaluated by identified modification patterns. Furthermore, patients in lower MREGS group, featured as prolonged lifespan and advantageous immune activation, was associated with high tumor mutation burden (TMB) and good prognosis with anti-PDL1 immune therapy.

Conclusion

Detecting the level of m⁶A methylation of exosomes in patients’ serum, MREGS could be utilized for assessing m⁶A-related exosome gene modification patterns of each patient and their own matching characteristics of TME cell infiltration, and next ascertain immunological classifications of tumors and assist clinicians to tailor optimal treatment for patients.

Cancer Biology

Physiology

Biological Chemistry

Colon cancer

m6A modification

Exosome

Tumor microenvironment

Immunotherapy

In the past three decades, exosome undoubtedly has been one of the most pathbreaking scientific achievement in cellular biology. Endosomal-originated vesicles (40–100 nm) exert an influence on cell management and intercellular crosstalk. Exosomes were first discovered in reticulocytes which were relevant with transferrin in 1983, and the size of vesicles was almost 50 nm [1]. Many experimental conclusions over the years show there are tremendous alterations exited in exosomes between normal and cancer cell, showing that both exosome release rate and content (most notably microRNAs) are pronouncedly enhanced. And it is also demonstrated that exosomes from cancer cell lines culture medium are more than normal ones [2].

Hitherto more than 160 RNA chemical modifications have been identified, which exert great influence on RNA structure [3]. Several types of RNA modifications of eukaryotic mRNAs have been widely reported, containing m¹A, m⁶A and m⁵C. And it is well known that m⁶A is the most commonly found in posttranscriptional modification of mRNA and long noncoding RNA in most eukaryotes [4]. Although m⁶A was first identified in the 1970s, there was a large amount of unexploited area in the functions and mechanisms of m⁶A modification. m⁶A methylation is also reversible which is similar to DNA and histone methylation in eukaryotes. Accumulated experimental conclusions elucidated that m⁶A methylation participated in miscellaneous respects of RNA metabolism, containing the processing of mRNA and non-coding RNA (ncRNA), mRNA decay, nuclear export, translation regulation and RNA-protein interactions. m⁶A methylation exert multifarious virtual influence on mammalian development, containing embryogenesis, pluripotency, circadian rhythms, and recently evidence has expounded that tumor genesis and progression are pronouncedly relevant with m⁶A methylation such as the aberrant expression of m⁶A regulators which regulated RNA metabolism in cancer cells. [2, 5].

MicroRNAs (miRNAs) are small non-coding RNAs with 20 to 25 nucleotides. It can be transported into targeted cells and aim to targeted mRNA by exosomes from tumors [6, 7]. Accumulated experimental results proved that the biogenesis of miRNA could be influenced by some situations in cells like cancerization, but the mechanisms are still unclear. And the characteristics of miRNA contained in exosomes hinge on their cell sources. Furthermore, increasing researchers have realized the efficiency and convenience of detecting exosomal content from blood which was called liquid biopsy in some cancer types [8–12]. Combined with tumor specificity of miRNAs, it could serve as diagnostic and prognostic biomarkers in years to come.

Immunotherapy is a heated-discussed topic in oncology. In recent years, the arrival of several innovative monoclonal antibodies used for immunotherapy has revolutionized cancer treatment for a wide range of solid tumors. In Western countries, colon cancer has been a disease with high incidence and high mortality all the way. And about one third of patients with colon cancer develop or develop metastatic disease. Despite many advances in systemic therapies, about 86 percent of patients with advanced stages die within five years of diagnosis. Therefore, immunotherapy has become the most hopeful and effective treatment for tumors at this time [13].

In the past 5 years, a large number of articles have been published on exosomes, but the relationship between m⁶A and exosomes is still unclear. To advance the field of exosome biology, m⁶A-related gene modification and cancer immunotherapy, we have compiled a comprehensive study on the biology and therapeutics of cancer exosomes. Therefore, we developed a scoring system to quantify the m⁶A-related exosome gene modification patterns in individual patients.

Colon Cancer Dataset Source And Preprocessing

Public gene-expression data and full clinical annotation were searched in Gene-Expression Omnibus (GEO) and the Cancer Genome Atlas (TCGA) database. Patients without survival information were removed from further evaluation. In total, 7 eligible CC cohorts (GSE17536, GSE29621, GSE33113, GSE37892, GSE38832 and GSE39582 and TCGA-COAD) were gathered in this study for further analysis. For microarray data from Affymetrix®, we downloaded the raw “CEL” files and adopted a robust multiarray averaging method with the affy and simpleaffy packages to perform background adjustment and quantile normalization. For microarray data from other platforms, the normalized matrix files were directly downloaded. As to datasets in TCGA, RNA sequencing data (FPKM value) of gene expression were downloaded from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) using the R package TCGAbiolinks, which was specifically developed for integrative analysis with GDC data. Then FPKM values were transformed into transcripts per kilobase million (TPM) values. Batch effects from non-biological technical biases were corrected using the “ComBat” algorithm of sva package. The somatic mutation data was acquired from TCGA database. The GSE39582 dataset from GEO was downloaded for Copy Number Variation (CNV) analysis. Data were analyzed with the R (version 3.6.1) and R Bioconductor packages.

Unsupervised clustering for 59 m ⁶ A-telated exosome gene.

Unsupervised clustering analysis was applied to identify distinct m⁶A-related exosome gene modification patterns based on the expression of 59 m⁶A-related exosome gene and classify patients for further analysis. The number of clusters and their stability were determined by the consensus clustering algorithm. We used the ConsensuClusterPlus package to perform the above steps and 1000 times repetitions were conducted for guaranteeing the stability of classification.

Statistical analysis

Correlations coefficients between the TME infiltrating immune cells and expression of m⁶A-related exosome genes were computed by Spearman and distance correlation analyses. One-way ANOVA and Kruskal-Wallis tests were used to conduct difference comparisons of three or more groups. Based on the correlation between MREGS and patients’ survival, the cut-off point of each dataset subgroup was determined using the survminer R package. The “surv-cutpoint” function, which repeatedly tested all potential cut points for finding the maximum rank statistic, was applied to dichotomize MREGS, and then patients were divided into high and low MREGS groups based on the maximally selected log-rank statistics to decrease the batch effect of calculation. The survival curves for the prognostic analysis were generated via the Kaplan-Meier method and log-rank tests were utilized to identify significance of differences. We adopted a univariate Cox regression model to calculate the hazard ratios (HR) for m⁶A-related exosome gene and m⁶A phenotype-related exosome genes. The independent prognostic factors were ascertained through a multivariable Cox regression model. Patients with detailed clinical data were eligible for final multivariate prognostic analysis. The forestplot R package was employed to visualize the results of multivariate prognostic analysis for MREGS in TCGA-COAD cohort. The specificity and sensitivity of MREGS were assessed through receiver operating characteristic (ROC) curve, and the area under the curve (AUC) were quantified using pROC R package. The waterfall function of maftools package was used to present the mutation landscape in patients with high and low MREGS subtype in TCGA-COAD cohort. The R package of RCircos was adopted to plot the copy number variation landscape of 59 m⁶A-related exosome gene in 23 pairs of chromosomes. All statistical P value were two side, with p < 0.05 as statistically significance. All data processing was done in R 3.6.1 software.

Other detailed information for the materials and methods is provided in the additional file 1.

Landscape of genetic variation of m⁶A-related exosome genes in colon cancer

Figure 1a summarized the relationship between m⁶A regulators and exosome-related genes. First of all, a colossal amount of differential expression genes were selected from exosomes in serum between colon cancer patients and normal people, which named as exosome related genes. And then, those associated to m⁶A regulators were sorted out and we obtained 59 m⁶A-related exosome genes. Initially, it is commonly acknowledged that somatic mutation and CNV (copy number variation) are linked with cancers [14, 15]. Therefore, we prefer to acquire a preliminary knowledge of frequency of somatic variation and CNV of 59 m⁶A-related exosome genes in CC patients. In 399 patients, mutations of m⁶A-related exosome genes occurred on 159 patients, accounting for 39.85%. The histogram displayed that RNF43 was ranked first on mutation frequency and ZNF423 subsequently followed. But there are not any mutations exhibiting in TFAM, ALKBH5, C12orf57, GOLGA8A, GOLGA8N, GTF2F2 and NPIPA1 (Fig. 1b). In the following analyses shown the prominent relationship of mutation co-occurrence relationship including ANKRD12 and RNF43, RBM39 and CUL2, besides THOC2 and C2CD5 (Fig. S1a). Through observing distinct CNV alteration of 59 m⁶A-related exosome genes, it is obviously found that copy number amplification occupied with a majority of alterations, whereas CNV deletion happened in LYSMD3, PGGT1B and SRFBP1 (Fig. 1c). Through the analysis of somatic mutation and CNV, it suggested that the chosen 59 genes, which were delivered by exosomes, probably derived from tumors. Subsequently, Fig. 1d exposed that the mutated site of m⁶A-related exosome gene on colon cancer patients’ chromosomes. Then, we enabled to tell CC patients and normal samples apart according to different expression of 59 m⁶A-related exosome genes (Fig. 1e). To assure whether expression of m⁶A-related exosome genes was affected by the aforementioned genetic variations in CC patients, we found that the changes of CNV could exert prodigious influence on perturbations of m⁶A-related exosome gene expression through studying mRNA expression of these 59 genes between normal and CC samples. In addition, m⁶A-related exosome gene with amplificated CNV in CC tissues showed a remarkedly higher expression versus normal colon tissues such as TOMM34 and RBM39 (Fig. 1c and f). No matter in the transcriptome or genomics, there was significant difference between patients and normal people in 59 m⁶A-related exosome genes, suggesting that the expression imbalance of m⁶A-related exosome gene had positive influence on the CC genesis and development.

m⁶A methylation modification patterns mediated by 59 m⁶A-related exosome gene.

We first enrolled six GEO datasets (GSE17536, GSE29621, GSE33113, GSE37892, GSE38832 and GSE39582) which contained clinical data and overall survival (OS) data into one meta-cohort. The comprehensive landscape of m⁶A-related exosome gene interactions for CC patients was depicted with the m⁶A-related exosome gene network (Fig. 2a). The relationship between those exosome related genes were shown in Fig. S1a, and HR value of the genes were shown in Fig. S1b. Then based on the expression of 59 m⁶A-related exosome genes, we were able to divide CC patients into three qualitative m⁶A-related exosome gene modification patterns using the R package of ConsensusClusterPlus. And three patterns (442 cases in pattern 1, 263 cases in pattern 2 and 361 cases in pattern 3) were distinguished by unsupervised clustering algorithm. Finally, these patterns were called as cluster 1-3, respectively (Fig. S2). The tremendous survival preponderance in cluster 1 was disclosed by prognostic analysis for three modification subtypes (Fig. 2b).

TME cell infiltration characteristics in distinct m⁶A-related exosome gene modification patterns.

For the reason that the biological behavior between the three m⁶A modification patterns was not completely understood, we conducted gene set variation analysis (GSVA) of hallmark gene sets in six GEO datasets (GSE17536, GSE29621, GSE33113, GSE37892, GSE38832 and GSE39582). As is shown in Fig. 2c, cluster 3 was markedly enriched in mismatch repair, pyrimidine metabolism and Toll-like receptor signaling pathway (Fig. 2c). However, we observed that innate immune cells such as NK cell, macrophage, MDSC and so on were conspicuously enriched in TME cell infiltration of cluster 2 (Fig. 3a). Previous studies defined that tumors are classified as three phenotypes including immune-excluded, immune-inflamed and immune-desert. And immune-excluded phenotype was characteristic of plentiful immune cells around the tumor cell nest, yet the tumor capsules showed powerful protectivity from penetration of immune cell. And it also demonstrated that poor efficiency on tumor penetration of immune cell may be caused by stromal activation [16]. Hence, we hypothesized that the anti-tumor function of immune cells in cluster 2 was restricted by stromal activation. And it is confirmed by subsequent analyses that the activation of epithelial-mesenchymal transition (EMT), transforming growth factor beta (TGFβ) and angiogenesis pathways, which is relevant to stromal activation in tumor, were remarkably increased in cluster 2 (Fig. 3b). On account of the mentioned analysis, cluster 2 which was featured with innate immune cell infiltration and stromal activation, corresponded with immune-excluded subtype. Cluster 3 which was featured with specific immune cell infiltration and immune activation, corresponded with immune-inflamed subtype; According to the comparison between cluster 1 and 3, it revealed that cluster 1 was not concentrating on antigen processing and presentation, chemokine signaling pathway and cytokine-cytokine receptor interaction which is associated with adaptive immune. For that cluster 1 which was featured with immunological ignorance, corresponded with immune-desert phenotype (Figs. 2c-d and 3a-b). Next, with the method of spearman’s correlation analysis, the mutual effect between each tumor-infiltrating immune cell type and each m⁶A-related exosome gene was demonstrated (Fig. S3).

m⁶A methylation modification patterns in GSE39582 cohort

Putting our eyes on the GSE39582 cohort is beneficial for us to make a comprehensive understanding of m⁶A-related exosome gene modification patterns in numerous clinical cases. Consequently, using unsupervised clustering algorithm, three patterns were distinctly classified in GSE39582 cohort (Fig. 3c and Figs. S4a-d). Also, among three different m⁶A-related exosome gene modification patterns, the prominent difference was shown on the Principal Component Analysis (PCA) scatter diagram. (Fig. 3d). Marisa et al. innovatively classified patients who were suffering from colon cancer into four dominant molecular subtypes including CSC (cancer stem cell), CIN (chromosome instability), KRASm (KRAS mutant) and dMMR (defective mismatch repair). And they concluded that CIN associated with up-regulation of the EMT pathways while dMMR associated with up-regulation of the immune pathways and cell proliferation. However, EMT is downregulated in KRASm subtype [17]. Consistent with the previous findings, CIN subtype patients were mostly divided into cluster 2 and 3, while dMMR subtype predominantly focused in cluster 1 (Fig. 3e). Furthermore, prognostic analysis depicted that cluster 1 had a noteworthy advantage of survival rate comparing to cluster2 and 3 (Fig. S4e). Also, we found that there was different expression of m⁶A regulators among three m⁶A methylation modification patterns (Fig. S4f).

Generation of m⁶A-related exosome gene signatures and functional annotation

With the purpose of discovering unknown biological functioning of m⁶A-related exosome gene modification patterns, we used limma package to identify 3787 m⁶A-related exosome genes phenotype-related DEGs (differential expression genes) (Fig. S4g). And then GO enrichment analysis for differential expression genes was exhibited by the clusterProfiler package. As a consequence, 3787 selected DEGs, which were regarded as m⁶A-related exosome gene signatures, were denoted as the pivotally salient indicator of three m⁶A-related exosome gene modification patterns. Obviously, these signatures performed significant abundance of biological processes that corresponded with the formation of exosome, m⁶A modification and RNA transport, which verified tumors can deliver m⁶A-methylated RNA by exosomes to target cells (Fig. 3f). Owing to the absence of validating to this supervision mechanism, we applied unsupervised clustering analysis which based on the selected 3787 signatures to divide patients into different genomic subtypes and obtained three newly-established distinguishing subtypes called m⁶A-related exosome gene cluster A-C respectively (Figs. S5a-S5d and Fig. 4a). This phenomenon ensured again that there were three distinguishing m⁶A-related exosome gene modification patterns existing in colon cancer. Patients with colon cancer in gene cluster A (213 patients) ranked first in the prognosis analysis. While the worst survival rate was shown in gene cluster B, with 85 patients classified. 259 patients were classified in gene cluster C, which were regarded as an intermediate prognosis (Fig. 4b). Furthermore, among three gene clusters (A-C), it was discovered that the prodigious difference in 24 m⁶A regulators expression was consistent with the anticipated outcome of m⁶A methylation modification patterns (Fig. 4c).

Characteristics of clinical and transcriptome traits in m⁶A-related exosome gene phenotypes.

For the sake of figuring out the function of m⁶A-related exosome gene phenotypes in the TME immunity moderation, we would choose some cytokines and chemokines of three groups, which were gathered from published articles online. First, TGRB1, SMAD9, TWIST1, CLDN3, TGFBR2, ACTA2, COL4A1, ZEB1 and VIM denoted as the transcripts of transforming growth factor (TGFβ)/ EMT pathway. Secondly, PD-L1, CTLA-4, IDO1, LAG3, HAVCR2, PD-1, PD-L2, CD80, CD86, TIGIT and TNFRSF9 were relevant to the transcripts of immune checkpoints. And third, TNF, IFNG, TBX2, GZMB, CD8A, PRF1, GZMA, CXCL9 and CXCL10 were to be correlated with the pathways of immune activation [18, 19]. We considered that gene cluster C was supposed to be stromal-activated group because of its increased mRNA expression of TGFβ/EMT pathway. However, transcriptional mRNAs correlating with immunity system activity were up-regulated in gene cluster A, which suggesting it probably was categorized as the immune-inflamed group (Figs. S5f-5h). The next step is to pick out a portion of known signatures in CC patients on the purpose of describing the function of m⁶A-related exosome signatures (Fig. S5e). As is expected, these box plots showed that gene cluster C was featured with higher stroma activity and cancer progression and metastasis such as increased EMT and WNT target, while cell cycle, DNA replication and mismatch repair were remarkably enhanced in gene cluster A.

In order to more accurately predict the patterns of m⁶A methylation modification in a single tumor, we should make use of an extra method to compensate for the error caused by the heterogeneity of each tumor. Hence, we developed a scoring scheme called the MREGS (m⁶A-related exosome gene score), which is based on the identified 3787 signature genes to calculate the score of individual CC patients. On account of the complicacy of m⁶A-related exosome gene modification, we used an alluvial diagram to observe the flow of data (Fig. 4d). Furthermore, by exploring the relationship between MREGS and some recognized signatures, we could have more knowledge of m⁶A-related exosome gene signatures. And it was found that MREGS was positively associated with EMT-, TGFβ-related pathway and angiogenesis, whereas was negatively related with DNA damage repair and mismatch repair (Fig. 4e). Also, there was prodigious difference exited in MREGS between m⁶A-related exosome gene clusters (cluster 1, 2 and 3) by Kruskal-Wallis test. It displayed that cluster 2 corresponded with high MREGS but cluster 1 represented lower MREGS relatively (Fig. 4f). Additionally, as the activation of fibroblasts are dependent on TGFβ secreted by immunocytes or cancer cells and IL-11 that secreted by TGFβ1-stimulated CAFs (cancer-associated fibroblasts) could enhance the survival rate and invasion ability of cancer cells[20, 21], we chose to analyze stromal activation-related and TGFβ pathways. And we found that increased stromal activation and TGFβ pathway activity in fibroblasts were dominantly corresponded with high MREGS (Fig. 4g). Moreover, dMMR subtype ranked last in MREGS, while CSC subtype was the highest (Fig. 5a). The mentioned description powerfully elucidated that low MREGS was associated with DNA damage repair and high MREGS was associated with stromal activation. Furthermore, we make an attempt to ascertain the significance of MREGS in forecasting patients’ prognosis. Patients with low MREGS demonstrated a prominent survival benefit (Fig. 5b), because their 5-year survival rate was two times higher than patients with high MREGS. Whether MREGS could act as a robust prognostic biomarker for colon cancer is the next question we were supposed to solve. So we made use of Multivariate Cox regression model analysis which contained the element of patients’ gender, age, stage, tumor location and MMR status, subtype and confirmed MREGS as a creditable and independent prognostic biomarker for evaluating patient outcomes (Figs. S6a-b). Due to adjuvant chemotherapy (ADJC) is a common treatment after colon cancer surgery, we designedly validated MREGS to predict the efficacy of ADJC in patients suffering from colon cancer. We discovered that patients with low MREGS-ADJC displayed remarkable adverse reactions compared with patients who did not received adjuvant chemotherapy. While ADJC plays a positive role on patients with high MREGS. The high MREGS-ADJC patients displayed pivotal therapeutic advantages compared with patients who did not received adjuvant chemotherapy. The other obtained results indicated that low MREGS patients invariably displayed the tremendous survival advantage regardless of ADJC treatment (Fig. 5c). Moreover, according to mutation patterns, tumors can be classified into dMMR and pMMR. DMMR (defective mismatch repair) had more mutation burden than pMMR (proficient mismatch repair)[22]. And then we discussed the correlation between MREGS and molecular subtype and discovered that pMMR subtype were distinctly related with higher MREGS. In addition, in the stage Ⅳ patients, there was a remarkable enhancement in MREGS comparing to three other groups (Fig. S6b). This is strongly consistent with the result that dMMR subtype was related to higher survival comparing with pMMR subtype. These phenomena illustrated MREGS could also indirectly assess a portion of clinical features such as MMR subtypes and clinical stage, etc.

Characteristics of m⁶A-related exosome gene modification in TCGA molecular subtypes and tumor somatic mutation

To comprehensively study the characteristics of m⁶A-related exosome gene modification patterns, we introduced TCGA project. Three phenotypes classified by TCGA project is comprised of chromosomal instability (CIN), invasive and microsatellite instability (MSI). Next, the difference of MREGS among three phenotypes was calculated. The higher MREGS was pronouncedly concentrated on CIN and had a shorter lifespan, whereas invasive phenotype is relevant with lower MREGS, which was related to better survival (Figs. 5d-e). Different stages indicated differential expression of m⁶A-related exosome genes in TCGA project（Fig. 5f）. Additionally, we used the maftools software package to analyze the differences in the distribution of high MREGS and low MREGS somatic mutations in the TCGA-COAD cohort. Low MREGS group presented more extensive tumor mutation burden than the high MREGS group (Figs. 5g-h). Accumulated evidence demonstrated a potential connection between enhanced survival rate whose receiving PD-1/PD-L1 immunological therapy and higher somatic tumor mutation burden (TMB) was discovered. Consequently, it was indirectly elucidated that clinical reactions to some immunological treatment like anti-PD-1/PD-L1 drugs may depend on different m⁶A-related exosome gene modification patterns in tumors. Also, we confirmed MREGS as a reliable method to predict prognostic outcomes after immunotherapy. In both clinical trials and preclinical studies, higher TMB patients with receiving immune checkpoint inhibitor treatment has a prominent superiority in survival rate and clinical response [23].

m⁶A-related exosome gene modification patterns in the role of anti-PD-1/L1 immunotherapy

With the purpose of testifying the reliability of MREGS and its prognostic value, MREGS signatures, obtained from GSE39582 cohort, was applied to five other colon cancer datasets (GSE17536, GSE29621, GSE33113, GSE37892, GSE38832; Figs. S8c-g). The ROC curve was used to evaluate the accuracy of a model, and higher the AUC value, the more accurate it is. And MREGS model displayed high accuracy of predictive superiority in patients with 3-year (AUC=0.785) and 5-year (AUC=0.754) colon cancer (Figs. S8h-i). Immunological therapy exemplified by anti-PD-L1/ PD-1 drugs have strikingly made immense progress in molecular targeted cancer therapy in recent years. Therefore, we selected two immunotherapy datasets (IMvigor210, GSE78220) to validate the predictive ability of MREGS to those patients who had been treated PD-L1/ PD-1 inhibitors. In anti-PD-L1 dataset (IMvigor210), through analyzing clinical response and survival rate, we exposed that low MREGS group surpassed high MREGS group (Figs. 6a-c). While in anti-PD-1 dataset (GSE78220), there was no significant difference between patients with low MREGS and high MREGS (Figs. 6d-f). The following analysis unraveled TME stroma and TGFβ pathway in fibroblasts were notably activated in high MREGS group, which mediated tumor migration (Fig. 6h). Tumor neoantigen burden (TNB) is associated with efficacy of immunological therapy and a major element in the judgement of clinical immunotherapies. And high TNB would be expected to be characterized by extensive T cell responses and to be specifically sensitive to immunotherapy[24]. And then we discovered that low MREGS group with high neoantigens displayed significant preponderance in survival rate (Fig. 6i). As aforementioned, MREGS, quantification of m⁶A-related exosome gene signatures, is anticipated to be an underlying and reliable biomarker for evaluating patients’ outcomes and clinical manifestation after treating immunological therapy (Fig. 6j). All in all, calculated from m⁶A-related exosome gene signatures, MREGS is available for forecasting patients’ clinical response whose receiving anti-PDL1 drugs and being a referable indicator for judgement of surgeons.

Increasing evidence illustrated that m⁶A modification acted as a monumental role in inflammation, innate immunity, and antineoplastic interaction with many m⁶A regulators. Because of a majority of researches individually paid attention to a kind of immunocyte in tumor environment or a functioning regulator, it was unknown about the characteristics of TME infiltration mediated by combined function of multiple m⁶A regulators. Exosome, a vector containing RNA and protein to promote tumors growth, can travel through the whole body. Also, it can deliver “information” from tumors, which can be learned more conveniently than needling biopsy [6, 25]. Thus, through decoding “information” from exosomes in the body fluid, ascertaining the status of distinguishing m⁶A-related exosome gene modification patterns will strengthen the cognitive knowledge of TME immunological reaction and enlighten more targeted immune therapy like personalized treatment.

In this study, in accordance with 59 m⁶A-related exosome gene, we indicated that there were distinguishing characteristics of TME cell infiltration in different m⁶A-related exosome gene modification patterns. Cluster 1, featured by the suppression of immunity, was divided into immune-desert group; Cluster 2, enriched by natural immunocytes and enhanced activity of stroma, was divided into immune-excluded group; Cluster 3, which was fulfilled of acquired immune cell in tumors, was divided into immune-inflamed group. As previous researches demonstrated, the scientists would name the tumors in immune-excluded and immune-desert group as non-inflamed tumors. On the other hand, the tumors in immune-inflamed group manifested mountainous immune cell infiltration in TME [16, 26, 27]. It is shown that the richness of immunocytes was high in immune-excluded phenotype. However, instead of penetrating the tumor parenchyma, immunocytes just were inclined to retain in the stroma which was encircled by tumor cells [28–30]. On the basis of newly reports, increased TGFβ and activation of EMT-related pathways in the tumor microenvironment would contribute to the mechanism of immune evasion that block the way to immune infiltration [31, 32]. Some specific molecular inhibitors aiming to TGFβ, which enable to renovate the tumor microenvironment and liberate T cells from tumors[32, 33]. Moreover, accordant with the above conclusions, cluster 2 showed a remarkable status of stromal activation, containing high expression of EMT and TGFβ pathways and angiogenesis factors, which were considered as suppression of T-cell activation. Connected with each cluster’s characteristics of cell infiltration in tumor microenvironment, it is ensured that the immune-phenotype classification for three modification patterns is considerably credible. In accordance with the conclusion we made, it is no question that a series of activated innate immunocytes existed in the cluster 2 but it displayed poorer survival rate.

Next, mRNA transcriptome differences between distinguishing m⁶A-related exosome gene modification patterns were testified to be remarkably relevant with the biology pathways of m⁶A modification, formation of exosomes and RNA transport. And these differentially expressed genes were named as m⁶A-related exosome gene signatures. Consistent with the clustering results of the m⁶A modification patterns (Cluster 1, 2, 3), three gene clusters (A, B, C) which were clustered by m⁶A-related exosome gene signatures, were also significantly correlated with stromal activation and immune response. Hence it will increase our knowledge of characteristic of TME cell-infiltrating through assessing of m⁶A-related exosome gene modification patterns comprehensively. Owing to difference of m⁶A modification in different tumors, every tumor is supposed to be quantified its personal m⁶A modification patterns. Therefore, we set a scoring role to individually quantify the level of m⁶A methylation in exosomes of sufferers with colon cancer–m⁶A-related exosome gene score (MREGS). As mentioned results showed, the immune-excluded-related modification pattern displayed a higher MREGS, whereas lower MREGS was exhibited in the immune-inflamed-related pattern. It indicated that MREGS is one of reliable indexes to evaluate m⁶A-related exosome gene modification pattern of individual tumor and further identify tumor immune phenotype.

Our data also illustrated that low MREGS showed a strong correlation with high tumor mutation burden (hTMB) which had developed a more powerful and sensitive biomarker to immune checkpoint inhibitors and a positive index to receive immunotherapeutic treatments like immune checkpoint inhibitor (ICPI) [34]. In addition, low MREGS group was related to dMMR (different mismatch repair) subtype. As previous reports said, metastatic colorectal cancer patients with dMMR could have tremendous benefit from immune checkpoint inhibitor such as PD-L1 monoclonal antibody pembrolizumab [35]. Moreover, we found that low MREGS group with high neoantigens had higher survival rate than high MREGS group with high neoantigens. Anagnostou et.al revealed that along with the more identification of neoantigens, T cells fighting against the tumors was increased by immune system thereby raising the efficacy of ICPI [9]. The MREGS with integrated some kinds of biomarkers like mutation load, neoantigen load, MMR status and stromal activation and TME immune phenotypes, could be more beneficial for designing strategies for immunotherapy. Also, in the cohort with anti-PD-L1 immune therapy, the predictive value of MREGS was confirmed and remarkable difference existed in MREGS between responders and non-responders. Interestingly, through analyzing the correlation between adjuvant chemotherapy (ADJC) and MREGS, high MREGS group was more suitable for ADJC than low MREGS group because ADJC took a negative effect on low MREGS group reflecting in survival rate. Combined with two fact that high MREGS was relevant with the augmented expression of TGFβ and EMT and higher possibility of tumor recurrence happened in patients with increased TGFβ- and EMT-related pathways [36], it’s not hard to draw a conclusion that ADJC can restrain the metastasis of residual cancer cells. On the other hand, ADJC also had a certain toxicity on patients as the lower survival rate in low MREGS with ADJC group versus with non-ADJC group. Consequently, MREGS is probably expected to be a biomarker to sift optimal sufferers that need to receive ADJC like microsatellite status and BRAF and KRAS mutations [37].

In sum, detecting the level of m⁶A methylation of exosomes in patients’ serum, MREGS could be utilized for assessing m⁶A-related exosome gene modification patterns of each patient and their own matching characteristics of TME cell infiltration, and next ascertain immunological classifications of tumors and assist clinicians to tailor optimal treatment for patients. We also expounded MREGS could evaluate patients’ clinicopathological features such as genetic variation, status of tumor cell infiltration, MMR status, clinical stages and TMB etc. We could find the correlation between MREGS and clinicopathological characteristics in our analysis graphics. Furthermore, predicting the clinical reaction to anti-PDL1 immune therapy and the therapeutic effect of ADJC will be one of MREGS’s greatest advantages. Although this approach had the advantages of good as convenient detection, it also had a number of unresolved challenges, for example extraction and purity of exosomes, sensitivity of detected m⁶A methylation and how to confront cancer metastasis which meant lesions were not just found in one tissue of one patient. From our perspective, the current solutions are continuously reforming exosome extraction technologies and platforms and deeply mining clinical data to explore tumor-specific exosome markers.

CC: Colon cancer; CNV: Copy number variation; DEG: Differential expression gene; TME: Tumor microenvironment; GSVA: Gene set variation analysis; EMT: Epithelial-mesenchymal transition; TGFβ: Transforming growth factor-β; PCA: Principal component analysis; CSC: Cancer stem cell; CIN: Chromosome instability; dMMR: Defective mismatch repair; pMMR: Proficient mismatch repair; MREGS: m⁶A-related exosome gene score; ADJC: Adjuvant chemotherapy; MSI: Microsatellite instability; CIMP: CpG island methylator phenotype; TMB: Tumor mutation burden.

Acknowledgements

We would like to thank the Shandong University for supporting this manuscript.

Funding

This work was supported by grants from the National Natural Science Foundation of China (31971061).

Authors' contributions

J.Z. analyzed the data and wrote the manuscript; C.Z., Y.J. and Y.L. analyzed the data; D.C. was responsible for interpretation of data and writing of the manuscript; F.X, J.L. conceptualized the project and funded the work.

Ethics approval and consent to participate

Our study is based on open source data, so there are no ethical issues and other conflicts of interest.

Consent for publication

All the co-authors agreed to publish the final version of this manuscript.

Competing interests

The authors declare no potential conflicts of interest.

Pan B, Johnstone R. Fate of the transferrin receptor during maturation of sheep reticulocytes in vitro: selective externalization of the receptor. Cell. 1983;33(3):967–78.
He L, Li H, Wu A, Peng Y, Shu G, Yin G. Functions of N6-methyladenosine and its role in cancer. Molecular Cancer 2019, 18(1).
Boccaletto P, Machnicka MA, Purta E, Piatkowski P, Baginski B, Wirecki TK, de Crecy-Lagard V, Ross R, Limbach PA, Kotter A, et al. MODOMICS: a database of RNA modification pathways. 2017 update. Nucleic Acids Res. 2018;46(D1):D303–7.
Liu ZX, Li LM, Sun HL, Liu SM. Link Between m6A Modification and Cancers. Front Bioeng Biotechnol. 2018;6:89.
Sun T, Wu R, Ming L. The role of m6A RNA methylation in cancer. Biomed Pharmacother. 2019;112:108613.
Milane L, Singh A, Mattheolabakis G, Suresh M, Amiji MM. Exosome mediated communication within the tumor microenvironment. J Control Release. 2015;219:278–94.
Falcone G, Felsani A, D'Agnano I. Signaling by exosomal microRNAs in cancer. J Exp Clin Cancer Res. 2015;34:32.
Li S, Yi M, Dong B, Tan X, Luo S, Wu K. The role of exosomes in liquid biopsy for cancer diagnosis and prognosis prediction. Int J Cancer. 2021;148(11):2640–51.
Anagnostou V, Smith KN, Forde PM, Niknafs N, Bhattacharya R, White J, Zhang T, Adleff V, Phallen J, Wali N, et al. Evolution of Neoantigen Landscape during Immune Checkpoint Blockade in Non-Small Cell Lung Cancer. Cancer Discov. 2017;7(3):264–76.
Ge Y, Mu W, Ba Q, Li J, Jiang Y, Xia Q, Wang H. Hepatocellular carcinoma-derived exosomes in organotropic metastasis, recurrence and early diagnosis application. Cancer Lett. 2020;477:41–8.
Yang KS, Ciprani D, O'Shea A, Liss AS, Yang R, Fletcher-Mercaldo S, Mino-Kenudson M, Fernandez-Del Castillo C, Weissleder R. Extracellular Vesicle Analysis Allows for Identification of Invasive IPMN. Gastroenterology. 2021;160(4):1345–58 e1311.
Guo D, Xu Y, Ding J, Dong J, Jia N, Li Y, Zhang M. Roles and Clinical Applications of Exosomes in Cardiovascular Disease. Biomed Res Int. 2020;2020:5424281.
Wrobel P, Ahmed S. Current status of immunotherapy in metastatic colorectal cancer. Int J Colorectal Dis. 2019;34(1):13–25.
Shao X, Lv N, Liao J, Long J, Xue R, Ai N, Xu D, Fan X. Copy number variation is highly correlated with differential gene expression: a pan-cancer study. BMC Med Genet. 2019;20(1):175.
Martincorena I, Campbell PJ. Somatic mutation in cancer and normal cells. Science. 2015;349(6255):1483–9.
Chen DS, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. 2017;541(7637):321–30.
Marisa L, de Reynies A, Duval A, Selves J, Gaub MP, Vescovo L, Etienne-Grimaldi MC, Schiappa R, Guenot D, Ayadi M, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. PLoS Med. 2013;10(5):e1001453.
Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, Schinzel AC, Sandy P, Meylan E, Scholl C, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. 2009;462(7269):108–12.
Zeng D, Li M, Zhou R, Zhang J, Sun H, Shi M, Bin J, Liao Y, Rao J, Liao W. Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol Res. 2019;7(5):737–50.
Kalluri R. The biology and function of fibroblasts in cancer. Nat Rev Cancer. 2016;16(9):582–98.
Calon A, Espinet E, Palomo-Ponce S, Tauriello DV, Iglesias M, Cespedes MV, Sevillano M, Nadal C, Jung P, Zhang XH, et al. Dependency of colorectal cancer on a TGF-beta-driven program in stromal cells for metastasis initiation. Cancer Cell. 2012;22(5):571–84.
Ganesh K, Stadler ZK, Cercek A, Mendelsohn RB, Shia J, Segal NH, Diaz LA. Jr.: Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol. 2019;16(6):361–75.
Samstein RM, Lee CH, Shoushtari AN, Hellmann MD, Shen R, Janjigian YY, Barron DA, Zehir A, Jordan EJ, Omuro A, et al. Tumor mutational load predicts survival after immunotherapy across multiple cancer types. Nat Genet. 2019;51(2):202–6.
Schumacher TN, Schreiber RD. Neoantigens in cancer immunotherapy. Science. 2015;348(6230):69–74.
Ruivo CF, Adem B, Silva M, Melo SA. The Biology of Cancer Exosomes: Insights and New Perspectives. Can Res. 2017;77(23):6480–8.
Turley SJ, Cremasco V, Astarita JL. Immunological hallmarks of stromal cells in the tumour microenvironment. Nat Rev Immunol. 2015;15(11):669–82.
Gajewski TF, Woo SR, Zha Y, Spaapen R, Zheng Y, Corrales L, Spranger S. Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment. Curr Opin Immunol. 2013;25(2):268–76.
Gajewski TF. The Next Hurdle in Cancer Immunotherapy: Overcoming the Non-T-Cell-Inflamed Tumor Microenvironment. Semin Oncol. 2015;42(4):663–71.
Joyce JA, Fearon DT. T cell exclusion, immune privilege, and the tumor microenvironment. Science. 2015;348(6230):74–80.
Salmon H, Franciszkiewicz K, Damotte D, Dieu-Nosjean MC, Validire P, Trautmann A, Mami-Chouaib F, Donnadieu E. Matrix architecture defines the preferential localization and migration of T cells into the stroma of human lung tumors. J Clin Invest. 2012;122(3):899–910.
Tauriello DVF, Palomo-Ponce S, Stork D, Berenguer-Llergo A, Badia-Ramentol J, Iglesias M, Sevillano M, Ibiza S, Canellas A, Hernando-Momblona X, et al. TGFbeta drives immune evasion in genetically reconstituted colon cancer metastasis. Nature. 2018;554(7693):538–43.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554(7693):544–8.
Panagi M, Voutouri C, Mpekris F, Papageorgis P, Martin MR, Martin JD, Demetriou P, Pierides C, Polydorou C, Stylianou A, et al. TGF-beta inhibition combined with cytotoxic nanomedicine normalizes triple negative breast cancer microenvironment towards anti-tumor immunity. Theranostics. 2020;10(4):1910–22.
Chalmers ZR, Connelly CF, Fabrizio D, Gay L, Ali SM, Ennis R, Schrock A, Campbell B, Shlien A, Chmielecki J, et al. Analysis of 100,000 human cancer genomes reveals the landscape of tumor mutational burden. Genome Med. 2017;9(1):34.
Le DT, Uram JN, Wang H, Bartlett BR, Kemberling H, Eyring AD, Skora AD, Luber BS, Azad NS, Laheru D, et al. PD-1 Blockade in Tumors with Mismatch-Repair Deficiency. N Engl J Med. 2015;372(26):2509–20.
Wang Y, Shi J, Chai K, Ying X, P Zhou B. The role of snail in EMT and tumorigenesis. Curr Cancer Drug Targets. 2013;13(9):963–72.
Auclin E, Zaanan A, Vernerey D, Douard R, Gallois C, Laurent-Puig P, Bonnetain F, Taieb J. Subgroups and prognostication in stage III colon cancer: future perspectives for adjuvant therapy. Ann Oncol. 2017;28(5):958–68.

Supplementarymaterial.docx

Download PDF

Editor assigned by journal
15 Sep, 2021
Submission checks completed at journal
14 Sep, 2021
Editor invited by journal
14 Sep, 2021
First submitted to journal
13 Sep, 2021

You are reading this latest preprint version

m⁶A related exosome gene siganture in colon cancer identifies distinct microenvironment characterization and immunotherapeutical responses

Status:

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

Introduction

Methods

Colon Cancer Dataset Source And Preprocessing

Statistical analysis

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1