Check for updates
Elucidating the clinical and immunological value of m6A regulator- mediated methylation modification patterns in adrenocortical carcinoma
WENHAO XU1,#; HAOMING LI2,#; YASIR HAMEED3; MOSTAFA A. ABDEL-MAKSOUD4; SAEEDAH MUSAED ALMUTAIRI4; AYMAN MUBARAK4; MOHAMMED AUFY5; WAEL ALTURAIKI6; ABDULAZIZ J. ALSHALANI6; AYMAN M. MAHMOUD7, *; CHEN LI8,*
1 Department of Urology, Urological Surgery Research Institute, First Affiliated Hospital, Army Medical University (Third Military Medical University), Chongqing, China
2 Department of Urology, Affiliated Hospital of Guilin Medical University, Guilin, China
3 Department of Applied Biological Sciences, Tokyo University of Science, Tokyo, Japan
4 Department of Botany and Microbiology, College of Science, King Saud University, Riyadh, Saudi Arabia
5 Department of Pharmaceutical Sciences, Division of Pharmacology and Toxicology, University of Vienna, Vienna, Austria
6 Department of Medical Laboratory Sciences, College of Applied Medical Sciences, King Saud University, Riyadh, Saudi Arabia
7 Department of Life Sciences, Faculty of Science and Engineering, Manchester Metropolitan University, Manchester, UK
8 Department of Biology, Chemistry, Pharmacy, Free University of Berlin, Berlin, Germany
Key words: m6A, Tumor microenvironment, Adrenocortical carcinoma, Immunotherapy, Prognosis
Abstract: N6-methyladenosine methylation (m6A) is a common type of epigenetic alteration that prominently affects the prognosis of tumor patients. However, it is unknown how the m6A regulator affects the tumor microenvironment (TME) cell infiltration in adrenocortical carcinoma (ACC) and how it affects the prognosis of ACC patients yet. The m6A alteration patterns of 112 ACC patients were evaluated, furthermore, the association with immune infiltration cell features was investigated. The unsupervised clustering method was applied to typify the m6A alteration patterns of ACC patients. The principal component analysis (PCA) technique was taken to create the m6A score to assess the alteration pattern in specific malignancies. We found two independent patterns of m6A alteration in ACC patients. The TME cell infiltration features were significantly in accordance with phenotypes of tumor immune-inflamed and immune desert in both patterns. The m6Ascore also served as an independent predictive factor in ACC patients. The somatic copy number variation (CNV) and patients prognosis can be predicted by m6A alteration patterns. Moreover, the ACC patients with high m6A scores had better overall survival (OS) and higher efficiency in immune checkpoint blockade therapy. Our work demonstrated the significance of m6A alteration to the ACC patients immunotherapy. The individual m6A alteration patterns analysis might contribute to ACC patients prognosis prediction and immunotherapy choice.
Introduction
The unusual endocrine tumor generated from the adrenocortical is termed adrenocortical carcinoma (ACC). According to the Netherlands Cancer Registry, between 1993 and 2010, there were about 0.5 and 2.0 incidences of
ACC per million individuals [1]. In addition to being uncommon, patients with ACC also have a poor 5-year overall survival (OS) rate of approximately 15% to 44% [2,3]. The prognosis of ACC is significantly affected by the disease stage upon diagnosis, with a 5-year survival rate of 80% for stage I and 13% for stage IV [4]. The sole therapy for early to mid-stage ACC is surgical operation, however, even after the tumor has been completely excised, generally 19%-34% local recurrence rate exists [5,6]. The advanced ACC patients who are unable to receive surgery or with metastases might benefit from other therapies, but it is achieved by oral chemotherapy-mitotane therapy currently.
*Address correspondence to: Ayman M. Mahmoud,
A.Mahmoud@mmu.ac.uk; Chen Li, chen.li.scholar@gmail.com
“These authors have contributed equally to this work Received: 17 February 2023; Accepted: 26 April 2023; Published: 21 July 2023
www.techscience.com/journal/or
CC
, BY
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Therefore, understanding the m6A alteration method is essential for providing a novel strategy for ACC patients.
The m6A alteration, insertion of a methyl group to the N6 position of adenylic acid, was initially identified in the 1970s [7], since m6A alteration was shown to be dynamically reversible and obesity-associated protein was discovered to contain m6A demethylase function [8]. The methyltransferase complex, demethylase, and the related readers, which are referred to as “writers,” “erasers,” and “readers,” control the associated alteration of m6A. Numerous studies have interpreted the role of m6A methylation alteration in tumor invasion. In breast cancer, writers (METTL3) and erasers (ALKBH5) might promote tumor cell replication and migration by controlling the degree of m6A methylation [9]. However, the effect of m6A alteration in ACC patients is not well known yet.
Tumor microenvironment (TME), comprised of extracellular matrix, various soluble substances, immune cells, stromal cells, endothelium, and cancerous cells [10], mitigates tumor deterioration and affects the consequences of immunotherapeutic inhibitors in clinic [11,12]. The m6A alteration has been reported with close relation to immune cell infiltration that controls T cell homeostasis by targeting the IL-7/STAT5/SOCS pathways [13], and METTL3 drives M1 macrophage polarization by directly methylating STAT1 mRNA [14].
The noteworthy limitation of previous studies is merely a small amount of m6A regulators were analyzed, the systematic investigation is required, thus, in current work, we integrated various m6A alternations associated with TME cell infiltration features in ACC patients. ACC subjects from the TCGA and GEO databases were combined, and m6A alteration patterns were assessed. Based on the comparison of m6A alteration patterns and TME immune cell infiltration, two distinct m6A alteration patterns were identified, more importantly, these two patterns were very consistent with immune- inflamed and immune-desert phenotypes, respectively. This highlights that m6A alteration serves a vital role in the TME of ACC patients. Subsequently, the complicated biological processes were determined by enrichment analysis. Finally, m6A alteration was quantified by the created m6A score model which might develop a precise approach to enhance therapeutic utility.
Materials and Methods
Source and processing of the dataset for adrenocortical carcinoma
The workflow was shown in Fig. 1. In the current study, TCGA-ACC (79 subjects with clinic survival information among total of 92 subjects) and GSE33371 from the GEO database (23 subjects with clinic survival information among total of 33 subjects) were analyzed. 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
TCGA-ACC (79 subjects with clinic survival information among total of 92 subjects)
GSE33371 (23 subjects with clinic survival information among total of 33 subjects)
Identification m6A modification patterns by unsupervised clustering (m6A cluster A and B)
TME infiltration estimation
1
m6A associated phenotype genes identification
DFGs KEGG and Go analysis
Identification m6A associated phenotype genes by Consensus clustering algorithm (Gene cluster A, B and C)
m6A methylation pattern quantification by PCA algorithm
I
Correlation between m6A score and TME infiltration
Correlation between m6A score and tumor somatic mutation
The role of m6A in the immunotherapy
transformed into transcripts per kilobase million (TPM) values. Batch effects from non-biological technical biases were corrected using the “ComBat” algorithm of “sva” package. To further analyze copy number variation (CNV), we obtained the TCGA database’s basic nucleotide variation data. The GSE33371 microarray data from the Genetic Code U133 Plus 2.0 Array from Affymetrix. For microarray data from Affymetrix®, the raw “CEL” files were downloaded and adopted a robust multiarray averaging method with the affy and simplify packages to perform background adjustment and quantile normalization.
Unsupervised clustering of 23 m6A regulators
23 well-known m6A regulators, including 13 readers, 8 writers, and 2 erasers were obtained from previous references [15]. The m6A-related genomic matrix from the combined sample according to 23 m6A regulators was extracted by ‘limma’ package. The unsupervised cluster analysis was applied to categorize the patients, we found several m6A alteration patterns based on the expression of 23 m6A regulators. The Euclidean distance was chosen as the clustering metric in the sample clustering, and the K- means clustering approach was used to accomplish the clustering goal.
Gene set variation analysis (GSVA)
The comparison of the variations across biochemical functions according to m6A alteration patterns was performed by the unsupervised, nonparametric “GSVA” R package. Specifically, “c2.cp.kegg.v7.4.symbols.” gene set was examined. The difference is obvious when adjusted p < 0.05. To create the differential path heatmap, the first 20 differential pathways are chosen. The ‘limma’, ‘GSEABase’, ‘GSVA’, and ‘pheatmap’ packages, were taken to complete the above analyses.
TME cell infiltration estimation
Single-sample gene set enrichment analysis algorithm (ssGSEA) was used to determine the relative abundance of each cell infiltrate in the sample, including activated B cells, CD56 Natural killer cells, eosinophil monocytes, and many others. A box graphic was used to illustrate the distinct immunological invading cells.
DEG screening for m6A alteration patterns
The differentially expressed genes (DEGs) according to distinct m6A alteration patterns were determined by the empirical Bayesian, p < 0.001 was considered as significance. The cellular component, biological process, and molecular function of DEGs were examined by the Gene Ontology (GO) database. Additionally, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to conduct the pathway analysis of DEGs.
The creation of the m6A gene signature
We developed a scoring method, termed as m6Ascore, to assess the m6A alteration pattern of individual ACC patient. The detailed steps were as follows. Cox regression analysis was utilized to identify DEGs with prognostic value, and then the unsupervised cluster was used to analyze different clusters with significant prognosis. Subsequently, m6A relevant gene signature was constructed by the principal component analysis (PCA), and both principal component 1 and 2 were selected. The advantage of this method is to focus on the score in the set with the largest block of well- correlated (or anti-correlated) genes, while down-weighting contributions from uncorrelated genes. Finally, the m6Ascore was defined by a method similar to Gene-gene interaction (GGI) [16,17]:
m6Ascore ⇒(PC1; + PC2i)
where, i stands for the expression of a gene associated to m6A.
The function of the m6A score
The associations of m6Ascore and immune infiltrating cells, tumor mutation burden (TMB), m6A alteration pattern, m6A gene cluster type, and genes associated with immunotherapy were investigated to demonstrate the clinical value of m6Ascore, p < 0.05 was regarded as significant.
Statistical analysis
The Kruskal-Wallis test was used for the nonparametric test, and one-way ANOVA was used for the parametric test. The association between m6Ascore and TMB was determined by Spearman correlation analysis. The survival was constructed by the Kaplan-Meier curve, the significance was evaluated by the log-rank test. The subjects were evaluated by PCA, paired with patient survival, and then separated into two groups using the optimal cutoff value given by the ‘Survminer’ R software program. The hazard ratios (HR) and independent prognostic variables of the m6A regulatory gene and m6A phenotypic associated gene were calculated by Cox regression. The ‘RCircos’ R software was used to create the position circle map of the m6A regulator CNV on
the 23 chromosomes. The waterfall was created by “maftools” package, and the box diagram was created by the ‘ggpubr’ package. R 4.0.5 was used to do all analyses, and p < 0.05 was regarded as significant.
Results
The genomic landscape of m6A regulatory mutation in ACC 23 m6A regulators total, eight writers, two erasers, and thirteen readers, were used in this study. The m6A regulatory gene was altered in 8 of the 92 samples from TCGA cohort, yielding an 8.7% mutation frequency (Fig. 2A). The nonsense mutation is the second most frequent mutation after missense mutations. ZC3H13, YTHDF3, METTL14, WTAP, RMB15, RMB15B, YTHDC2, YTHDF1, IGFBP1, IGFBP3, and FTO were among the m6A regulatory genes with mutations. The univariate Cox regression analysis revealed 18 m6A regulators were related to the prognosis of ACC patients. The network node diagram showed the relationship between m6A regulators and prognosis of ACC patients (Fig. 2B and Suppl. Table S1). According to Fig. 2C, the CNV mutation of m6A regulator could be found on 23 pairs of chromosomes. Amplification accounted for the majority of CNV alterations of 23 m6A regulatory genes, while CNV deletions were in METTL14, METTL16, RMB15, RMB15B, LRPPRC, and YTHDF1 (Fig. 2D).
23 regulators mediate the patterns of m6A alteration
The clinical characteristics of the subjects from TCGA-ACC and GSE33371 that were included in this work are shown in Suppl. Table S2. Based on the expression levels of m6A- related genes, the survival curves of ACC patients are shown in Fig. 3, there exists a strong association between writers, erasers, and readers as well as between the expression of the m6A regulator in the same functional category. Unsupervised cluster analysis was used to categorize, and the expression of 23 m6A regulators was used to distinguish between two distinct m6A alteration patterns, 98 Pattern A cases and 14 Pattern B. These two patterns are termed m6A cluster A and m6A cluster B. (Figs. 4A-4E and Suppl. Table S3). The prognostic study of the two types of m6A clusters revealed that m6A cluster B had a considerable survival advantage (Fig. 4F).
TME cell infiltration variations under two m6A alteration patterns
GSVA enrichment analysis was applied to examine the biological differences between the two types of m6A clusters. The signal transduction pathways RNA degradation, aminoacyl tRNA biosynthesis, basal transcription factors, and neurotrophic signaling pathways were all considerably enriched in the m6A cluster A, as shown in Fig. 5A and Suppl. Table S4. The arachidonic acid metabolism, drug metabolism, cytochrome p450 metabolism, and linoleic acid metabolism pathways were all considerably concentrated in m6A cluster B. The m6A cluster B had much more immune cell infiltration than m6A cluster A, particularly activated CD8 T cells, CD56 bright
(A)
(B)
Altered in 8 (8.7%) of 92 samples.
1776
YTHOF!
YTHDF3
YTHOFI
YTHDC2
0
2
HNUAPC
· erasers
ZC3H13
0
2%
THDCI
YTHDF3
METTL14
2%
1%
FMRI
· readers
WTAP
RBM15
1%
· writers
1%
RBM15B
1%
ALKBHS
YTHDC2
YTHDF1
1%
LREPR
IGFBP1
1%
1%
IGFBP3
1%
FTO
FTO
METTL3
1%
METTL16
0% 0%
VIRMA
YTHDC1
0% 0%
YTHDF2
HNRNPC
0%
FMR1
0%
0%
MISB
LRPPRC
RNPA2B1
0%
0%
IGEBPZ
. Risk factors
IGFBP2
RBMX
0%
0%
· Favorable factors
ALKBH5
0%
METIL3
· C>T = T>A
ZC3H13
· C>G = T>C
METTLI
C>A . T>G
WIAP
VIRMA
· Missense_Mutation = Multi_Hit Nonsense_Mutation
Postive correlation with P<0.0001
Negative correlation with P <0.0001
Cox test, p value
1e-04
0.001
0.01
0.05
1
(C)
(D)
×
<
-
10
. GAIN . LOSS
21
5
30
2
19
RBMX
FMR1
YTHDF2
RBM15
8
18
YTHDF1
LRPPRC
IGFBP2
CNV.frequency(%)
17
RBM158
3
ALKBH5
6
16
METTL16
FTO
15
YTHDC1
METTL14
4
METTL3
HNRNPC
4
14
ZC3H13
VIRMA
YTHDF3
HNRNPA2B1 WTAP
YTHDC2
5
13
IGFBP3
SFBP
2
12
6
~
0
1
YTHDC2
METTL3
HNRNPC
ALKBH5
FMR1
WTAP
FTO
RBMX
IGFBP2
METTL16
YTHDF2
RBM15
YTHDC1
IGFBP
IGFBP3
VIRMA
ZC3H13
LRPPRC
METTL14
HNRNPA2B1
YTHDF3
YTHDF1
RBM15B
0
@
00
NK cells, monocytes, neutrophils, and type 17 T helper cells (Fig. 5B). The immune-inflamed phenotype was defined by the presence of T cells that express CD4 and CD8 in the tumor parenchyma, in conjunction with myeloid cells and monocytes [18]. Additionally, individuals with the immune- inflamed phenotype have a better prognosis since they have better reaction to immunotherapy, which is supported by our findings. The transcriptome profiles of the m6A alteration pattern were then subjected to PCA, and the results revealed a substantial difference between the m6A cluster A and m6A cluster B in TME cell infiltration characterizations (Fig. 5C). Cluster B was designated as an immunological-inflamed phenotype, which is defined by immune activation and adaptive immune cell infiltration, whereas Cluster A was designated as an immune-desert phenotype, which is characterized by the suppression of immunity.
Genetic function evaluation of DEGs and creation of m6A gene clusters
Total of 8874 DEGs related to the m6A alteration pattern was identified by “limmma” package, followed by GO and KEGG enrichment (Suppl. Tables S5 and S6). The DEFs were enriched in the catabolic process of mRNA, the cell- substrate junction, and transcription coregulator activity biological processes, which demonstrate the non-negligible effect of m6A alteration in CC, BP, and MF regulation (Fig. 6A). The KEGG analysis of DEGs revealed variations in the mTOR signaling route, the insulin signaling pathway, the mRNA surveillance pathway, and other pathways (Fig. 6B). Further unsupervised cluster analysis of DEGs was conducted, the ACC patients were divided into three distinct genomic subtypes, dubbed m6A gene clusters A-C, as shown in Figs. 6C-6F. The majority of the follow-up for patients with m6A gene cluster A focused on mortality. The
m6A METHYLATION PATTERNS IN ADRENOCORTICAL CARCINOMA
823
(A)
(B)
(C)
ALICOHS
high
low
FMR1
high
low
HNRNPA201
high
low
1.00
1.00
1.00+
0.75
0.75
0.75
0.50
0.50
0.50
0.25
p=0.039
0.25
p=0.045
0.25
p<0.001
Overall survival
Overall survival
Overall survival
0.00
0.00
0.00
0
1
2
5
4
5
6
7
8
$
10
11
12
3
14
15
0
1
2
4
5
6
3
8
9
10
Y
17
7
14
15
0
1
5
4
5
6
8
9
10
11
12
13
14
15
Time(years)
Time(years)
Time(years)
74
54
O
nian
A
9
O
9
O
co
AS
C
jo
A
S
68
La
&N
O
Se
a
NN
EN
INN
0
1
2
3
4
5
6
7
8
9
10
11
12
3
14
15
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7
8
9
10
11
12
1
3
14
15
FMR1
ALKBH5
Time(years)
Time(years)
Time(years)
HNRNPA281
(D)
HNANPC
high
low
E)
IGFBP1
high
F
IGFBP2
high
low
1.00-
1.00
1.00-
0.75
0.75
0.75
0,50
0.50
0.50
0.25
p<0.001
0.25
p=0.007
0.25
p=0.004
Overall survival
Overall survival
Overall survival
0.00
0.00
0.00
0
1
A
3
4
5
7 8
9
10
11
12
13
14
15
0
1
2
3
4
5
7
8
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7 8
9
10
11
12
3
4
14
15
Time(years)
Time(years)
Time(years)
E
A
Y
LE
g
33
9
es
9
9
9
O
3
O
0
O
o
0
1
2
3
4
5
6
7
B
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7
1
8
9
10
11
12
13
14
15
AS
BU
Pour
O
IGFBP2
IGFBP1
Time(years)
Time(years)
Time(years)
HNRNPC
G)
IGFBP3
high
kow
(H)
LRPPRC
high
(I)
low
METTLE
high
kow
1.00
1.00
1.00+
0.75
0.75
0.75
0.50
0.50
0.50
0.25
p<0.001
0.25
p=0.004
0.25
p=0.001
Overall survival
Overall survival
Overall survival
0.00
0.00
0.00
0
1
2
S
4
5
6
7
9
9
10
1
12
13
14
15
0
2
3
4
5
&
8
9
10
11
12
13
14
15
0
2
3
4
5
8
8
9
10
11
12
13
1
15
Time(years)
Time(years)
Time(years)
P
02.00
D
A
d
S
MA
A
Pa
1
O
I
2
359
2-
0G
6
GN
NO
83
GM
COCA
3300
56
8
80
85
CHU
GH
0
1
2
1
3
4
5
8
0
10
11
2
13
14
15
0
1
2
3
4
5
1
7
8
9
10
11
12
13
14
15
0
1
2
3
4
5
6
8
9
10
11
12
13
14
15
IGFBP3
METTL3
LRPPRC
Time(years)
Time(years)
Time(years)
(J)
METTL14
K)
RBM15
(L)
high
low
high
low
RBM15B
high
low
1.00
1.00
1.00
0.75
0.75
0.75
.50
0.50
0.50
0.25
p=0.031
0.25
p<0.001
0.25
p=0.013
Overall survival
Overall survival
Overall survival
0.00
0.00
0.00
0
1
2
3
4
5
6
7
8
9
10
1
12
13
14
15
0
1
2
6
7
8
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7
8
9
10
11
12
1
13
14
15
En
o
%
Time(years)
Time(years)
Time(years)
X
e
C
DEN
9
5
2
O
0
5
3
C
قرب العـ
16
DO
00
0
1
2
3
4
6
6
7
8
9
10
1
11
12
13
14
15
0
1
2
A
£
1
TN
ZO
5
8
a
o
3
6
10
11
12
13
14
15
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
AN
1
A
15
T
I
RBM15
RBM158
Time(years)
Time(years)
Time(years)
METTL14
(M)
WTAP
(N)
low
YTHDC1
high
ow
(O)
high
YTHDC2
high
low
1.00-
1.00-
1.00-
0.75
0.75
0.75
0.50
0.50
0.50
0.25
p=0.001
0.25
p=0.006
0.25
p=0.022
Overall survival
Overall survival
Overall survival
0.00
0.00
0.00
FIGURE 3. The prognostic analysis of ACC
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
0
1
A
5
$
3
8
9
10
11
12
13
1
15
0
2
3
4
3
$
1
8
9
10
11
12
13
14
15
Time(years)
Time(years)
Time(years)
patients was based on the expression levels
e
high
low
BE
54
Z
9
99
3
14
10
2
Pun
S
1
0
0
of m6A-related genes Kaplan-Meier curves
C
IN
C
3
O
no
SE
cod
NN
NY
~00
26
an
SE
81
8
BS
La
0
1
2
3
4
5
6
9
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
0
2
3
4
5
6
7
8
10
11
12
13
14
15
WTAP
YTHDC2
YTHDC1
Time(years)
Time(years)
Time(years)
(P)
(Q)
for patients with high or low ALKBH5
YTHOF1
high
low
YTHOF2
high
low
YTHOF3
high
low
(A), FMR1 (B), HNRNPA2B1 (C), IGFBP3
1.00
1.00
1.00-
(D), LRPPRC (E), METTL3 (F), YTHDF3
@ 0.75
0.75
0.75
(G), YTHDF2 (H), YTHDF1 (I), HNRNPC
0.50
0.50
0.50
(J), IGFBP1 (K), IGFBP2 (L), METTL14
0.25
p=0.032
0.25
p<0.001
0.25
p=0.047
Overall survival
Overall survival
Overall survival
(M), RBM15 (N), RBM15B (O), YTHDC2
0.00
0.00
0.00
0
1
2
3
5
6
7
0
8
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7
8
9
10
11
12
13
4
15
0
1
2
3
4
5
£
8
9
,
10
11
12
13
14
15
Time(years)
Time(years)
Time(years)
(P), YTHDC1 (Q), WTAP (R) expression
in the merge cohort (TCGA-ACC and
high
3
A
9
12
h
5
2
0
0
O
14
A
4
U
2
I
C
A
A
DJ
0
Ou
26
0-
2
98
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
0
1
2
3
Q
4
5
6
7
8
9
10
11
12
13
14
15
0
1
2
3
4
5
6
7
/
8
10
11
12
3
14
15
GSE33371).
YTHDF3
YTHDF2
Time(years)
YTHDF1
Time(years)
Time(years)
features of patients with m6A gene cluster B were comparable
Among three m6A gene clusters, we found there was a
to those of patients with m6A gene cluster A. Contrarily, the
significant difference 23 m6A regulators expressions. FTO,
majority of the follow-up for patients with m6A gene cluster C
FMR1, and METTL14 were highly expressed in m6A gene
characteristics focused on live monitoring because they were
cluster B, while the majority of the m6A-related genes were
almost all in the m6A cluster B pattern (Fig. 6G). The
highly expressed in m6A gene cluster A (Fig. 6I).
altered genomic phenotype of m6A was significantly
correlated with the OS rate of ACC patients, according to
Correlation between m6A score and phenotypes associated with
the Kaplan-Meier curve. The prognosis for patients with
m6A
m6A gene cluster A was the worst, whereas the prognosis
The m6A mutation pattern in ACC patients was determined
for patients with m6A gene cluster B was the best (Fig. 6H).
by the created the m6Ascore system. The association
(A)
(B)
(C)
(D)
consensus matrix k=2
consensus matrix k=3
consensus matrix k=4
Delta area
0
3
0
O
-
A
relative change in area under CDF curve
:
=
0
3
3
8
0
2
3
4
5
6
7
8
9
(E)
(F)
k
Fustat
Fustat
m6Acluster
Gender
Project
Alive
1.00-
A
môAcluster
4
Dead
Gender
B
METTL3
2
Female Male
METTL14
Project
0.75
WTAP
0
GSE33371
Survival probability
VIRMA
TOGA
m6Acluster
ZC3H13
A
RBM15
B
0.50
RBM158
-4
YTHDC1
YTHDC2
0.25
YTHOF1
p=0.033
YTHOF2
YTHOF3
HNRNPC
0.00
FMR1
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
LRPPRC
Time(years)
HNRNPA2B1
IGFBP1
m6Acluster
Number at risk
IGFBP2
A
88
75
56
43
30
24
15
11
8
7
5
3
3
1
1
1
IGFBP3
B
14
14
11
8
7
6
6
4
3
2
1
1
1
0
0
0
FTO
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ALKBH5
Time(years)
between the m6A cluster, m6A gene cluster, m6Ascore, and the survival status of specific patients was shown in the Sankey diagram. The m6Ascore was created according to the m6A alteration pattern of DEGs in the PCA (Fig. 7A and Suppl. Table S7). Additionally, the association between m6Ascore immune cell infiltration in TME was investigated, and most immune cells, including activated B cells, active CD8 T cells, and type 17 T helper cells, were substantially associated with m6Ascore (Fig. 7B). In the m6A altered model with the high abundance of immune cell infiltration, the higher m6Ascore was associated with inflammation in the immune system, while the lower m6Ascore was associated with the immunological desert in the low abundance of immune cell infiltration model. The Kruskal- Wallis test revealed the median m6Ascore of m6A gene cluster B was higher than that of m6A gene cluster A. Additionally, there were notable variations in m6Ascore amongst m6A gene clusters (Fig. 7C). The three-stage grading of patients could be made by the m6Ascore, that grade the prognosis of patients more precisely. The median m6A score of the m6A cluster B was higher than m6A cluster A (Fig. 7D). When the population was divided by gender, the men subgroup showed higher m6Ascore than the women (Figs. 7E and 7F). The low m6Ascore subgroup had higher percentage of deaths (Figs. 7G and 7H). Additionally, the patient prognosis was accessed by m6Ascore, high m6Ascore patients had longer survival
(Fig. 7I). Immune cell infiltration was more pronounced in patients with high m6Ascores which might lead to longer survival. This phenotype also existed in the population divided by gender (Figs. 7J and 7K).
m6A alteration in tumor somatic mutation features
The association between m6Ascore and tumor somatic mutations was assessed by “maftools” package, the substantial difference in TMB between the groups with high and low m6Ascores (Fig. 8A). Additional assay revealed the percentage of TMB decreased as the m6Ascore increased that a negative association between the m6Ascore and TMB (Fig. 8B). As shown in the waterfall chart (Fig. 8C), 38 cases altered among 44 samples (86.36%) in low m6Ascore group, and 19 cases altered among 35 samples (54.29%) in high m6Ascore group (Fig. 8D). ACC patients with low TMB have longer survival (Fig. 8E), furthermore, low TMB with high m6Ascore group has the best prognosis, while high TMB with low m6Ascore has the worst prognosis (Fig. 8F).
m6A alteration pattern’s function in immunological checkpoints
Clinical immunotherapy, including PD-1/L1, CTLA4, and NIK blocker is crucial in malignancies management. The responses of groups with high and low m6Ascores to immune checkpoint blockade therapy were assessed in the merging cohort (TCGA-ACC and GSE33371). The ACC
(A)
(B)
m6Acluster
m6Acluster
Project
2
A
KEGG_ARACHIDONIC_ACID_METABOLISM
B
KEGG_NEUROACTIVE_LIGAND_RECEPTOR_INTERACTION
1
Project
3
KEGG_OLFACTORY_TRANSDUCTION
0
GSE33371
TCGA
KEGG_LINOLEIC_ACID_METABOLISM
KEGG_RETINOL_METABOLISM
-1
2
KEGG_METABOLISM_OF_XENOBIOTICS_BY_CYTOCHROME_P450
-2
KEGG_DRUG_METABOLISM_CYTOCHROME_P450
1
KEGG_TASTE_TRANSDUCTION
m6Acluster
KEGG_MATURITY_ONSET_DIABETES_OF_THE_YOUNG
PC2
A
KEGG_SNARE_INTERACTIONS_IN_VESICULAR_TRANSPORT
0
..
B.
B
KEGG_NEUROTROPHIN_SIGNALING_PATHWAY
KEGG_ENDOCYTOSIS
KEGG_AMINOACYL_TRNA_BIOSYNTHESIS
-1
KEGG_LYSINE_DEGRADATION
KEGG_UBIQUITIN_MEDIATED_PROTEOLYSIS
-2.
KEGG_BASAL_TRANSCRIPTION_FACTORS
KEGG_PROTEIN_EXPORT
KEGG_RNA_POLYMERASE
-5
0
5
10
KEGG_RNA_DEGRADATION
PC1
15
KEGG_SPLICEOSOME
(C)
m6Acluster
A
B
1.00
ns
ns
ns
ns
ns
ns
ns
ns
ns
ns
ns
**
ns
ns
*
ns
ns
ns
ns
*
ns
0.75
…
…
A
Immune infiltration
0.50
0.25
0.00
Activated.B.cellna
Activated.CD4.T.cellna
Activated.CD8.T.cellna
Activated.dendritic.cellna
CD56bright.natural.killer.cellna
CD56dim.natural.killer.cellna
Eosinophilna
Gamma.delta.T.cellna
Immature .. B.cellna
Immature.dendritic.cellna
MDSCna
Macrophagena
Mast.cellna
Monocytena
Natural.killer.T.cellna
Natural.killer.cellna
Neutrophilna
Plasmacytoid.dendritic.cellna
Regulatory.T.cellna
T.follicular.helper.cellna
Type.1.T.helper.cellna
Type.17.T.helper.cellna
Type.2.T.helper.cellna
Patients with high m6Ascores had substantially increased levels of PD-L1 and CTLA4 expression (p < 0.01), suggesting a possible therapeutic response (Figs. 9A and 9B). However, the therapeutic potential of the NIK blocker might be weaker, since the high m6Ascore group had much lower NIK expression than the low m6Ascore group (Fig. 9C).
Discussion
On RNA adenine, including mRNA and lncRNA, the methylation alteration known as m6A may occur [7]. Similar to another epigenetic control of DNA and histone alteration, the m6A alteration is dynamically reversible in mammalian cells [19]. The m6A RNA alteration affected RNA at many points during its life cycle, including RNA
processing, nuclear export regulation, translation, and RNA degradation [20]. Evidence suggests that m6A alteration is intimately linked to the development of tumors. The expression level of m6A-related protein, a key regulator of tumor formation and progression, often defines the pathological course of the tumor [21]. m6A alteration has been shown to promote or inhibit cancer in several tumor types, including breast cancer [9], lung cancer [22], and acute myeloid leukemia [23], while preventing glioblastoma [24]. The impact of ACC m6A alteration mode on a patient’s prognosis is unclear since it is an uncommon endocrine malignant tumor. We combined the effects of 23 m6A regulators, in contrast to most studies that concentrated on one or two m6A regulators, to shed light on the function of two distinct m6A alteration patterns.
(A)
(B)
ncRNA metabolic process
ribonucleoprotein complex biogenesis
Herpes simplex virus 1 infection
proteasomal protein catabolic Process
Amyotrophic lateral sclerosis
RONA SplICING
Salmonella infection
ubiquitin-dependent protein catabolic process
Endocytosis
ncANA processing-
6
Shigellosis-
RNA catabolic process
Human T-cell leukemia virus 1 infection-
MRNA catabolic process-
Protein processing in endop aprile transport!
ribosome biogenesis macroautophagy
Count
RNA transport
·
150
Ubiquitin mediated proteolysis
qvalue
mitochondrial matrice-
.
200
Autophagy ” animal
nuclear speck
O
250
Spliceosome
mTOR signaling pathway
2.50-07
cell-substrate junction
300
Cellular senescence
5.00-07
focal adhesion-
350
Ribosome-
7.50-07
spindle
8
Cell cycle
Insulin signaling pathway
1.00-06
somal region
transferase complex, transferring
qvalue
Lysosome
Count
phosphorus-containing groups spliceosomal complex
9.6119050-50
2079618e-19
Apoptosis
ribosomal subunit
FoxO signaling pathway-
100
Neurotrophin signaling pathway
transcription coregulator activity
4.159635e-19
MRNA surveillance pathway
200
ubiquitin-like protein transferase acavy1
6.2394536-19
RNA degradation-
300
ubiquitin-protein activity
Chronic myeloid leukemia-
catalytic activity, acting on RNA-
0.31927De-19
Pancreatic cancer-
DNA-binding transcription factor binding
ubiquitin-like protein ligase binding
5
Mitophagy = animal
ubiquitin binding-
Renal cell carcinoma
RNA polymerase II-specific DNA-binding transcription factor binding
Nucleotide excision repair
transcription coactivator activity-
DNA replication
Base excision repair
histone binding-
Autophagy = other
·
0.02
0.03 0.04
(C)
(D)
GeneRatio
0.02
0.04
0.06
0.08
GeneRatio
(E)
consensus matrix km2
gensenaus matrix kad
consensus matrix. k=4
F)
Delta area
::
I
11
-
1
=
relativt change in area under COF cutve
3
2
0.2
2
1
6
$
.
7
®
9
(G)
(H)
₡
1.00-
geneCluster
- A
B
val probability
0.75
C
☏ 0.50
Survival
0.25
p<0.
001
0.00
0
N
6
7
B
9
10
11
12
13
14
15
Time(years)
geneCluster
Number at risk
OD2
w
O
10-0
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Time(years)
(I)
geneCluster II
A
B
C
**
7.5
Gene expression
5.0
2.5
0.0
METTL3
METTL14
WTAP
VIRMA
ZC3H13
RBM15
RBM15B
YTHDC1
YTHDC2
YTHDF1
YTHDF2
YTHDF3
HNRNPC
FMR1
LRPPRC
HNRNPA2B1
IGFBP1
IGFBP2
IGFBP3
FTO
ALKBH5
This could help us to understand TME anti-tumor immune response clearly and develop efficient immunotherapy.
We integrated the m6A-related data from the ACC in the TCGA database and the GEO database, formed two clusters using unsupervised cluster analysis, and found that the two m6A clusters had significantly different OS and TME infiltration characteristics. Here, we identified two unique m6A methylation alteration patterns based on 23 m6A regulators. The TME cell invasion in these two patterns was
noticeably different. Cluster B was defined by the activation of adaptive immunity, which corresponds to an immune- inflamed phenotype, whereas cluster A was characterized by the suppression of immunity, which corresponds to an immune-desert phenotype. The immune-desert phenotypes could be considered non-inflamed tumors. The immune- inflamed phenotype, referred to as a hot tumor, is characterized by a significant immune cell infiltration in the TME. In particular, m6A cluster B had greater immune cell
(A)
(B)
(C)
geneCluster
PARBE
C
Activated. B.cellna
Activated CD4.T.cellna
Activated CDB.T.cellna
Activated dendritic.cellna
Eosinophilna
Common delta T.celina
Immature. B.celina
Immature.dendritic.celina
Natural killer. Tcellna
mCAscore
Mastcellna
Monocytena
Natural killer cellna
Neutrophilna
Regulatory. T.celina
T.folicularhelper. cellna
Type.1. Thewelcome
Type.17.T.helpercena
Type.2. Thelper.cellna
3.6e-08
MOSCna
2.4e-13
A
High
5.8e-14
Alive
meAscore
.
1
Activated.B.celina
Activated.CD4.T.celina
0.8
200
Activated.CDB.T.celina
Activated dendritic.celna
A
CDS6bright.natural klier.celna
0.6
CDSEdim.natural killer.celna
m6Ascore
Eosinophilna
Gamma.delta.T.cellna
.
0.4
Immature.8.celina
100
0.2
Immature.dendritic.celina
B
MDSCna
Macrophagena
0
Low
Mast celina
Monocytena
-02
Dead
Natural killer. T.celina
.
Natural kiler cellina
0
Neutrophilna
-0.4
Plasmacytoid dendritic.celna
B
Regulatory.T.celna
-0.6
c
Tfollicular.helper. celna
Type.1.T.helper. celna
Type.17.T.helper.celna Type.2.Thelper.celna
-0.8
.
A
B
C
m6Acluster
geneCluster
mfAscore
Fustat
-1
geneCluster
(D)
(E)
(F)
(G)
m6Acluster A B
100
Gender
Male
Female
100
3.3e-09
0.036
200
200
39%
75
54%
75
70%
Gender
*
78%
m6Ascore
Percent weight
Percent weight
Fustat
100
50
Female
m6Ascore
100
50
Alive
Male
Dead
61%
0
25
46%
0
25
30%
22%
0
0
A
B
m6Acluster
Low
High
Male
Female
Gender
Low
High
m6Ascore
m6Ascore
(H)
(I)
(J)
Patients with Male
(K)
Fustat
Dead
Alive
1.00
m6Ascore
High
Low
Patients with Female
m6Ascore
0.00067
Low
1.00-
m6Ascore
High 1
Low
200
Survival probability
0.75
High
Survival probability
1.00-
0,50
0.75
Survival probability
0.75
*
m6Ascore
0.50
100
0.50
0.25
p<0.001
0.25
p<0.001
0.25
p=0.003
0.00
0.00
0.00
0
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Time(years)
0
1
2
3
4
5
6
7
8
9
0
1
2
3
4
5
6 7 8 9 10 11 12 13 14 15
Time(years)
Time(years)
m6Ascore
Number at risk
m6Ascore
Low
ligh
56
44
28
21
12
9
5
3
3
3
2
1
1
m6Ascore
6
3
1
1
1
High Low
35 32 39 15 13 12 19 2 2 3 3 3 3 9 9 9
46
45
39
30
25
21
16
12
8
4
Dead
0
0.
Alive
3
0
High-
21
21
19
15
12
NO
NO
4
لى حـ
0
3
56 7 8 9 10 11 12 1
14
15
Fustat
0
1
2
3
1 2
4
5
13
4
5
6
7
8
9
10
11
12
13
14
15
Low
12
5
4
A
1
Time(years)
0
1
2
3
4
5
6
7
8
9
Time(years)
Time(years)
infiltration than m6A cluster A, including activated CD8 T cells, CD56 NK cells, monocytes, neutrophils, and type 17 T helper cells. TME immune cell infiltration was widespread in ACC, according to our hypothesis, immune pathway activity was affected by m6A-related regulators, resulting in a range of cell infiltration characteristics. Through death receptor- mediated or mechanism-independent apoptosis, CD56 NK cells might promote the death of tumor cells [25]. The development of early immunological responses, adaptive responses (IFN-c), and the control of NK cells (IL-10) may be significantly influenced by CD56 bright NK cells [26]. T
cell receptor-mediated cytotoxicity affects T cells and NK cells, moreover, the higher degree of T cell infiltration also improves protection against cancer [27]. Interestingly, the m6A cluster survival curve supports our hypothesis that patients with high levels of immune cell infiltration had a better prognosis. The immune cell infiltration of TME might prevent the emergence and progression of ACC, which offers evidence for immunotherapy utilization.
The prognostic variations across various m6A alteration patterns were strongly associated with immune-related biological pathways according to ssGSEA, GO, and KEGG
(A)
(B)
(C)
12
Altered in 38 (86.36%) of 44 samples.
m6aScore
Low
High
R =- 0.42, p = 0.00012
-
943
12.5
0.00054
0
12
Q
TP53
CTNNB1
27%
8
25%
MUC16
11%
10.0
Tumor Burden Mutation
MUC4
TIN
23%
Tumor Burden Mutation
18%
geneCluster
CNTNAP5
PKHD1
16%
14%
7.5
NF1
A
ASXL3
11%
HMCN1
Sa
4
B
PCDH15
GAR 9%
c
DST
CMYA5
9%
5.0
PRKARIA
7%
9%
FBN2
SVEPA
5%
9%
CCDC168
MEN1
0%
2.5
11%
0
ANK2
LRP1
7%
9%
m6Ascore
0.0
. Missense Mutation
Frame_Shift_Del
Frame_Shift_Ins
In_Frame_Del
m6Ascore
High
Low
High
0
100
Nonsense_Mutation ” Multi_Hit
200
Low
(D)
(E)
m6Ascore
(F)
Altered in 19 (54.29%) of 35 samples.
290
1.00
1.00-
- H-TMB+H-m6Ascore
- H-TMB+L-m&Ascore
- L-TMB+H-m&Ascore
0
8
H-TMB
- L-TMB+L-m6Ascore
0-
TP53
CTNNB1
3% 6%
Survival probability
0.75
L-TMB
0.75
Survival probability
MUC16
17%
MUC4
TIN
CNTNAP5
2 0%
0.50
0.50
PKHD1
NF1
3%
ASXL3
3%
3%
HMCN1
9%
PCDH15
3%
0.25
p<0.001
0.25
DST
CMYA5
3%
6%
p<0.001
PRKARTA
FBN2
6%
SVEP1
6%
3%
0.00
CCDC168
0%
0.00
MEN1
3%
0
1
2
3
4
5
6
7
8
9
10
11
12
ANK2
Time(years)
0
1
2
3
4
5
6
7
8
9
10
11
12
LRP1
5%
0%
Time(years)
m6Ascore
Number at risk
Number at risk
. Missense Mutation Frame_Shift_Del
Frame_Shift_Ins
m6Ascore
H-TMB-
· Multi_Hit
24
22
13
9
4
3
1
0
High
0
0
0
0
0
H-TMB+H-mBAscore
H-TMB+L-mBAscore
4
4
3
3
1
0
0 2 NOOO
0
0
0
ONOO
Low
L-TMB
O
55
53
45
35
26
21
15
11
8
7
4
2
2
-TMB+H-m6Ascore
2
18
10
29
0
1
2
3
4
5
6
7
8
9
10
11
12
L-TMB+L-m&Ascore
24
22
16
14
8
6
3
A
1
O
Time(years)
0
1
2
3
4
5
6
7
8
9
10
11
12
Time(years)
analyses. DEGs, which were regarded as m6A-related characteristic genes, clustering outcome was comparable to the m6A cluster. The establishment of three gene clusters based on the m6A signature gene and their strong association with the OS of ACC patients further confirmed the link between m6A alteration and the progression of ACC. Therefore, analyzing the m6A alteration pattern will deepen the understanding of immune cell infiltration features in TMEs. The m6A alteration pattern of a single ACC patient was assessed by m6Ascore, considering the individual heterogeneity of m6A alteration. The higher m6Ascore, correlated with an immune-inflamed phenotype, was present in the m6A altered model that was characterized by high abundance immune cell infiltration. The lower m6Ascore in the m6A altered mode was associated with immunological desert phenotype and low quantity immune cell infiltration. Further investigation revealed that ACC patients with high m6Ascores had longer OS, demonstrating the predictive value of the m6A alteration, and the inhibitory effect of high immune cell infiltration in TME. In summary, the m6Ascore model could be applied to predict the TME infiltration features and ACC patients survival accurately.
TMB is a detectable biomarker [28] that can be used to identify patients with good responses to immune checkpoint
inhibitor therapy [29]. The higher TMB has also been shown to correlate with clinical benefits from ICI therapy within multiple tumors, in our study, TMB is also a novel antigens generation that activates anti-tumor immunity [30- 32]. These immunotherapies have shown remarkable anti- tumor effects in prior clinical trials [27,33], and have received approval for a wide range of tumor types with a variety of indications, including metastatic castration- resistant prostate cancer [34]. Our findings revealed a substantial inverse relationship between m6Ascore and TMB, despite ACC patients with high TMB having a poorer OS, it is an immunoreactive tumor and a “hot tumor,” which has significant therapeutic implications since these patients might be more sensitive to immune medicine.
We also verified that there was a difference in NIK and PD-L1 expression between the two m6Ascore groups. Numerous costimulatory substances have an effect on NIK downstream proteins. According to a melanoma study, the loss of NIK resulted in both an increase in tumor burden and a reduction in tumor-infiltrating T cells, indicating that NIK is essential for both T cell survival and anti-tumor immunity. NIK-targeted therapy offers a perspective for investigating immunotherapy, as it identifies a novel regulatory component of T cell metabolism. The NIK- targeted treatment in combination with immune checkpoint
(A)
(B)
(C)
m6Ascore
Low
High
m6Ascore
Low
High
m6Ascore
Low
High
0.092
2.5-
0.06
0.002
2.5
·
·
·
·
·
·
·
4
PD-L1 expression
2.0
2.0
CTLA-4 expression
.
NIK expression
3
1.5
·
1.5
·
2
1.0
1.0
1
0.5
·
Low
High
Low
High
Low
High
m6Ascore
m6Ascore
m6Ascore
inhibitors could be considered for ACC patients who missed the best chance for surgery. There are limitations in the current study, and bench experiments will be conducted to support our findings.
As for the validation study, the Kaplan-Meier curve and Log-rank test of low and high m6Ascore patients in the combined cohort to validate the accuracy of our model. Although these validated data are not from the same source as the data we modeled, in multiple datasets, m6Ascore has a credible prognostic value, demonstrating the reliability of our conclusions.
Conclusions
In conclusion, TME immune cell infiltration features and m6A methylation alteration pattern could be assessed by the m6Ascore model. Furthermore, it was also applied to determine the TMB, OS, and m6A genotyping in ACC patients. More importantly, the effectiveness of novel immune checkpoint blockades (PD-1/L1 & CTLA4) based on the m6Ascore in the clinic was predicted, such as NIK- targeted treatment with immune checkpoint inhibitors, offering some new insights for ACC immunotherapy.
Acknowledgement: None.
Funding Statement: The authors extend their appreciation to the Researchers Supporting Project Number (RSPD2023R725) King Saud University, Riyadh, Saud Arabia.
Author Contributions: Wenhao Xu, Ayman Moawad Mahmoud, and Chen Li conceived and designed the study. Wenhao Xu, Haoming Li, and Yasir Hameed analyzed the bioinformatics data. Mostafa A. Abdel-Maksoud, Saeedah Musaed Almutairi, Ayman Mubarak, Mohammed Aufy, Wael Alturaiki, and Abdulaziz J. Alshalani joined the data analysis. All the authors contributed to the manuscript draft and approved the final version.
Availability of Data and Materials: All data used in this work can be acquired from the TCGA database (https://portal.gdc. cancer.gov/), GEO database (https://www.ncbi.nlm.nih.gov/ geo/).
Ethics Approval: None.
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.
References
1. Kerkhofs, T., Verhoeven, R., van der Zwan, J., Dieleman, J., Kerstens, M. et al. (2013). Adrenocortical carcinoma: A population-based study on incidence and survival in the Netherlands since 1993. European Journal of Cancer, 49(11), 2579-2586. https://doi.org/10.1016/j.ejca.2013.02.034
2. Icard, P., Goudet, P., Charpenay, C., Andreassian, B., Carnaille, B. et al. (2001). Adrenocortical carcinomas: Surgical trends and results of a 253-patient series from the french association of endocrine surgeons study group. World Journal of Surgery, 25(7), 891-897. https://doi.org/10.1007/s00268-001-0047-y
3. Schulick R. D., Brennan M. F. (1999). Long-term survival after complete resection and repeat resection in patients with adrenocortical carcinoma. Annals of Surgical Oncology, 6(8), 719-726. https://doi.org/10.1007/s10434-999-0719-7
4. Grubbs, E., Lee, J. J. C. (2009). Limited prognostic value of the 2004 international union against cancer staging classification for adrenocortical carcinoma: Proposal for a revised TNM classification. Cancer, 115(24), 5847-5848. https://doi.org/10. 1002/cncr.24693
5. Gonzalez, R., Tamm, E., Ng, C., Phan, A., Vassilopoulou-Sellin, R. et al. (2007). Response to mitotane predicts outcome in patients with recurrent adrenal cortical carcinoma. Surgery, 142(6), 867-875. https://doi.org/10.1016/j.surg.2007.09.006
6. Icard, P., Chapuis, Y., Andreassian, B., Bernard, A., Proye, C. J. S. (1992). Adrenocortical carcinoma in surgically treated patients: A retrospective study on 156 cases by the french association of endocrine surgery. Surgery, 112(6), 972-980.
7. Dubin D. T., Taylor R. H. (1975). The methylation state of poly A-containing messenger RNA from cultured hamster cells. Nucleic Acids Research, 2(10), 1653-1668. https://doi.org/10. 1093/nar/2.10.1653
8. Zhao B. S., Roundtree I. A., He C. (2017). Post-transcriptional gene regulation by mRNA modifications. Nature Reviews Molecular Cell Biology, 18(1), 31-42. https://doi.org/10.1038/ nrm.2016.132
9. Cai, X., Wang, X., Cao, C., Gao, Y., Zhang, S. et al. (2018). HBXIP-elevated methyltransferase METTL3 promotes the
progression of breast cancer via inhibiting tumor suppressor let- 7G. Cancer Letters, 415, 11-19. https://doi.org/10.1016/j.canlet. 2017.11.018
10. Ren, F., Zhao, Q., Zhao, M. H., Zhu, S. G., Liu, B. et al. (2021). Immune infiltration profiling in gastric cancer and their clinical implications. Cancer Science, 112(9), 3569-3584. https://doi.org/10.1111/cas.15057
11. Topalian, S., Hodi, F., Brahmer, J., Gettinger, S., Smith, D. et al. (2012). Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. The New England Journal of Medicine, 366(26), 2443-2454. https://doi.org/10.1056/NEJMoa1200690
12. Wang, Q., Hu, B., Hu, X., Kim, H., Squatrito, M. et al. (2018). Tumor evolution of glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. Cancer Cell, 32(1), 42-56.e6. https://doi. org/10.1016/j.ccell.2017.06.003
13. Li, H., Tong, J., Zhu, S., Batista, P., Duffy, E. et al. (2017). mA mRNA methylation controls T cell homeostasis by targeting the IL-7/STAT5/SOCS pathways. Nature, 548(7667), 338-342. https://doi.org/10.1038/nature23450
14. Lu, X., Gao, C., Liu, C., Zhuang, J., Su, P. et al. (2019). Identification of the key pathways and genes involved in HER2-positive breast cancer with brain metastasis. Pathology, Research and Practice, 215(8), 152475. https://doi.org/10.1016/j. prp.2019.152475
15. Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L. et al. (2020). m6A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. BioMed Central, 19(1), 53. https://doi.org/10.1186/ s12943-020-01170-0
16. Sotiriou, C., Wirapati, P., Loi, S., Harris, A., Fox, S. et al. (2006). Gene expression profiling in breast cancer: Understanding the molecular basis of histologic grade to improve prognosis. Journal of the National Cancer Institute, 98(4), 262-272. https://doi.org/10.1093/jnci/dj052
17. Zeng, D., Li, M., Zhou, R., Zhang, J., Sun, H. et al. (2019). Tumor microenvironment characterization in gastric cancer identifies prognostic and immunotherapeutically relevant gene signatures. Cancer Immunolgy Resarch, 7(5), 737-750. https:// doi.org/10.1158/2326-6066.Cir-18-0436
18. Chen, D., Mellman, I. J. N. (2017). Elements of cancer immunity and the cancer-immune set point. Nature, 541(7637), 321-330. https://doi.org/10.1038/nature21349
19. Liu, X., Qin, J., Gao, T., Li, C., Chen, X. et al. (2020). Analysis of METTL3 and METTL14 in hepatocellular carcinoma. Sedentary Life and Nutrition, 12(21), 21638-21659. https://doi.org/10. 18632/aging.103959
20. Tang, M., Lv, Y. (2021). The role of Nº-methyladenosine modified circular RNA in pathophysiological processes. International Journal of Biological Sciences, 17(9), 2262-2277. https://doi.org/10.7150/ijbs.60131
21. Qing, Y., Su, R., Chen, J. J. B. (2021). RNA modifications in hematopoietic malignancies: A new research frontier. Blood, 138(8), 637-648. https://doi.org/10.1182/blood.2019004263
22. Liu, J., Ren, D., Du, Z., Wang, H., Zhang, H. et al. (2018). mA demethylase FTO facilitates tumor progression in lung
squamous cell carcinoma by regulating MZF1 expression. Biochemical and Biophysical Research Communications, 502(4), 456-464. https://doi.org/10.1016/j.bbrc.2018.05.175
23. Bansal, H., Yihua, Q., Iyer, S., Ganapathy, S., Proia, D. et al. (2014). WTAP is a novel oncogenic protein in acute myeloid leukemia. Leukemia, 28(5), 1171-1174. https://doi.org/10.1038/ leu.2014.16
24. Dixit, D., Xie, Q., Rich, J., Zhao, J. C. (2017). Messenger RNA methylation regulates glioblastoma tumorigenesis. Cancer Cell, 31(4), 474-475. https://doi.org/10.1016/j.ccell.2017.03.010
25. Zwirner, N., Domaica, C., Fuertes, M. B. (2021). Regulatory functions of NK cells during infections and cancer. Journal of Leukocyte Biology, 109(1), 185-194. https://doi.org/10.1002/JLB. 3MR0820-685R
26. Poli, A., Michel, T., Thérésine, M., Andrès, E., Hentges, F. et al. (2010). CD56bright natural killer (NK) cells: An important NK cell subset. Immunology, 126(4), 458-465. https://doi.org/10. 1111/j.1365-2567.2008.03027.x
27. Dar, A., Patil, R., Chiplunkar, S. V. (2014). Insights into the relationship between toll like receptors and gamma delta T cell responses. Frontiers in Immunology, 5, 366. https://doi.org/10. 3389/fimmu.2014.00366
28. Hainesworth, J. (2021). Researchers strive to refine TMB. Cancer Discovery, 11(6), 1314. https://doi.org/10.1158/2159-8290. CD-ND2021-0107
29. Goodman, A., Kato, S., Bazhenova, L., Patel, S., Frampton, G. et al. (2017). Tumor mutational burden as an independent predictor of response to immunotherapy in diverse cancers. Molecular Cancer Therapeutics, 16(11), 2598-2608. https://doi. org/10.1158/1535-7163.MCT-17-0386
30. Grizzi, G., Caccese, M., Gkountakos, A., Carbognin, L., Tortora, G. et al. (2017). Putative predictors of efficacy for immune checkpoint inhibitors in non-small-cell lung cancer: Facing the complexity of the immune system. Expert Review of Molecular Diagnostics, 17(12), 1055-1069. https://doi.org/10.1080/ 14737159.2017.1393333
31. Gubin M., Artyomov M., Mardis E., Schreiber R. D. (2015). Tumor neoantigens: Building a framework for personalized cancer immunotherapy. The Journal of Clinical Investigation, 125(9), 3413-3421. https://doi.org/10.1172/JCI80008
32. Juergens, R., Zukotynski, K., Singnurkar, A., Snider, D., Valliant, J. et al. (2016). Imaging biomarkers in immunotherapy. Imaging Cancer, 1(2), e190031. https://doi.org/10.4137/BIC.S31805
33. Robert, C., Ribas, A., Wolchok, J., Hodi, F., Hamid, O. et al. (2014). Anti-programmed-death-receptor-1 treatment with pembrolizumab in ipilimumab-refractory advanced melanoma: A randomised dose-comparison cohort of a phase 1 trial. The Lancent, 384(9948), 1109-1117. https://doi.org/10.1016/ S0140-6736(14)60958-2
34. Graff, J., Liang, L., Kim, J., Stenzl, A. (2021). KEYNOTE-641: Phase III study of pembrolizumab plus enzalutamide for metastatic castration-resistant prostate cancer. Future Oncology, 17(23), 3017-3026. https://doi.org/10.2217/ fon-2020-1008
Supplementary Materials
Table S1: Spearman correlation analysis of the 23 m6A modification regulators
Table S2: Basic information of datasets included in this study for identifying distinct m6A methylation modification patterns Table S3: Estimating relative abundance of tumor microenvironment cells in ACC patients by the Single-Sample Gene-Set Enrichment Analysis (ssGSEA)
Table S4: The activation states of biological pathways by GSVA enrichment analysis
Table S5: Prognostic analysis of m6A phenotype-related genes using a univariate Cox regression model Table S6: Functional annotation for m6A phenotype-related genes (Gene Ontology-Biological process) Table S7: The changes of m6A clusters, survival time and m6Ascore