Landscape Of M6a RNA Methylation Regulators and Tumor Microenvironment Cell-Inltration Characterization In Gastroesophageal Adenocarcinomas

Background: RNA N6-methyladenosine (m 6 A) modication plays a nonnegligible role in shaping individual tumor microenvironment (TME) characterizations. However, the landscape and relationship of m 6 A modication and TME cell inltration remain unknown in gastroesophageal adenocarcinomas (GEA). Methods: We systematically examined the TME of GEA focusing on the distinct m 6 A modication patterns from the public databases. Intrinsic patterns of m 6 A modication were evaluated for associations with clinicopathological characteristics, underlying biological pathways, tumor immune cell inltration, oncological outcomes and treatment responses. We generated a single-cell transcriptome atlas of the GEA sample inhouse to validate the role of the m 6 A modication pattern on TME. Results: We identied and validated the landscape of m 6 A regulators and tumor-inltrating immune cells in GEA. Then, two distinct m 6 A modication patterns of GEA (cluster1/2 subgroup) were dened, and we correlated two subgroups with TME cell-inltrating characteristics. Cluster2 subgroup correlates with a poorer prognosis, down-regulated PD-1 expression, higher risk scores and distinct immune cell inltration. Additionally, PPI and WGCNA network analysis were integrated to identify key module genes closely related to immune inltration of GEA to nd immunotherapy markers. And COL4A1 and COL5A2 in brown module were signicantly correlated to the prognosis, PD-1/L1 and CTLA-4 expression of GEA patients. Interesting, low COL5A2 expression was linked to an enhanced response to anti-PD-1 immunotherapy. Finally, a prognostic risk score was constructed using three m 6 A regulator-associated signatures that represented an independent prognosis factor for GEA. The single-cell transcriptome atlas of GEA sample validated the role of m 6 A modication pattern on TME and revealed that three m 6 A regulators are highly expressed in CD4 + T cells, CD8 + T cells, Tregs and Macrophages. Conclusions: Our work revealed m 6 A RNA methylation regulators are a type of


Background
Gastroesophageal adenocarcinomas (GEA) is still a major cause of cancer-related mortality worldwide [1].Currently, the development of effective targeted therapeutics for GEA patients lags behind other cancers.Despite recent improvements in the multidisciplinary and multimodality treatment, in fact, the overall prognosis for GEA patients remains poor, with a global 5-year survival rate lower than 30% for gastric cancer (GC) and about 19% for esophagus adenocarcinoma [2].For purpose of the high heterogeneity and complicated disease processes of GEA, there is still a lack of effective prognostic markers in this disease.Therefore, identifying molecular biomarkers and novel potential therapies is critical to predict GEA patients' prognosis and determine the personalized treatment.
Notably, N6-methyladenosine (m 6 A), the most abundant modi cation on mRNAs in eukaryotes, is closely related to stem cell differentiation, immune response, embryonic development, and microRNA (miRNA) editing; they also play an essential role in the progression of various cancers [3][4][5][6][7].The m 6 A methylation levels in tumors mainly depend on the expression of m 6 A methylation regulators.m 6 A is modulated by methyltransferase complex ("writers"), demethylases ("erasers") and RNA-binding proteins("readers"), which performs a series of biological functions [8].The aberrant expression of m 6 A regulators plays a vital regulatory role in tumor progression, prognosis, and radio-resistance.Yongsheng Li et al. [9] have showed the characteristics of m 6 A RNA methylation across 33 types of cancer and speculated that the mechanism of m 6 A RNA modi cation might be associated with the activation or depression of some oncogenic pathways such as PI3K-AKT-mTOR signaling, KRAS and P53 pathways.However, given the limited knowledge of the role of m 6 A methylation in GEA, studying the precise correlation between m 6 Arelated regulator genes and its clinical prognosis is in high demand.
Immunotherapy represented by immunological checkpoint blockade (ICB, PD-1/L1and CTLA-4) has demonstrated surprising clinical e cacy in a small number of patients with durable responses.In September 2017, the U.S. Food and Drug Administration (FDA) granted accelerated approval for pembrolizumab for the treatment of patients with recurrent, locally advanced or metastatic, gastric adenocarcinoma or GEA, whose tumors expressed PD-L1 and with disease progression on or after 2 or more systemic therapies.Disappointingly, response rates of immune checkpoint inhibitor monotherapy in GEA are around 10%-25% depending on the number of previous lines of chemotherapy and PD-L1 status [10].Hence, it is important and necessary to understand the complexity of the tumor microenvironment (TME) and identify subclasses of tumor immune microenvironment existing in the patients' tumor to predict and administer corresponding immunotherapy.Notably, several studies have indicated the special relationship between TME in ltrating immune cells and m 6 A modi cation.For instance, Dali et al. [11] showed that loss of YTHDF1 in classical dendritic cells enhanced the cross-presentation of tumor antigens and the cross-priming of CD8 + T cells in vivo and YTHDF1 may be a potential therapeutic target in anticancer immunotherapy.Bo Zhang et al. [12] determined three distinct m 6 A modi cation patterns in gastric cancer and found that the TME cell-in ltrating characteristics under these three patterns were highly consistent with the three immune phenotypes of tumors.Seungwon Yang et al. [13] suggested that m 6 A demethylation by fat mass and obesity-associated protein (FTO) increases melanoma growth and decreases response to anti-PD-1 blockade immunotherapy.Na Li et al. [14] showed that Alkbh5 regulated the composition of tumor-in ltrating Treg and myeloid-derived suppressor cells and sensitized tumors to cancer immunotherapy.However, until now, the role of m 6 A regulators in the malignancy and prognosis of GEA has not been comprehensively clari ed.Therefore, research focusing on m 6 A regulators is warranted to elucidate the potential regulatory mechanism of m 6 A methylation in TME, which may reveal the potential mechanism and targets of immunotherapy.
In this study, we systematically evaluated the role of m 6 A modi cation, and correlated the m 6 A modi cation with the TME cell-in ltrating characteristics in GEA.Two GEA subtypes (cluster1/2) were determined via the consensus clustering for m 6 A regulators that strati ed the prognosis of patients, different TIICs and PD-1 expression.After WGCNA analysis, low COL5A2 expression was found to be linked to enhanced response to anti-PD-1 immunotherapy.Risk score developed from three m 6 A regulatorbased signatures was an independent prognostic indicator of patients with GEA.The m 6 A regulatorbased risk signatures were signi cantly related to the immune cell in ltration levels of patients with GEA.
We employ the single-cell RNA sequencing (scRNA-seq) to validate m 6 A regulatory genes in distinct cell populations of TME in GEA.Therefore, this study sought to provide insights into the regulatory mechanisms associated with TME and the strategies for GEA immunotherapy.

Data processing
The mRNA (RNA sequencing) Fragments Per Kilobase of transcript per Million Fragments standardized expression data and corresponding clinicopathological features of TCGA-STAD&ESCA cohorts were retrieved for 159 GEA tissues and 39 adjacent non-tumor tissues from The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/) and 121 GEA tissues from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/).Patients without prognostic information were excluded from the analysis.The dataset of GSE96669 was obtained using the GPL10558 platform (Illumina HumanHT-12 V4.0 expression beadchip).We utilized the limma package to conduct the normalization process, deleting the normal or repeated samples for subsequent analysis.Then, the clinicopathological parameters for included samples also were download from the TCGA database.The relevant data TCGA and GEO provided is publicly available and open source; hence, approval by a local ethics committee was not required.

Evaluation of tumor-in ltrating immune cells (TIICs)
CIBERSORT algorithm was applied to calculate the fractions of the 22 types of TIICs [15], which is considered better than previous deconvolution methods for the analysis of unknown mixture content and noise.We used this algorithm to statistically estimate the relative proportions of cell subpopulations from complex tissue expression pro les, making it a useful tool to estimate the abundances of special cells in mixed tissue.In this research, we used the R package "CIBERSORT" to estimate the fraction of immune cells of TCGA samples, which followed by quality ltering that tumor samples with P < 0.05 were selected for the following analysis.
Generation of immune score, stromal score, and ESTIMATE score ESTIMATE algorithm was exploited to infer the fraction of immune and stromal cells in tumor tissues based on gene expression signature, including microarray expression data sets, new microarray, as well as RNA-seq transcriptome pro les.The R script of the ESTIMATE algorithm was downloaded from the public source website (https://sourceforge.net/projects/estimateproject/).Then, we calculated the immune scores, stromal scores and ESTIMATE scores for each sample of the TCGA dataset, respectively.And the higher the respective score, the larger the ratio of the corresponding component in TME.After we got three scores from the ESTIMATE method, we could classify the samples into high-and low-level groups according to the median score, respectively.

Unsupervised clustering of m 6 A methylation regulators
In order to further investigate the function of m 6 A RNA methylation regulators in GEA, we clustered the GEA patients into different groups by using the R package ConsensusClusterPlus (50 iterations, resample rate of 80%, and Pearson correlation, http://www.bioconductor.org/) based on the expression of the 21 m 6 A RNA methylation regulators [16].The number of clusters and their stability were determined by the consensus clustering algorithm.Principal Components Analysis (PCA) was used with the R package for R v3.6.3 to study the gene expression patterns in different GEA groups.

Differentially expressed genes (DEGs)
We used R package "limma" with log 2 | fold-change (FC) | > 1 and adjusted P-value < 0.05 to perform differentiation analysis of the gene expression, and DEGs were generated by the comparison between GEA samples vs. adjacent non-cancerous samples in TCGA and GSE96669 datasets.Venn online software (http://bioinformatics.psb.ugent.be/webtools/Venn/)was used to identify the overlapping DEGs between tumor and normal samples.

Weighted gene co-expression network analysis (WGCNA) of DEGs
WGCNA is a useful tool to establish the co-expression network between gene pattern and clinical traits using WGCNA package in R based on the RNA-seq data from TCGA database [17].In the rst step, we calculated a similarity matrix using bi-weight mid-correlation, as it is more robust to outliers.After that, a weighted adjacency matrix was de ned by raising the co-expression similarity to appropriate soft-thresholding power.The best power (β-value) was chosen based on the criterion of approximate scalefree topology.Next, we transformed the adjacency into a topological overlap matrix (TOM) and calculated the corresponding dissimilarity to minimize the effects of noise and spurious associations.Hierarchical clustering was used to produce a hierarchical clustering tree and dynamic tree cut method to assign coexpressed genes to each module.Modules were constructed with a minimum module size of 20 genes, and highly similar modules were combined using a dissimilarity threshold of 0.25.

Screening signi cant modules and functional enrichment analysis
In order to identify the signi cance of each module, gene signi cance (GS) was calculated using linear regression by log10 conversion of the p value between gene expression and clinical features.And module eigengenes (MEs) were de ned as the rst principal component of each gene module and adopted as the representative of all genes in each module.Next, we calculated the correlation between gene modules and clinical traits by WGCNA package in R and draw a heatmap.After obtaining these, genes in the gene modules from Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed to observe the function of selected signi cant gene modules using the clusterPro ler package in R. Enriched terms and pathways with adjusting P value < 0.05 were selected.

Protein-protein interaction (PPI) network and Hub genes identi cation
For clarifying the drivers of inducing carcinogenesis in a more reliable way, PPI analysis was performed necessarily.The Retrieval of Interacting Genes (STRING) database (http://string-db.org)online tool was used to evaluate interactive relationships and generate PPI networks among the DEGs in selected gene modules.Interaction score 0.7 served as cutoff value prior to visualization.Then Cytoscape software (http://cytoscape.org/development_team.html)was selected to visualize the results of PPI networks.Furthermore, CytoHubba app identifying hub objects from complex interaction in Cytoscape software was used to nd top hub genes.Subsequently, top hub genes were selected and ranked by the Maximal Clique Centrality (MCC) method.Afterward, to select key genes that affect prognosis, survival data including living status and survival time was extracted from the TCGA database.Kaplan-Meier survival curves were built to screen for genes signi cantly associated with prognosis.

Construction of m 6 A related genes signature
Univariate Cox regression analysis of the expression of 21 m 6 A RNA methylation regulators was conducted to determine the candidate genes associated with overall survival (OS).After that, regulators associated with OS in univariate analyses were subsequently selected for the least absolute shrinkage and selection operator (LASSO) Cox regression to construct a m 6 A-related risk signature for clinical prognosis [18].Finally, three m 6 A RNA methylation regulators with their corresponding coe cients were determined by the minimum mean cross-validated error, choosing the optimal penalty parameter λ related to the minimum 10-fold cross validation within the training set.The risk score of each patient with GEA in the TCGC cohort was calculated using the following formula: Where Xi is the standardized expression value of each selected m 6 A RNA methylation regulator, and Coe is the corresponding coe cient of the gene.All patients were divided into low-and high-risk groups based on the median value of the risk scores.Survival curves in the high-risk and low-risk groups were estimated using the Kaplan-Meier method.Additionally, the receiver operating characteristic (ROC) curves and area under the ROC curves (AUC values) were applied to access sensitivity and speci city.And AUC > 0.5 was considered as a signi cant diagnostic model.

Gene set enrichment analysis (GSEA)
GSEA is a computational method usually used to determine whether a set of basically de ned gene sets exhibits statistically signi cant differences between two biological states.GSEA was provided by the JAVA program with MSigDB v7.1 and downloaded from the website of Broad Institute [19].According to the median value of RS, the samples were divided into two groups, and "c2.cp.kegg.v7.1.symbols.gmt"gene set enrichment analysis was carried out, with p-value < 0.05 and q-value < 0.05 as indicative of statistical signi cance.The enrichment pathway was visualized using the R packages "ggplot2" and "clusterPro ler".

Patients and sample information
We totally collected 16 non-neoplastic and neoplastic samples from GEA patients who underwent surgical treatments in the Gastrointestinal Surgery Department of Jinan Central Hospital A liated to Shandong University from 2018 to 2020.Fresh samples from two patients were immediately used for single-cell RNA-seq analysis, and the rest of fresh tumor and non-neoplastic tissues were frozen and stored at -80℃ that used for PCR analysis.Clinical characteristics of included patients was shown in supplementary table S1.This research was approved by the Medical Ethics Committee of Jinan Central Hospital A liated to Shandong University and the sample acquisition and usage was performed in accordance with the approved guidelines.Informed consent was acquired from each involved patient.

Quantitative real-time polymerase chain reaction (qRT-PCR)
For evaluating the expression levels of three signature regulators and hub genes, we extracted total RNA from clinical GEA samples by using RNA trizol reagent (CWBIO).According to the instructions of the manufacturer, cDNA synthesis was carrying out by using reverse transcription kit (CWBIO).The qRT-PCR analysis was conducted on The LightCycler 480 Real-Time PCR System.Related gene expression levels were calculated using the 2 −△△CT method and the related GAPDH mRNA expression was used as an endogenous control.Primer sequences are presented in supplementary Table S2.

Tissue dissociation and preparation of single-cell suspensions
Two tumors and One adjacent normal mucosa were processed immediately after being obtained from GEA patients.Fresh samples were then washed with phosphate-buffered saline (PBS, CORNING).Tissues were dissociated into single cells in dissociation solution (Tumor Dissociation Kit, human; Miltenyi Biotec) in 37 ℃ water bath with shaking for 20 min at 100 rpm.Digestion was terminated with 1× PBS containing 10% FBS(CORNING).Cell suspensions were ltered using a 70 µm lter and then centrifuged at 300g for 5 min at 4°C to pellet dead cells and red blood cells.The cells were washed twice and resuspended in PBS with 0.04% bovine serum albumin (BSA, Sigma).The overall cell viability was con rmed by trypan blue exclusion and single cell suspensions were counted using a Countstar Rigel S2 Counter.

Library preparation and sequencing
Single-cell suspensions were loaded to 10× Chromium to capture 18070 single cell according to the manufacturer's instructions of 10× Genomics Chromium Single-Cell 3' kit (V3).The cDNA ampli cation and library construction steps were conducted via following the standard protocol.Libraries were sequenced on an Illumina NovaSeq 6000 sequencing system (paired-end multiplexing run,150bp) by LC-Bio Technology co.ltd (HangZhou,China) at a minimum depth of 20000 reads per cell.

Single-cell gene expression quanti cation, quality control, and cell type determination
The sequencing data from 10x Genomics were aligned and quanti ed using the CellRanger software package (version 3.1) against the human reference genome (GRCh38.p13).Raw gene expression matrices were imported and processed using the Seurat R package (version 3.1.1).Standard scRNA-seq ltering excludes low-quality cells with > 10 ~ 25% mitochondrial UMIs, we ltered low-quality cells by two measurements: 1) less than 500 expressed genes.2) over 25% UMIs derived from the mitochondrial genome.Gene expression matrices of the remaining 18070 cells were normalized to the total cellular UMI count.sing the LogNormalize method of the "Normalization" function of the Seurat software to calculated the expression value of genes; PCA (Principal component analysis) analysis was performed using the normalized expression value.Within all the PCs, the top 10 PCs were used to do clustering and t-SNE analysis; To nd clusters, selecting weighted Shared Nearest Neighbor (SNN) graph-based clustering method.Marker genes for each cluster were identi ed with the Wilcoxon rank-sum test via the FindAllMarkers function in Seurat.This selects markers genes which are expressed in more than 10% of the cells in a cluster and average log (Fold Change) of greater than 0.26.

Statistical analysis
Data were analyzed using R software (version 3.6.3)and GraphPad Prism (version 6).The Wilcoxon's test was used to compare the expression of m 6 A RNA methylation regulators between cancer and normal tissues.Spearman correlation analysis was performed using "corrplot" package in R. The distributions of age, sex, histological grade and TNM stage between clusters and between risk subgroups were analyzed using Chi-square test.Wilcoxon rank sum or Kruskal-Wallis rank sum test as the signi cance test depending on the number of clinicopathological features and immunotherapy response for comparison.Survival curves were plotted by using the "survival" package in R. The ROC analysis was performed for the evaluation of the AUC value in the follow-up period with the "survivalROC" package.Log-rank test was used to assess statistical signi cance.All statistical results with p < 0.05 were regarded to be statistically signi cant.

Results
The Landscape of m 6 A Methylation Regulators and TIICs in GEA To explore the important biological functions of each m 6 A RNA methylation regulator in tumorigenesis and development, we rst compare 21 m 6 A methylation regulators expression in tumor and normal samples.And the result strongly indicates that the expression of most m 6 A RNA methylation regulators was related to the occurrence of GEA (Fig. 1A-D).We speculate that the change of m 6 A RNA methylation regulators may be an intrinsic trait that can characterize individual differences.Correlation analysis was also employed to investigate the relationship between the expression level of m 6 A RNA methylation regulators of GEA.The relationship between the 21 m 6 A RNA methylation regulators is positively correlated, including YTHDF3 and KIAA1429 are most relevant in the TCGA dataset and the strongest positive correlation consisted between METTL3 and RBM15B in GSE96669 (Fig. 1E, F).The analyses presented above suggested that the high heterogeneity of expressional alteration landscape in m 6 A regulators between normal and GEA samples, indicating that the expression imbalance of m 6 A regulators played a crucial role in the GEA occurrence and progression.
Next, the difference between GEA tissues and adjacent tissues in 22 immune cells was analyzed by using the CIBERSORT algorithm in TCGA.Obviously, the proportion of immune cells in GEA tumor tissues was signi cantly different from that in normal tissues (Supplementary Figure S1A, B).We speculate that the change in the correlation of immune cells may be an internal characteristic that can re ect the external differences.Next, we investigated mutual relationship between 22 immune cells in GEA samples, and the results showed that most of the relationships between the immune cells were a negative correlation, and the Macrophages M2 and the B cells naive are most negative relevant (supplementary Figure S1C).Meanwhile, the positive correlation between NK cells resting and T cells CD4 memory activated is the most signi cant (supplementary Figure S1C).

Correlation of TME Components with Clinicopathological Characteristics and m 6 A Methylation Regulators
For determining the relationship between the proportion of immune and stromal components in TME with the clinicopathological characteristics, we analyzed the corresponding clinical information of GEA cases from the TCGA database.Stromal score was positively correlated to TMN-T stages (P = 0.009), tumor grades (P < 0.01) and tumor stages (P = 0.037) (supplementary Figure S2A); immune scores were associated with advanced tumor grades (P = 0.021), higher TNM-N level (P = 0.042), and females have higher immune scores than males (P = 0.045) (supplementary Figure S2B); ESTIMATE score showed the positive correlation with N, T classi cation of TMN stages, tumor grade and sex(P < 0.05) (supplementary FigureS2C).Therefore, these results indicated that the ratio of immune and stromal components was related to the progress of GEA, such as invasion and metastasis.Next, to explore the correlation between the high/low ratio of immune and stromal components in TME with m 6 A regulators, we found that most m 6 A regulators were highly expressed in samples with low immune and stromal scores, which indicate a special connection between TME components and m 6 A regulators (supplementary Figure S2D, E).

Consensus Clustering for m 6 A RNA Methylation Regulators Correlated with Distinct Survival, Immune Cell In ltration
As patients have a very poor prognosis, we tried to investigate GEA subtypes based on the expression of m 6 A RNA methylation regulators to explore its possible pathogenesis.We removed 39 normal samples and grouped 159 cancer tissues using the ConsensusClusterPlus package.According to the expression similarity of m 6 A RNA methylation regulators, k = 2 was the most optimum with clustering stability datasets increasing from k = 2-9 (Fig. 2A-C, supplementary Figure S3).Hence, GEA samples from the TCGA dataset were pre-classi ed into two groups (100samples in one group labeled as cluster1 and 59 samples in another group labeled as cluster2) through consensus cluster analysis.PCA was performed to elucidate the difference in transcriptional pro les between the Cluster1 and Cluster2 subgroups.And the results showed that there was a clear distinction between these two subgroups and the results of our classi cation by m 6 A RNA methylation regulators are correct (Fig. 2D).Kaplan-Meier survival analysis for the clustered samples revealed a noticeable decrease in the OS of cluster2 compared with cluster1, suggesting that the 21 methylation regulators could classify the GEA samples at prognostic level (Fig. 2E).Moreover, we discovered that cluster2 have lower PD-1 expression and most TIICs fraction was signi cantly higher in cluster 1, such as T cells CD8, T cells follicular helper, Monocytes, Macrophages M0 and Mast cells resting (P < 0.05, Fig. 2F, G).The clustering results revealed that m 6 A modi cation patterns were strongly affect TME in ltrating immune cells fraction and immunotherapy, which further con rmed our hypothesis.

Differentially Expressed Genes (DEGs) Screening and WGCNA Analysis
For exploring potential biomarkers associated with genome, m 6 A modi cation and TME cell-in ltration, we rstly identi ed 1341 DEGs in the GSE96669 dataset (supplementary table S3) and 6360 DEGs in the TCGA dataset (supplementary table S4) between GEA samples and adjacent normal samples (Fig. 3A, B).The overlapping 492 DEGs were used to further construct co-expression network (Fig. 3C, supplementary table S5).Subsequently, WGCNA was performed to construct a gene co-expression network and correlate gene modules with m 6 A cluster, immune scores and stromal scores, which were selected as trait data of WGCNA.In the case of a scale-free network and topological a hierarchical clustering tree based on dynamic hybrid cutting is established after the outlier samples are eliminated (Fig. 3D, supplementary Figure S4A).To ensure a scale-free network, we selected β = 4 (scale-free R 2 = 0.90) as a soft threshold (supplementary Figure S4B).Finally, we identi ed six gene modules correlated to m 6 A cluster, immune scores and stromal scores, respectively.
The Module trait relationships were estimated by the correlation between modules and phenotypes, which made it easier to identify highly correlated modules and phenotypes.Figure 3E showed that the brown module had signi cant relevant to immune scores (cor = 0.65, P = 2e-19).Scatter diagrams of gene signi cance was shown in Fig. 3F, which could more directly observe the signi cance of genes in brown modules related to immune scores.For the signi cant modules and key genes, GO and KEGG pathway enrichment analyses were performed.And GO analysis showed that the genes in the brown module are mainly enriched in extracellular matrix organization, extracellular structure organization and so forth (Fig. 4A).KEGG pathway enrichment analysis indicated that the genes in the brown module were mainly associated with protein digestion and absorption, PI3K-Akt signaling pathway and so forth (Fig. 4B).Based on the above analysis, we found that most enrichment pathways are associated with tumor progression.

Hub Identi cation and Its Role in Immunotherapy
The PPI network among genes in the brown module (40 nodes and 52 edges) was established by using the STRING database.Based on the MCC scores, the top ten highest-scored genes in the brown module were selected as hub genes for further analysis (Fig. 3G).Furthermore, two hub genes (COL4A1, COL5A2) in the brown module were signi cantly negatively related to the prognosis of patients with GEA (P < 0.05; Fig. 5A, B).Consistent with the above bioinformatic results, PCR analysis also revealed that COL4A1 and COL5A2 were signi cantly high expressed in tumor (supplementary Figure S5).Immunotherapies represented by PD-1/L1 and CTLA-4 blockade has undoubtedly emerged a major breakthrough in cancer therapy.Obviously, COL4A1 was signi cantly related to PD-L1 and CTLA-4 expression (Fig. 5C, D).And COL5A2 had signi cantly correlation with PD-1/L1 and CTLA-4 expression (Fig. 5E-G).The above results suggested that COL4A1 and COL5A2 may be potential biomarkers for predicting the effect of immunotherapy.Next, in the anti-PD-1 cohort (GSE78220), the signi cant clinical response to anti-PD-1 immunotherapy in patients with low COL5A2 expression compared to those with high COL5A2 expression was con rmed (Fig. 5H).In summary, our work indicated that m 6 A methylation modi cation patterns were signi cantly correlated with tumor immune phenotypes and response to immunotherapy, and COL5A2 may be a potential biomarker contribute to predicting the response to anti-PD-1 immunotherapy.

Construction and Validation of Prognostic Signatures for m 6 A RNA Methylation Regulators
For the purpose of the prognostic value of the 21 m 6 A regulators in GEA, univariate Cox regression analysis was performed based on the expression levels of the regulators from TCGA (Fig. 6A).The results indicated that 3 out of 21 tested genes including KIAA1429, HNRNPA2B1 and FMR1 are signi cantly correlated with OS (P < 0.05) and three genes are risky genes with HR > 1.For ensuring the reliability of results, PCR results also suggested that three regulators were signi cantly high expressed in tumor samples (supplementary Figure S5).To identi ed the most powerful prognostic m 6 A RNA methylation regulators, the last absolute shrinkage and selection operator (LASSO) Cox regression analysis to the 3 prognosis-related genes was conducted (Fig. 6B, C) and the coe cient of each independent prognostic gene was shown in supplementary table S6.The LASSO results showed that three regulators (KIAA1429, HNRNPA2B1 and FMR1) were the powerful prognostic factors and a risk signature was constructed.Next, the risk score was estimated according to the coe cients obtained from the LASSO analysis.To explore the prognostic role of the three-gene risk signature, we separated the GEA patients in the TCGA dataset into low and high-risk groups based on the median risk score, the results indicated that the high-risk group has a worse survival in patients with GEA (Fig. 6D).To evaluate the predictive accuracy of the risk score model, we performed the ROC curve and calculated the AUC.The AUC was 0.943 in the 5-year ROC curve for the prognostic model (Fig. 6E).Then, the univariate and multivariate Cox regression analyses results suggested the risk signature is an independent prognostic indicator (Fig. 6F, G).In addition, three m 6 A regulators were signi cantly higher expressed in cluser2 (Fig. 6H, I).This result was consistent with the fact that cluster 2 had a higher risk score, lower immune cells in ltration and poorer prognosis than cluser1.
Next, we evaluated the relative abundance of 22 TIICs for each patient within two risk groups using CIBERSORT.It was revealed that T cell follicular helper was signi cantly higher in the high-risk group compared with the low-risk group; however, Monocytes and T cells CD8 were signi cantly higher in the low-risk group compared with the high-risk group (Fig. 6J).Moreover, to screen for the possible signaling pathways and mechanisms that are signi cantly altered within high-and low-risk groups, GSEA analysis was performed with data from the TCGA cohort.As shown in Fig. 6K, RNA modi cation related pathway and cancer associated pathway were more enrichment in high risk group.The pathways enrichment analysis provided evidence of molecular mechanisms affected by the risk signature.

Correlating Public Datasets to
To investigate the cell populations and the associated molecular characteristics during the occurrence and development of GEA, we took 2 cases of GEA and 1 case of adjacent normal tissue for single-cell sequencing analysis.After quality control, 18070 single cells were clustered into 13 major clusters.
Cluster-speci c genes were used to annotate cell types with classic markers described in the supplementary table S7.Consistent with previous public dataset analysis, CD8 + T cells were enriched in normal tissue, while CD4 + T cells, Macrophages and Treg were more concentrated in tumor tissue (Fig. 7A).In addition, it seems that the proportion of CD8 + T cells and CD4 + T cells exhibited signi cant changes, indicating that CD8 + T cells may differentiate into CD4 + T cells.To further investigate this ongoing process, we performed trajectory analysis on CD8 + T cells and CD4 + T cells.Similar phenomenon was observed in the computational pipeline (Fig. 7B).Meanwhile, we noticed that CD69, T cell activation marker, was downregulated in this differentiation process, while CXCL8 and CXCL13 were upregulated (Fig. 7C).Previous work had reported that CXCL13/CXCR5 signaling axis makes pivotal contributions to the development and progression of human cancers [20].CXCL8 is mainly secreted by macrophages and contributes to the immunosuppressive microenvironment by inducing PD-L1 + macrophages in gastric cancer [21].These results indicated that the immunosuppressive microenvironment is gradually forming in this differentiation process.Then, we found that three regulators (KIAA1429, HNRNPA2B1 and FMR1) were higher expressed in tumor than normal tissue, which veri ed the reliability of the previous analysis results again (Fig. 7D).Strikingly, three m 6 A regulators originated mostly from T cells such as CD4 + T cells and Tregs in GEA (Fig. 7E), which suggested that m 6 A modi cation may participate in sustaining immunosuppressive microenvironment via control T cell differentiation.COL5A2, closely related to the expression of KIAA1429 and FMR1, also proved to be highly expressed in tumor tissues in single-cell analysis results (Fig. 7F, G).

Discussion
In this present research, we attempted to demonstrate the expression patterns, prognostic values, and effects on the TME of m 6 A regulators in GEA.Differential expression analysis found that the majority of m 6 A RNA methylation regulators were signi cantly differently expressed between adjacent normal and tumor samples.This should preliminarily suggest that those m 6 A regulators have a distinct role of in uencing GEA initiation and could be applied for the subsequent analysis.Compare to normal tissues, GEA is locally in ltrated with higher immune cell subgroups, including B cells naive, B cells memory, Plasma cells, T cells gamma delta, Macrophages M0, Macrophages M1, Dendritic cells resting, Dendritic cells active.Meanwhile, ESTIMATE algorithm-derived immune scores, stromal scores and ESTIMATE scores were applied to facilitate the quanti cation of the non-tumor components in malignancy [22].
Stromal, immune and ESTIMATE scores for tumor tissue were found to be signi cantly associated with the clinicopathologic features of the tumor such as age, differentiation grade and TNM stage.
Importantly, most m 6 A regulators' expression was signi cantly associated with immune/stromal scores.Furthermore, we constructed a comprehensive single-cell transcriptome atlas of 18070 cells from 2 samples of GEA and 1 sample of adjacent normal mucosa and revealed m 6 A regulators exert an important role in maintaining aberrant immune in ltration.Therefore, these results demonstrated that aberrant immune in ltration and m 6 A regulators expression in GEA as a tightly regulated process might play important roles in tumor development and this process had clinical importance.
We characterized the effects of distinct m 6 A methylation modi cation on different GEA subtypes by clustering m 6 A regulators.The two subtypes showed signi cant differences in patient prognosis, PD-1 expression, immune cell in ltration and RS.This suggests that two subtypes' differences are essential and re ect the heterogeneity of the immune microenvironment of GEA, which is worthy of further study.In order to investigate the expression characteristics of m 6 A methylation regulators in tumor, many studies clustered the tumor samples into different groups using consensus clustering analysis.For instance, Jing Chen et al. [23] identi ed two clusters of clear cell renal carcinoma with signi cant differences in OS and tumor stage between them based on the expression pattern of m 6 A RNA methylation regulators by means of consensus clustering.Similarly, Lilan Yi et al. [24] showed that two molecular subtypes were identi ed by consensus clustering for 15 m 6 A regulators, and two subtypes were distinct in prognosis, PD-L1 expression, immunoscore, and immune cell in ltration.However, so far, the expression of m 6 A regulators has remained elusive for typing research by consensus clustering analysis in GEA.Our research had once again identi ed the special relationship between m 6 A modi cation and tumor immune cell in ltration, which has important clinical value.
To elucidate the potential association between the genome, m 6 A modi cation pattern and TME in GEA, we systematically clustered the co-expressed genes by WGCNA analysis.This approach allowed us to identify gene modules most related to cancer immunological phenotypes.COL4A1 and COL5A2, the two hub genes in the module, in the collagen family can encode pro-alpha collagen chains that are assembled into collagen.They have a relatively high expression in tumor cells and related to the prognosis and clinicopathological factors of patients.Désert et al. [25] reported that elevated expression of COL4A1 was signi cantly correlated with tumor stage and worse overall survival in patients with hepatocellular carcinoma.Zhang et al. [26] demonstrated the abnormally high expression of COL4A1 in GC, and high expression of COL4A1 was closely correlated with primary tumor size, lymph node metastasis and distant metastasis, with the silencing of COL4A1 signi cantly inhibiting cell proliferation of GC cells in vitro.Meanwhile, elevated COL4A1 gene expression has been found to be associated with trastuzumab resistance in GC [27].Several studied had reported that COL5A2 might function as crucial role in the initiation and progression of tumor by using bioinformatics technologies [28,29].More importantly, COL5A2, was correlated with stromal scores in GC, promoted the recruitment of circulating monocytes into the TME and facilitated their differentiation into tumor-associated macrophages [30].Similarly, in our research, we found that COL4A1 and COL5A2 were signi cantly related to prognosis of GEA patients and regulate TME in ltration characteristics.Intriguingly, the expression of COL5A2 and COL4A1 was signi cantly correlated with ICB (PD-1/L1and CTLA-4) expression.Notably, COL5A2 expression was also linked to enhanced response to anti-PD-1 immunotherapy.Above results suggest that distinct m 6 A modi cation patterns of COL5A2 may establish totally different TME, which in uence immunotherapy response.
Whether m 6 A RNA methylation regulators have prognostic value in cancer is of great signi cance [31].We performed univariate and LASSO Cox regression analyses to construct a prognostic related risk signature with three m 6 A RNA methylation regulators, including KIAA1429, HNRNPA2B1 and FMR1, which divided the GEA patients into low-and high-risk groups.In the m 6 A methyltransferase complex, KIAA1429 acts as a scaffold in bridging the catalytic core components of methyltransferase complex and RNA substrates, which affect the installation of m 6 A at speci c locations [32].Ran Miao et al. [33] found that KIAA1429 could serve as an oncogene in gastric cancer by stabilizing c-Jun mRNA in an m 6 A-independent manner.HNRNPA2B1 is a nuclear reader of the m 6 A mark and has important effects on primary microRNA processing and alternative splicing.Barceló et al. [34] reported that HNRNPA2B1 acts as a regulator of KRAS-dependent tumorigenesis through the critical pancreatic ductal adenocarcinoma cells signaling pathway PI3K/AKT.FMR1 gene and consequently lack of synthesis of FMR protein (FMRP) is associated with fragile X syndrome, and FMRP plays a critical role in chromatin dynamics, RNA binding, mRNA transport, and mRNA translation [35,36].Jianhao Li et al. [37] indicated that high expression of KIAA1429 and HNRNPA2B1 were signi cantly associated with poor prognosis in osteosarcoma, and m 6 A regulators might be involved in osteosarcoma progression through humoral immune response.Francesca Zalfa et al. [38] reported that there was an association between FMRP levels and the invasive phenotype in melanoma.In accordance with previous results, we found the three-gene risk signature, KIAA1429, HNRNPA2B1 and FMR1, showed a good performance for predicting GEA patients' prognosis and immune cells in ltration characteristics.Importantly, our further study rstly revealed that the three m 6 A regulators' expression in T cells and macrophages may play a crucial role in tumor progression.Moreover, COL5A2 expression was signi cantly related to KIAA1429 and FMR1 expression.Therefore, we speculate that KIAA1429 and FMR1 regulate the number of macrophages and T cells in TME by affecting the m 6 A modi cation of COL5A2, which in turn affects the clinicopathological characteristics of patients.
The tumor microenvironment plays an essential regulatory role in tumorigenesis, and its heterogeneity can lead to multiple dimensions, including patient prognosis and therapeutic response [39][40][41].Here, we generated a single-cell transcriptome atlas of GEA sample to validate the role of m 6 A modi cation pattern on TME. 13 different cell types were identi ed in the GEA microenvironment.Notably, CD8 + T cells mostly originated from normal mucosal tissues, while Macrophages, CD4 + T cells and Treg cells were enriched in GEA tissues.Therefore, we indicated that the downregulated immunogenicity of cancer cells potentially contributes to the formation of an immunosuppressive microenvironment.Hanjie Li et al. [42] reported that a large population of CD8 + T cells showing continuous progression from an early effector "transitional" into a dysfunctional T cell state and the intensity of the dysfunctional signature is related to tumor reactivity.Moreover, CD4 + T cell help in the form of IL-21 may potentially be harnessed to bolster the formation of protective cytolytic CX3CR1 + CD8 + T cells and improve control over tumor progression [43].Our monocle analysis veri ed again that CD4 + T cells regulate CD8 + T cell differentiation in tumor progression of GEA.m 6 A RNA modi cation controls the differentiation of naive T cells and sustains the suppressive functions of Tregs [5,44].Then, our work revealed that three m 6 A regulators are highly expressed in CD4 + T cells, CD8 + T cells, Tregs and Macrophages, which was consistent with previous works.In short, we rst discovered that KIAA1429, HNRNPA2B1 and FMR1 regulate T cell differentiation in the GEA microenvironment, which may provide new targets to optimize immunotherapy.

Conclusions
In summary, this study systematically evaluated prognostic value, role in TME, correlation of immunotherapy, and potential regulatory mechanisms of m 6 A RNA regulators in GEA.A notable result was the identi cation of three m 6 A RNA methylation regulatory genes risk signature, KIAA1429, HNRNPA2B1 and FMR1, as the essential determinant of GEA patients' prognosis and immune cells in ltration characteristics.The distinct m 6 A modi cation patterns of COL5A2 were found to be linked to enhanced response to anti-PD-1 immunotherapy.The m 6 A RNA methylation may be involved in regulation of TME via multiple RNA and cancer-related processed such as PI3K/AKT signaling pathway.
Our work highlights the importance of CD8 + T cells differentiation in antitumor immunity and ICB-based immunotherapy for GEA at single-cell level.The information from this study makes important contributions to our understanding of m6A RNA regulators and TME of GEA and may help the development of a new generation of immune therapeutics and precision treatment in GEA.Abbreviations m 6 A: RNA N6-methyladenosine; TME: Tumor Microenvironment; GEA: Gastroesophageal Adenocarcinomas; TCGA: The Cancer Genome Atlas; GEO: Gene Expression Omnibus; DEGs: Differentially Expressed Genes; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; WGCNA: Weighted Gene Co-expression Network Analysis; PPI: Protein-Protein Interaction; TIICs: Tumor-In ltrating Immune Cells; scRNA-seq: single-cell RNA sequencing.

Declarations
Ethics approval and consent to participate All experimental protocols were approved by the Ethics Committee of the Jinan Central Hospital (Jinan, China) and performed in accordance with the relevant guidelines and regulations.Written informed consent was obtained from all patients.Consent for publication Not applicable.Availability of data and material The datasets analysed during the current study are available in the GEO repository (https://www.ncbi.nlm.nih.gov/geo/) and TCGA repository (https://portal.gdc.cancer.gov/).The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.Competing interests No potential con icts of interest were disclosed.Funding This study was nancially supported by the National Natural Science Foundation of China (Nos.31671468 and 31970728), the Academic Promotion Programm of Shandong First Medical University (No. 2019QL024), the Shandong Provincial Natural Science Foundation of China (Nos.ZR2016HM15, ZR2020LZL005, ZR2020MH050).Authors' Contributions DL contributed to conceptualization, resources, data curation, validation, methodology, writing-original draft.JZ contributed to investigation, writing-review and editing.ZW contributed to data curation, formal analysis, methodology.YF contributed to resources, formal analysis, methodology.MY contributed to software, formal analysis.XM contributed to data curation, visualization.YJ contributed to conceptualization, supervision, writing-review and editing, funding acquisition and validation.YW contributed to conceptualization, supervision, funding acquisition, validation, visualization.