Original Article N7-methylguanosine regulatory genes well represented by METTL1 define vastly different prognostic, immune and therapy landscapes in adrenocortical carcinoma

Fangshi Xu1,2, Danrui Cai3, Shanshan Liu4, Kaini He5, Jing Chen4, Li Qu4, Tie Chong1, Xueyi Li4, Bincheng Ren4

1Department of Urology, Second Affiliated Hospital of Xi’an Jiaotong University, No. 157, West Five Road, Xi’an 710004, Shaanxi, China; 2Department of Urology, Shaanxi Provincial People’s Hospital, No. 256, Friendship West Road, Xi’an 710068, Shaanxi, China; 3Department of Ophthalmology, Second Affiliated Hospital of Xi’an Jiaotong University, No. 157, West Five Road, Xi’an 710004, Shaanxi, China; 4Department of Rheumatology and Immunology, Second Affiliated Hospital of Xi’an Jiaotong University, Xi’an, Shaanxi, China; 5Department of Gastroenterology, Second Affiliated Hospital of Xi’an Jiaotong University, Xi’an, Shaanxi, China

Received May 16, 2022; Accepted January 30, 2023; Epub February 15, 2023; Published February 28, 2023

Abstract: Although N7-methylguanosine (m7G) is one of the most frequent RNA modifications, it has received little attention. Adrenocortical carcinoma (ACC) is a highly malignant and easily metastatic tumor, eagerly needing for novel therapeutic strategy. Herein, a novel m7G risk signature (METTL1, NCBP1, NUDT1 and NUDT5) was con- structed using the Lasso regression analysis. It possessed highly prognostic value and could improve the predictive accuracy and clinical making-decision benefit of traditional prognostic model. Its prognostic value was also success- fully validated in GSE19750 cohort. Through CIBERSORT, ESTIMATE, ssGSEA and GSEA analyzes, high-m7G risk score was found to be closely associated with increased enrichment of glycolysis and suppression of anti-cancer immune response. Therapeutic correlation of m7G risk signature was also investigated using tumor mutation bur- den, the expressions of immune checkpoints, TIDE score, IMvigor 210 cohort and TCGA cohort. m7G risk score was a potential biomarker for predicting the efficacy of ICBs and mitotane. Furthermore, we explored the biofunctions of METTL1 in ACC cells through a series of experimentations. Overexpression of METTL1 stimulated the proliferation, migration and invasion of H295R and SW13 cells. Immunofluorescence assays revealed that the infiltrating levels of CD8+ T cells was lower and that of macrophages was higher in clinical ACC samples with high METTL1 expres- sion compared to that in low expression ones. Silencing METTL1 could significantly inhibited tumor growth in mouse xenograft model. Western blot assays showed that METTL1 positively regulated the expression of glycolysis rate- limiting enzyme HK1. Finally, miR-885-5p and CEBPB were predicted as the upstream regulators of METTL1 through data mining of the public databases. In conclusions, m7G regulatory genes well represented by METTL1 profoundly affected the prognosis, tumor immune, therapeutic outcomes, and malignant progression of ACC.

Keywords: N7-methylguanosine, adrenocortical carcinoma, risk signature, prognosis, METTL1

Introduction

Adrenocortical carcinoma (ACC) is a highly malignant and aggressive urologic cancer with a rare incidence of 2.0 per million [1]. Due to its atypical incipient symptoms, patients common- ly present advanced or metastatic disease at the time of initial diagnosis, with a five-year overall survival rate (OSR) < 15% [2]. Radical adrenalectomy is the only option to achieve ACC cure. However, cases suitable for surgery account for only approximately 30% of all ACC patients [3]. Apart from radical resection, other

current therapeutic approaches have not achieved satisfactory effects. As the mainstay adjuvant therapy of ACC, mitotane combined with EDP (etoposide, doxorubicin, and cisplat- in), can produce an objective response rate (ORR) of less than 30%, and progression-free survival (PFS) of patients who receive this treat- ment is only 5.6 months [4]. It is thus necessary to develop novel therapeutic strategies so as to improve prognostic evaluation system for ACC.

RNA epigenetic modulation is a current topic in oncology. One prominent example is N6-methy-

Table 1. Three core m7G gene sets from MSigDB database
NamesGene countsDescription
GOMF m7G 5-PPPN Diphosphatase Activity12Catalysis of the reaction: 7-methylguanosine 5'-triphospho-5'-polynucleotide + H20 = 7-methylguanosine 5'-phosphate + polynucleotide
GOMF RNA 7-Methylguanosine Cap Binding13Binding to a 7-methylguanosine group added cotranscriptionally to the 5' end of RNA molecules tran- scribed by polymerase II
GOMF RNA Cap Binding20Binding to a 7-methylguanosine (m7G) group or derivative located at the 5' end of an RNA molecule

ladenosine (m6A) modification, which is closely involved in the prognosis, immune response, and development of ACC [5-7]. N7-methy- lguanosine (m7G) is a further prevalent pattern of RNA modification, however, its roles in can- cer are so far unclear. m7G refers to the gua- nine methylation at the 5’-cap of RNA, com- monly occurring at position 46 (G46) in the vari- able region of the tRNA loop [8]. The functional complex consisting of methyltransferase 1 (METTL1) and the WD repeat domain 4 (WDR4) is responsible for this guanosine methylation process [9]. RNA exhibits higher stability after m7G modification, which has attracted research interest in oncology. METTL1/WDR4-mediated m7G tRNA modification promotes the progres- sion of lung cancer [10]. m7G modification has been a focus of cancer research, however, the association of m7G regulator genes and cancer prognosis, cancer treatment, and the tumor immune microenvironment (TIM) are so far unclear.

In view of above context, we constructed a novel m7G risk signature for ACC clinical assessment through Lasso regression analy- sis. Its prognostic potential, immune effects, metabolic impacts, mutation features, and therapeutic correlations were comprehensively investigated. More importantly, we confirmed the oncogenic abilities of METTL1, the most critical regulator in m7G modification, during ACC progression for the first time through in vitro experiments. Our findings provide new insights regarding ACC treatment and assess- ment options.

Materials and methods

Data source

Clinical and transcriptomic data were retrieved from the TCGA (https://portal.gdc.cancer.gov/)

and GEO (https://www.ncbi.nlm.nih.gov/geo/) public databases. No normal samples were available in the TCGA-ACC project, thus we used 128 normal adrenal tissue samples from the GTEx database (https://xenabrowser.net/ datapages/) to screen differentially express- ed genes (DEGs). All transcriptome data was standardized through log2 (FPKM+1) transfor- mation. The clinical characteristics of the TCGA and GSE19750 cohorts are shown in Supplementary Table 1.

We reviewed studies on m7G modification and three pivotal gene sets in the Molecular Signatures Database (MSigDB) to establish an m7G-related gene set, which comprised 34 m7G regulators (Supplementary Table 2). Three MSigDB gene sets included ‘GOMF m7G 5-PPPN Diphosphatase Activity’, ‘GOMF RNA 7-Methylguanosine Cap Binding’ and ‘GOMF RNA Cap Binding’. Respective detailed descrip- tions are shown in Table 1. To further confirm the reliability of our m7G gene set, we con- structed its protein-protein interaction (PPI) network and conducted the corresponding bio- logical function analyses using the Metascape online tool (http://metascape.org/) [11].

m7G-related DEGs were identified using the ‘Limma’ package in R software (version 4.1.2). The following screening criteria thresholds were used: adjusted p-value < 0.05 and the ab- solute value of log2 FC >1 (2-fold difference in gene expression). Next, we identified prognos- tic m7G genes through Cox univariate regres- sion analysis. The intersection between DEGs and prognostic genes was obtained through a Venn diagram. Finally, we established a novel

m7G risk signature of ACC using Lasso regres- sion analysis.

Evaluation of the prognostic value

The optimal cutoff value of the m7G risk score was calculated using the Cutoff Finder online tool (http://molpath.charite.de/cutoff) [12]. According to this cutoff value, 79 ACC sampl- es were assigned to high- and low-m7G risk groups. Then, survival differences between the risk groups were determined through Kaplan- Meier analyses. Cox univariate and multivariate analyses were used to identify the independ- ent prognostic factors. The predictive accuracy of the m7G risk signature was evaluated through a receiver operating characteristic curve (ROC). Decision curve analysis (DCA) was applied to assess whether the m7G risk score could improve the traditional prognostic model based on clinical stage. The clinical subgroup analyses were conducted to ensure the appli- cable scope of the m7G model in prognostic analyses. Due to insufficient samples in N1 stages (n = 10), we did not perform survival dif- ference analyses in this subgroup. We utilized a nomogram comprising TNM-staging and m7G risk scores to predict the overall survival rate of individuals at 2, 3, and 5 years. Its prognostic accuracy was assessed through calibration curves. Further, the prognostic value of the m7G risk signature was validated in the GSE19750 dataset.

Immune analyses

The infiltration levels of 22 immune cell sub- types in each ACC sample were calculated using the CIBERSORT algorithm [13]. The activi- ties of 10 immune-related pathways were quan- tified using single-sample gene set enrichment analysis (ssGSEA) [14]. The R package ‘Limma’ was applied to determine differences in infiltra- tion levels of immune cells and the activities of immune-related pathways between different m7G risk groups. The ESTIMATE algorithm was employed to compare differences in stromal, immune, and ESTIMATE scores between high- and low-risk groups [15]. The corresponding tumor purity of each ACC sample was quanti- fied through the same algorithm. The TIMER database offers a comprehensive resource for systematical analysis of immune infiltrates across diverse cancer types (https://cistrome. shinyapps.io/timer/) [16]. The correlations

between the somatic copy-number alterations (SCNAs) of m7G signature genes and the abundance of six core immune cells were ana- lyzed using the ‘SCNA’ module in the TIMER database.

Gene set enrichment analysis (GSEA)

GSEA was utilized to assess the impacts of m7G risk scores on multiple metabolic pro- cesses, including glycolysis, nucleotide metab- olism, amino acid (AA) metabolism, and fatty acid (FA) metabolism. Analytical gene sets were obtained from the MSigDB database, and their basic information is presented in Supplementary Table 3. Phenotype labels were set as high-m7G risk samples versus low-m7G risk ones. The number of permutations was 1,000, and gene symbols were not collapsed.

Calculation of the tumor mutation burden (TMB) and mutational analyses

The TMB of each ACC sample was calculated as the total mutation frequency divided by 38. The corresponding somatic mutation data were obtained from the TCGA database. The ‘Data Type’ and ‘Workflow Type’ were ‘Masked somatic mutation’ and ‘VarScan’, respectively. Mutational data were visualized using the R package ‘maftools’. The cBioPortal database (http://cbioportal.org) was used to acquire information on somatic mutation frequency and patterns of m7G signature genes across two ACC projects (n = 184 samples).

Therapeutic correlation analyses

The TCGA-ACC cohort was used to compare dif- ferences in m7G risk score between responding and non-responding patients treated with radi- ation and mitotane therapy. Then, we explored potential linkages between the efficacy of ICBs and the m7G risk score from four perspectives, including TMB, tumor immune dysfunction and exclusion (TIDE) scoring, immune checkpoints (ICs) expression, and the IMvigor 210 cohort. TMB is considered a promising biomarker for predicting the efficacy of immune checkpoint blockades (ICBs) [17, 18], thus TMB differences between high- and low-risk groups were deter- mined. The TIDE scoring system is a crucial method to predict patient responses to anti- PD-1/L1 and anti-CTLA4 treatments based on the estimation of T cell dysfunction and tumor

immune evasion [19]. Using an online tool (http://tide.dfci.harvard.edu/login/), the TIDE score of each ACC patient was calculated, and differences in TIDE scores between the high- and low-m7G risk groups were determined. Expression levels of ICs can reflect the poten- tial to benefit from ICBs [20], thus correlations of expression of six ICs and m7G risk scores were tested. We then used the IMvigor210 cohort that reported the therapeutic outcomes of PD-1 inhibitor atezolizumab and the corre- sponding transcriptomic data [21] to confirm differences in m7G risk score between patients responding and not responding to therapy.

Analysis in upstream regulatory mechanism of METTL1

Three miRNA databases were employed to predict potential miRNA responsible for nega- tively regulating METTL1, including miRDB (http://www.mirdb.org/) [22], TargetScanHum (Ver 8.0, http://www.targetscan.org/vert_80/) [23] and ENCORI (https://starbase.sysu.edu. cn/) [24]. The minimum free energy of predict- ed miRNAs was calculated using RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/rna- hybrid/) online tool [25]. The binding site between miRNA and METTL1 was predicted using TargetScanHum database. Using Gene- Cards (https://www.genecards.org/) [26], ALG- GEN (https://alggen.lsi.upc.es/cgi-bin/promo_ v3/) [27] and hTFtarget (http://bioinfo.life.hust. edu.cn/hTFtarget/) [28] databases, we also investigated the potential regulatory transcrip- tion factors (TFs) of METTL1. The promoter sequences of METTL1 was obtained from UCSC genome database (http://genome.ucsc. edu/) [29]. The motif sequence of candidate TF and the prediction of binding site were derived from JASPAR database (https://jaspar.genereg. net/) [30].

Cell culture and transfection

Two adrenocortical cancer cell lines (H295R and SW13) were used for in vitro experiments. All cells were purchased from Procell Life Science & Technology Company (Wuhan, China). H295R cells were cultured in Dulbe- cco’s Modified Eagle Medium (DMEM) with 10% fetal bovine serum (FBS), termed DMEM/F12, ITS-G (an insulin, transferrin, and selenium solution) and 1% penicillin-streptomycin (P/S). SW13 cells were cultured in DMEM medium

containing 10% FBS and 1% P/S. sh-METTL1 and amplification plasmids (OE-METTL1) were purchased from HanHeng Biotechnology (Shanghai, China). Their respective sequences are shown in Supplementary Table 4. The cells were transfected using Lentiviruses (Hanheng Biotechnology, Shanghai, China).

Clinical samples and RT-qPCR

To confirm ectopic expression of METTL1 in ACC, we collected 10 pairs of ACC and adjacent normal tissues from the Department of Urology, Second Affiliated Hospital of Xi’an Jiaotong University to conduct RT-qPCRs. All patients provided written informed consent. The study protocol was approved by the Ethics Commit- tees of the Second Affiliated Hospital of Xi’an Jiaotong University.

Total RNA was extracted using TRIzol Reagent (TakaRa Bio, Shiga, Japan). RNA purity was measured based on the A260/A280 ratio us- ing a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Reverse transcription was performed using a PrimeScript RT Reagent Kit (TaKaRa Bio). Amplification was traced using SYBR-Green PCR Reagent (TaKaRa Bio) and an ABI Prism 7900 sequence detection system (Applied Biosystems, Foster City, CA, USA). GAPDH was used as an internal reference. The relative gene expression was calculated according to the 2-AACT method. Primer sequences are shown in Supplementary Table 5.

Western blot

The associations of METTL1 with four gly- colysis rate-limiting enzymes (PKM, PFKFB3, HK1 and HK2) were analyzed using Western blot. The experimental procedures were per- formed similar to previous study [31]. Briefly, transfected cells were lysed by RIPA buffer (Beyotime, China). After centrifugation, the supernatant was collected. Protein concentra- tion was measured by BCA kit (Phygene Life Sciences Company, Fuzhou, China). Sample proteins were separated by 10% SDS-PAGE. After electrophoresis, protein samples were transferred to PVDF membranes (BestBio, Shanghai, China). The PVDF membranes were blocked by 5% skim milk at 37℃ for 2 h. After washing by TBST buffer (BIOSIC, Nanjing, China) for three times, the membranes were incubat-

ed with the primary antibody (Abcam, UK) over- night at 4℃ and were incubated with the sec- ondary antibody (Abcam, UK) for 1 h at room temperature, respectively. Protein blots were exposure using ECL reagent (Abcam, UK) and detected by BioRad imaging system (BioRad, USA).

Colony formation assay

ACC cells at a density of 5 x 103 cells/well were seeded into six-well plates. When colonies were visible, they were washed using PBS, fixed with 4% paraformaldehyde, and were stained using Giemsa. Colonies were counted using a micro- scope at 20-fold magnification, with five ran- dom fields.

Transwell migration and invasion assays

For these assays, 24-well transwell chambers (Corning, NY, USA) were used. The experiment was conducted as described previously [5]. For transwell migration assays, DMEM/F12 or DMEM medium containing 0.1% FBS was added to the upper chambers, and medium containing 10% FBS was added to the lower chambers. After incubation for 24 h, migrated cells that adhered to the lower surface of the membrane were fixed by paraformaldehyde for 20 min and were stained with 0.1% crystal violet for 20 min. Cells in five random visual fields were counted at 20-fold magnification. When conducting the invasion assays, the upper chambers were precoated with Matrigel (Corning).

Immunofluorescence

We used 4-mm tissue sections of ACC clinical tissues for immunofluorescence assays as described previously [32]. Through immunoflu- orescence staining, the nucleus, CD8/CD163, and METTL1 were stained blue, red, and green respectively. The slides were analyzed using an automatic fluorescent microscope with a 40 x objective lens (Olympus BX53, Olympus, Tokyo, Japan).

Xenograft assay

We used six-weeks-old female BALB/c nude mice to conduct the tumor xenograft experi- ments. H295R cells that were stably transfect- ed with sh-METTL1 and sh-vector were injected

subcutaneously into the flanks of each mouse. The injection concentration and volume were 5 x 106 cells/mL and 100 uL, respectively. The tumor volume was calculated as 0.5 x tumor length x (tumor width)2. Tumor length and width were measured using a Vernier caliper every three days. After two weeks, all mice were killed, and xenograft tumors were collected. mRNA levels of METTL1 and P53 in xenograft tumors were evaluated by qPCR. This experi- ment was approved by the ethics committee of the Second Affiliated Hospital of Xi’an Jiaotong University.

Statistical analyses

All statistical analyses were performed using R software (version 4.1.2) and GraphPad Prism (version 8.0). Unpaired t-tests were used to test differences in continuous variables be- tween multiple experimental groups. Kolmo- gorov-Smirnov tests were used to assess the relationships between m7G risk scores and the clinicopathological characteristics of ACC. Survival analyses were based on the Kaplan- Meier method. Cell experiments were con- ducted using three independent replicates. Statistical significance is reported at P < 0.05.

Results

m7G is an important RNA modification

A flowchart of this study is shown in Figure 1. We established a reliable m7G-related gene set, by which a novel m7G risk signature was established. We then assessed its various roles in ACC clinical assessment and treat- ment. The main mechanism of m7G process is visualized in Figure 2A. m7G is most frequently located at position 46 in the tRNA variable region, termed G46 [8]. This methylation pro- cess is driven and catalyzed by the m7G func- tional complex that consists of two subunits, namely METTL1 and WDR4 [33]. The former exhibits methyltransferase activity, while the latter provides the molecular scaffolds for methylation reaction [8]. m7G can ultimately result in improving the stability of various modi- fied RNAs, including tRNA, mRNA, rRNA, and miRNA, which is strikingly different from m6A modification [9]. Further, m7G profoundly affects cancer progression, immune response, and drug resistance through modifying the

m7G and METTL1 in ACC

Figure 1. Flow chart of the present study. m7G, N7-methylguanosine; DEGs, differentially expressed genes; PCA, principal component analysis; ROC, receiver operating characteristic curve; DCA, decision curve analysis; SCNA, somatic copy number alteration; TMB, tumor mutation burden.

m7G related genes(n=34)

Validation cohort GSE19750 (n=22)

DEGs (n=16)

Prognositc genes (n=19)

Lasso Regression Analysis

m7G risk signature

Prognostic value

Immune effects

Metabolic reprogramming

Mutational information

Therapeutic correlation

PCA

Immune Cell abundance

Glycolysis

Mutation frequency

ICB

Survival difference

Amino acid metabolism

Mitotane

ESTIMATE score

Mutation pattern

Radiotherapy

ROC and DCA

Tumor purity

fatty acid metabolism

TMB

Independent factor

Immune pathway

Nucleotide metabolism

Mutational genes

SCNA

Clinical subgroups

mIHC

Expression

qPCR

Nomogram

Immune correlation

Immunohistochemistry

m7G risk signature

METTLI

Proliferation

Colony formation

Migration

Transwell assay

Invasion

Transwell assay

expressive status of pivotal regulatory genes [9, 10, 34].

Based on the regulatory mechanisms of m7G, we identified 34 core m7G-related genes from the MSigDB database. A PPI network of these m7G genes is shown in Figure 2B. Next, the hub module in m7G PPI network was identified (Figure 2C), in which METTL1 and WDR4 were included. Through biological function analyses, these genes were shown to be closely involved in tRNA methylation, RNA decapping, and the regulation of translation, which confirmed the reliability of our m7G gene set (Figure 2D).

A novel m7G risk signature for ACC assess- ment

Nearly half of m7G regulatory genes (16/34, 47.06%) were differentially expressed in ACC samples, compared to normal samples (Figure

3A). Up to 55.9% of the m7G genes were able to affect the prognosis of ACC (Figure 3B), and most of them were associated with unfavorable survival outcomes. Eight intersection genes were used for the Lasso regression analysis (Figure 3C), and a novel m7G risk signature was constructed as follows (Figure 3D-F): m7G risk score = 0.496 * (METTL1 relative expres- sion) + 0.714 * (NCBP1 relative expression) + 0.863 * (NUDT1 relative expression) + 0.576 * (NUDT5 relative expression). According to the optimal cutoff value of the m7G risk score (8.27), ACC patients in the TCGA cohort were assigned to high- and low-risk groups. PCA result indicated that the m7G risk score explained approximately 70% of the prognostic variance, confirming the capacity of our m7G model to stratify prognostic risks (Figure 3G). Further, patients with a high m7G risk were more likely to be in the late clinical, M, and T stages (Figure 3H).

Figure 2. m7G modification can improve the stability of multiple RNAs. A. Main mechanism of m7G modification. B. PPI network of 34 m7G regulatory genes. C. The hub module in the m7G PPI network. D. Biological function analyses of 34 m7G regulatory genes. PPI, protein-protein interaction.

A

B

tRNA

5’

3’

D-arm

T-arm

G46

CH3

m7G

METTL1

H2N

WDR4

Improve the stability

tRNA

RNA

mRNA

rRNA

miRNA

Cancer progression

Immune response

Drug resistance

+

C

D

NUDT1

CYFIP1

GO:0034655: nucleobase-containing compound catabolic process

GO:0006417: regulation of translation

NCBP2

R-HSA-8953854: Metabolism of RNA

METTL1

NCBP1

GO:0110154: RNA decapping

GO:0051030: snRNA transport

R-HSA-1169408: ISG15 antiviral mechanism

GO:0006446: regulation of translational initiation

EIF4E

WDR4

EIF4E

R-HSA-2393930: Phosphate bond hydrolysis by NUDT proteins

GO:0030488: tRNA methylation

hsa01521: EGFR tyrosine kinase inhibitor resistance

GO:1901655: cellular response to ketone

AGO2

EIF4E3

0.0

2.5

5.0

7.5

10.0

12.5

15.0

17.5

20.0

-log10(P)

m7G risk signature presents considerable prognostic value

The risk plots of the m7G signature are shown in Supplementary Figure 1. The proportion of death events in the high-risk group were sub- stantially higher than that in the low-risk group. Similarly, high m7G risk was associated with poor survival outcomes (HR = 12.78, P < 0.001; Figure 4A). With regard to prediction accuracy, the m7G risk score was the best indi- cator (AUC = 0.876) for prognostic assess- ment, compared to other traditional clinico- pathological characteristics of ACC (Figure 4B). Further, the m7G risk signature had the highest predictive accuracy for the three-year OSR of ACC patients (AUC = 0.953; Figure 4D). More importantly, applying m7G risk score to the prognostic model based on clinical stage great- ly increased its net benefit when making clini- cal decisions (Figure 4C). Further, combining the clinical stage and the m7G risk score also

improved previous predictive accuracy (AUC = 0.891; Figure 4E). Although clinical stage, T stage, M stage, and m7G risk score were asso- ciated with ACC prognosis (Figure 4F), only the m7G risk score was an independent progno- stic factor of ACC (HR = 4.103; Figure 4G). To determine the applicable scope of m7G risk sig- nature, we conducted clinical subgroup analy- ses. The results showed that m7G risk signa- ture could distinguish the survival differences of patients with each stage of ACC disease (Figure 4H-M). For the sake of clinical practice, we constructed a nomogram consisting of TNM-staging and m7G risk score to predict the 2, 3, 5-year survival rates of ACC patients (Figure 5A). The calibration plots indicated that the predicted probabilities were close to the actual survival rates (Figure 5B-D). Taken together, these results confirmed that the m7G risk signature is highly promising for prognostic assessment of ACC.

m7G and METTL1 in ACC

A

1

LSM1

0

5

10

15

20

C

m7G DEGs

Prognostic genes

D

E

Hazard ratio

77777

5

5 555

555

5554

4 4 4 3

33

20

332

7

5

5

4

Partial Likelihood Deviance

1.0

9.5

0.8

Coefficients

8

8

11

9.0

0.6

0.4-

8.5

0.2

0.0

8.0

-5

-4

-3

-2

-5

-4

-3

-2

Log(2)

Log(2)

Type EIF4E4Type NBpvalueHazard ratio
2TMETTL10.0042.606(1.361-4.991)
EIF4A1WDR4<0.0016.176(2.918-13.074)
NUDT4B0AGO2<0.0016.810(3.072-15.094)
CYFIP10.0372.057(1.045-4.053)
NUDT1-2DCPS<0.0015.898(2.613-13.313)
NUDT16EIF3D0.0251.863(1.081-3.211)
IFITS-4EIF4A1<0.0017.541(2.522-22.554)
EIF4E2<0.0014.559(2.514-8.265)
NCBP1LSM10.0132.493(1.208-5.144)
NUDTSNCBP1<0.0014.381(2.255-8.511)
METTL1NCBP20.0392.205(1.042-4.667)
NCBP2L0.0010.124(0.035-0.440)
NSUN2NCBP30.0273.648(1.161-11.462)
GEMIN5DCP2<0.0013.342(1.648-6.780)
NCBP2NUDT1<0.0013.532(2.128-5.864)
NUDT30.0093.510(1.373-8.976)
LARP1NUDT40.0400.558(0.320-0.972)
NUDT16L1NUDT5<0.0013.862(2.037-7.323)
EIF3DNUDT7<0.0010.402(0.254-0.636)
Figure 3. A novel m7G risk signature of ACC. A. Heatmap of m7G DEGs. B. Identification of prognostic genes through Cox univariate analysis. C. Overlap of m7G DEGs and m7G prognostic genes. D, E. Lasso regression analysis. F. Coefficients of 4 m7G signature genes. G. PCA results of m7G risk signature. H. Relationships between clinicopathological characteristics of ACC and m7G risk levels.

F

H

1.0

0.863

0.8

Coefficients

0.714

m7G Risk

Status

Stage

M

N

T

0.6

0.576

0.496

00000 00000

0.4

High (n = 28)

0.2

0.0

METTL1 NCBP1 NUDT1 NUDT5

G

3

2

PC2 (22.9%)

Low (n = 49)

1

Group

0

S Low-Risk

C High-Risk

- 1

p = 3.3e-08 p = 0.00075 p = 0.00028 p =0.28 p=

-2

p = 0.00055

-3

-4

-2

0

2

Alive Dead

II

IV

MO

M1

NO

N1

T1

T2

T3

T4

PC1 (46.3%)

A

B

TCGA-ACC cohort

C

m7G risk score

0.06

1.0

1.0

Model A

Low

Survival probability

High

Improved model A

0.8

0.8

0.04

All positive

0.6

Sensitivity (TPR)

Net Benefit

All negative

0.6

0.02

0.4

0.2

0.4

Gender (AUC = 0.506)

0.00

HR = 12.78 (4.99-32.77)

0.0

P < 0.001

Clinical stage (AUC = 0.787)

T (AUC =0.791)

0

2.5

5

7.5

10

12.5

0.2

M (AUC = 0.701)

-0.02

Overall survival time (Year)

N (AUC = 0.557)

Low

49

35

20

m7G risk score (AUC = 0.876)

8

4

2

0.0

-0.04

High

30

12

2

0

0

0.0

0.2

0.4

0.6

0.8

1.0

1-Specificity (FPR)

0.00

0.25

0.50

0.75

1.00

Threshold Probability

D

E

F

1.0

m7G risk score

Clinical stage+m7G risk score

1.0

0.8

0.8

Sensitivity (TPR)

Sensitivity (TPR)

0.6

0.6

0.4

0.4

0.2

1-Year (AUC = 0.855)

0.2

3-Year (AUC = 0.953)

model

AUC: 0.891

0.0

5-Year (AUC = 0.847)

0.0

Cl: 0.813-0.970

0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

1-Specificity (FPR)

1-Specificity (FPR)

0

5

10

G

H

Clinical stage I-II

Clinical stage III-IV

1.0

Low

1.0

Low

Survival probability

0.8

High

Survival probability

0.8

High

0.6

0.6

0.4

0.4

0.2

HR = 15.11 (2.88-79.28)

0.2

HR = 8.73 (1.98-38.55)

0.0

P = 0.001

0.0

P = 0.004

0

2.5

5

7.5

10

12.5

0

2.5

5

7.5

10

Overall survival time

(Year)

Overall survival time (Year)

Low

36

25

17

7

3

2

Low

12

9

3

1

1

High

10

4

2

0

0

High

19

0

0

0

0

2

4

6

J

T 1-2

K

T 3-4

L

MO

M

NO

1.0

Low

1.0

1

Low

1.0

Low

1.0

Low

Survival probability

0.8

High

Survival probability

0.8

High

Survival probability

0.8

High

Survival probability

0.8

High

0.6

0.6

0.6

0.6

0.4

0.4

0.4

0.4

0.2

HR = 19.53 (3.87-98.50)

0.2

0.2

HR = 5.05 (1.14-22.29)

HR = 12.26 (3.76-40.00)

0.2

HR = 11.76 (4.19-33.01)

0.0

P < 0.001

0.0

P = 0.032

0.0

P < 0.001

0.0

P < 0.001

0

2.5

5

7.5

10

12.5

0

1

2

3

4

5

6

0

2.5

5

7.5

10

12.5

0

2.5

5

7.5

10

12.5

Overall survival time

(Year)

Overall survival time

(Year)

Overall survival time

(Year)

Overall survival time (Year)

Low

40

29

19

8

4

2

Low

8

8

5

4

3

1

0

Low

45

31

20

8

4

2

Low

44

30

18

7

3

2

High

11

5

2

0

0

High

18

16

6

3

0

0

High

17

7

2

0

0

High 24

10

2

0

0

CharacteristicsHR (95% CI)P value
Gender1.043(0.473-2.299)0.917
Clinical stage2.902(1.844-4.569)<0.001
T3.363(2.098-5.393)<0.001
M6.038(2.664-13.683)<0.001
N2.058(0.774-5.472)0.148
m7G risk score4.410(2.726-7.134)<0.001
CharacteristicsHR (95% CI) P value
Clinical stage1.377(0.328-5.776)1 F 0.662
T2.063(0.747-5.693) N 0.162 I
M0.290(0.057-1.485) I 0.137
m7G risk score4.103(2.253-7.475) 1 <0.001

Figure 4. Prognostic value of the m7G risk signature. A. Overall survival difference between high- and low-m7G risk groups. B. Accuracy of m7G risk score and clinical characteristics of ACC for predicting OSR. C. The DCA results. Model A (blue line) represents the prognostic model based on clinical stage. Improved model A (red line) represents model A with m7G risk score added. D. Time-dependent accuracy of m7G risk score for predicting OSR. E. Predictive accuracy of the combination of clinical stage and m7G risk score. F, G. Identification of ACC-independent prognostic factors through Cox univariate (blue) and multivariate (red) analyses. H-M. Clinical subgroup analyses. OSR, overall survival rate.

The prognostic value of m7G risk signature in a validation cohort

Taking a further step, we tested the prognostic value of m7G risk signature using GSE19750 cohort. As expected, high m7G risk scores were

associated with an unfavorable prognosis in the GSE19750 cohort (Figure 5E). The predic- tive accuracy of the m7G risk signature in the GSE19750 cohort was approximately 0.80, which was slightly lower than that in the TCGA cohort (Figure 5F). In addition, the m7G risk

m7G and METTL1 in ACC

A

Points

0

20

40

60

80

100

B

1.00

C

x

1.0

Actual 1-Year OS(proportion)

0.95

0.9

T

T2

T3

Actual 2-Year OS(proportion)

0.8

T1

T4

0.90

M1

0.7

M

0.85

MO

0.6

N

N1

0.80

0.5

NO

0.75

0

m7G Risk

High

0.88

0.90

0.92

0.94

0.96

0.98

1.00

0.6

0.7

0.8

0.9

1.0

Low

Nomogram-Predicted Probability of 1-Year OS

Nomogram-Predicted Probability of 2-Year OS

Total Points

D

9

14

0

40

80

120

160

200

240

280

0.9

Linear Predictor

Actual 3-Year OS(proportion)

0.8

13

-3

-2

-1

0

1

2

3

m7G risk score

0.7

2-year Survival Probability

12

0.8

0.6 0.4 0.2

0.6

3-year Survival Probability

0.5

11

Spearman

0.8

0.6 0.4 0.2

5

r = 0.062

P = 0.770

5-year Survival Probability

:

10

0.8

0.6 0.4 0.2

0.5

0.6

0.7

0.8

0.9

0

500

1000

1500

2000

Nomogram-Predicted Probability of 3-Year OS

Tumor weight (g)

E

GSE19750 cohort

F

GSE19750 cohort

G

GSE19750 cohort

H

GSE19750 cohort

1.0

Low

1.0

14

*

14

ns

Survival probability

0.8

High

0.8

0.6

Sensitivity (TPR)

m7G risk score

13

m7G risk score

13

0.4

0.6

12

12

0.2

0.4

-

HR = 4.51 (1.42-14.28)

0.0

P = 0.011

11

0

5

10

15

0.2

11

Overall survival time (Year)

1-Year (AUC = 0.682)

3-Year (AUC = 0.890)

Low

11

7

4

2

0

0.0

5-Year (AUC = 0.791)

10

10

High

11

2

0

0

0

0.0

0.2

0.4

0.6

0.8

1.0

1-Specificity (FPR)

Stage I-II

Stage III-IV

Tumor size<5cm Tumor size≥5cm

J

K

Non-Secretory function

F

ns

I

4

Grade 3-4

ns

Secretory function

Grade 1-2

10

11

12

13

14

10

11

12

13

14

m7G risk score

m7G risk score

Figure 5. Validation of the prognostic value of m7G risk signature. A. Nomogram consisting of TNM-staging and m7G risk levels. B-D. Calibration plots for evaluating the predicting accuracy of m7G nomogram. E. Survival difference between high- and low-m7G risk groups in the GSE19750 cohort. F. Time-dependent predictive accuracy of m7G risk scores in the GSE19750 cohort. G. Difference in m7G risk scores between ACC patients at clinical stages I-II and III-IV. H. Difference in m7G risk scores between tumor sizes. I. Correlation of tumor size and m7G risk score. J. Relationship between ACC secretory function and m7G risk score. K. Relationship between histological grade and m7G risk score. * P < 0.05; NS, not statistically significant.

scores in patients at clinical stage III-IV were markedly higher than those of patients at clini- cal stage I-II (Figure 5G). Nevertheless, the m7G risk score was not correlated with tumor size, secretory function, and histological grade (Figure 5H-K).

High m7G risk implies the suppression of anti- tumor immune responses

The abundances of 21 immune subtypes in each ACC sample were variable (Supplement- ary Figure 2). High m7G risk was associated with decreased infiltration levels of CD8 T cells, resting CD4 memory T cells, activated NK cells, M2 macrophages, and resting mast cells. By contrast, higher infiltration levels of follicular helper T cells, M0 macrophages, and eosino- phils appeared in the m7G high-risk group com- pared with the low-risk group (Figure 6A). According to previous immunological studies, the above alterations of immune abundances are commonly detrimental to the anti-tumor immune process (Table 2). Furthermore, anti- gen presentation cell functions, check-point, cytolytic activity, and type-II IFN response were suppressed in the m7G high-risk group (Figure 6B). The immune score showed similar trends as abundances of immune cells and activities of immune pathways. Stromal score, immune score, and ESTIMATE score were markedly higher in the low-risk than in the high-risk gro- up (Figure 6C). By contrast, tumor purity was significantly higher in the high-risk than in the low-risk group (Figure 6D). Taken together, as shown in an immune heatmap (Supplementary Figure 3), different m7G risk levels were as- sociated with substantially different immune microenvironments.

The m7G risk level is associated with glycolysis and nucleotide metabolism

Metabolic reprogramming is a critical hallmark of tumor biology. Especially, glycolysis, which is a less efficient form of energy supply than oxi- dative phosphorylation, can drive tumor growth

and confer tumor cells drug resistance [35]. GSEA analyses showed that glycolysis was markedly enriched in ACC samples with high m7G risk (Figure 6E, 6F), which was conducive to ACC progression from a metabolism per- spective. Moreover, ‘Biosynthetic process’, ‘Nucleotide metabolism’, and ‘DNA replication’ were also enriched in the high-risk group (Figure 6G-I). Considering that active biosyn- thesis and nucleotide metabolism promote the occurrence and progression of cancers [36], these observations confirmed the correlations between high m7G risk scores and ACC pro- gression. Interestingly, there were no differenc- es in enrichments of FA and AA metabolisms between high- and low-risk groups (Figure 6J, 6K).

To go a step further, we analyzed the expres- sive correlations between m7G risk score and four glycolysis rate-limiting enzymes (PKM, PFKFB3, HK1 and HK2) using TCGA data. As shown in Figure 7A, HK1 expression in high m7G risk group was significantly higher than that in low m7G risk group. However, HK2 pre- sented the opposite trend, HK2 expression was lower in high m7G risk group. Besides, there were no differences in PKM and PFKFB3 expressions between two risk groups. From correlation view, HK1 and PKM expressions were positively correlated with m7G risk score (Figure 7B, 7C), whereas HK2 held negative correlation (Figure 7D). PFKFB3 expression was not correlated with m7G risk score (Figure 7E). These findings revealed that m7G risk score may herald the expressive alteration of glycolysis rate-limiting enzymes, which was the possible reason for high enrichment of glycoly- sis metabolism in high m7G risk group (Figure 6E-G).

Considering the critical roles of METTL1 in m7G modification and m7G risk score, we explored the effects of METTL1 on the expressions of above glycolysis enzymes through Western blot. The results showed that only HK1 ex- pression varied with the METTL1 expression

m7G and METTL1 in ACC

A

The immune abundances of cells

*

0.6

**

0.5

**

*

score

0.4

ns

ns

0.3

:

High

ns

Low

*

0.2


ns

ns

ns

**

ns

ns

*

ns

ns

ns

$

ns

0.1

ns

ns

0.0

B cells naive B cells memory

Plasma cells

T cells CD8

T cells CD4 naive

T cells CD4 memory resting T cells CD4 memory activated

T cells follicular helper T cells regulatory (Tregs)

T cells gamma delta

NK cells resting

NK cells activated

Monocytes Macrophages MO

Macrophages M1

Macrophages M2

Dendritic cells resting

Dendritic cells activated

Mast cells resting Mast cells activated

Eosinophils

Neutrophils

B

The activities of Immune-related pathways

C

1.0

4000

Immune Score

0.9

ns

ns

High

ns


Low

0.8

ns

2000

score

0.7

score

0.6

0

0.5

High

0.4

Low

-2000

APC co-inhibition

APC co-stimulation

Check-point

Cytolytic activity

Inflammation-promoting

Parainflammation

T cell co-inhibition

T cell co-stimulation

Type-I IFN Response

Type-II IFN Response

StromalScore ImmuneScore ESTIMATEScore

Figure 6. The effects of m7G risk score on TIM and metabolomics of ACC. A. Differential abundances of 22 immune cells between m7G high- and low-risk groups. B. Differences in activities of 10 immune signaling pathways between m7G high- and low-risk groups. C. Effects of m7G risk levels on the immune score. D. Dif- ferences in tumor purity between m7G high- and low-risk groups. E-K. Effects of m7G risk levels on the enrichments of glycolysis, biosynthetic process, nucleotide metabolism, FA metabolism, and AA metabolism. TIM, tumor immune microenvironment; FA, fatty acid; AA, amino acid; APC, antigen-presenting cell; IFN, interferon; *P< 0.05, ** P < 0.01, *** P < 0.001.

E

Enrichment plot: HALLMARK_GLYCOLYSIS

F

Enrichment plot: GO_GLYCOLYTIC_PROCESS

G

Enrichment plot: BIOSYNTHETIC_PROCESS

D

0.5

Enrichment score (ES)

0.4

P=0.017

Enrichment score (ES)

0.5

0.4

P=0.035

Enrichment score (ES)

P<0.001

Tumor Purity

0.4

0.3

0.3

0.3

1.1


0.2

0.2

0.2

1.0

0.1

0.1

0.1

0.0

0.0

0.0

0.9

-0.1

-0.1

0.8

Ranked list metric (Signal2Noise)

Ranked list metric (Signal2Noise)

Ranked list metric (Signal2Noise)

score

0.7

I (positively comelated)

If (positively comelated)

H (positively comelated)

2

2

2

4

0.6

Zero eress at 21855

4

Zero eress at 21855

1

Zero cross at 29855

0

C

D

0.5

1

‘L’ (negatively correlated)

‘L’ (negatively correlated)

L’ (negatively comelated)

0

10.000

20.000

30.000

40.000

50.000

0

10.000

20.000

30.000

40.000

50.000

0

10,000

20,000

30,000

40,000

50,000

T

T

Rank in Ordered Dataset

Rank in Ordered Dataset

Rank in Ordered Dataset

High-Risk

Low-Risk

Enrichment profile

Hits

Ranking metric scores

Enrichment profile

- Hits

Ranking metric scores

Enrichment profile - Hits

Ranking metric scores

H

Enrichment plot: WP_NUCLEOTIDE_METABOLISM

Enrichment plot: KEGG_DNA_REPLICATION

J

Enrichment plot: HALLMARK_FATTY_ACID_METABOLISM

K

Enrichment plot:

0.8

0.9

Enrichment score (ES)

AMINO_ACID_AND_DERIVATIVE_METABOLIC_PROCESS

0.7

P=0.004

Enrichment score (ES)

0.8

0.6

0.7

P<0.001

Enrichment score (ES)

0.3

P=0.439

Enrichment score (ES)

0.3

0.6

P=0.227

0.5

0.5

0.2

0.2

0.4

0.4

0.3

0.3

0.1

0.1

0.2

0.2

0.0

0.0

0.1

0.1

-0.1

0.0

0.0

-0.1

0.2

0.2

Ranked list metric (Signal2Noise)

Ranked list metric (Signal2Noise)

Ranked list metric (Signal2Noise)

Ranked list metric (Signal2Noise)

H (positively comelated)

H (positively comelated)

2

2

If (positively comelated)

2

H’ (positively cotrelated)

2

1

Zero oross at 20855

1

Zero oross at 20856

1

Zero eross at 20855

1

Zero oross at 20855

0

0

0

0

1

“L’ (negatively correlated)

“L’ (negatively comelated)

1

‘L’ (negatively correlated)

“L’ (negatively comrelated)

0

10.000

20.000

30.000

40.000

50.000

0

10.000

20.000

30.000

40.000

50.000

0

10,000

20,000

30,000

40,000

50,000

1

0

10,000

20,000

30.000

40,000

60.000

Rank in Ordered Dataset

Rank in Ordered Dataset

Rank in Ordered Dataset

Rank in Ordered Dataset

Enrichment profile - Hits

Ranking metric scores

Enrichment profile Hits

Ranking metric scores

Enrichment profile - Hits

Ranking metric scores

Enrichment profile Hits

Ranking metric scores

m7G and METTL1 in ACC

Table 2. Effects of m7G high risk on the immune microenvironment of ACC
Immune cellsChanging trendBasic functionFinal effect on anti-tumor immune
T cells CD8DecreasedCD8+ T cells can eradicate tumor cells by recognizing tumor-associated antigens.Unfavorable
T cells CD4 memory restingDecreasedMemory CD4 T cells can rapidly enhance anti- tumour activity of CTLS.Unfavorable
T cells follicular helperIncreasedTFH cells can secrete immune-protective factors but are exclusive with cytotoxic process.Uncertain
NK cells activatedDecreasedNK cells can kill tumor cells by cytotoxicity and producing IFN-y.Unfavorable
Macrophages M0IncreasedInfiltration of macrophages in solid tumors is associated with poor prognosis and correlates with chemotherapy resistance in most cancers.Unfavorable
Macrophages M2DecreasedMacrophages M2 promote tumor growth by inhibiting the functions of CD8+ T cells.Beneficial
Mast cells restingDecreasedMast cells exert the pro-oncogenic roles through releasing angiogenic factors, such as VEGF.Beneficial
EosinophilsIncreasedEosinophils can secrete pro-angiogenic and unique granule proteins, the latter factors possess anti-tumor capacities.Uncertain

m7G, N7-methylguanosine; ACC, adrenocortical carcinoma; CTLs, cytotoxic T lymphocytes; TFH, T cells follicular helper; NK, natural killer; IFN, interferon; VEGF, vascular endothelial growth factor.

Figure 7. Associations of four glycolysis rate-limiting enzymes with METTL1 and m7G risk score. A. The differences in glycolysis enzymes expressions between high- and low-risk groups (TCGA cohort). B-E. Expressive correlations between m7G risk score and four glycolysis enzymes (TCGA cohort). F. The effects of METTL1 on the expressions of four glycolysis enzymes (H295R cells).

A

TCGA cohort

B

C

Low

10

10

Relative expression

6

High

9

9

m7G risk score

m7G risk score

8

8

4

7

7

2

6

6

Spearman

5

R= 0.317

Spearman

P= 0.005

5

R = - 0.276

0

P= 0.014

PKM

HK1

HK2

PFKFB3

2

3

4

5

0

HK2 relative expression

2

4

6

HK1 relative expression

D

E

F

Control

OE-METTL1

sh-METTL1

10

10

H295R

9

9

m7G risk score

m7G risk score

HK1

102kDa

8

8

HK2

102kDa

7

7

PKM

58kDa

6

6

Spearman

Spearman

5

R= 0.238

P= 0.035

5

R = 0.206

PFKFB3

60kDa

P= 0.068

6

7

8

1

2

3

4

5

PKM relative expression

PFKFB3 relative expression

GAPDH

36kDa

change, but no expressive alterations of PKM, HK2 and PFKFB3 were observed (Figure 7F).

METTL1 overexpression could increase HK1 expression, whereas its deletion decreased

HK1 expression. Altogether, high m7G risk was associated with active glycolysis metabolism (Figure 6E-G) and METTL1, the core member in m7G risk signature, could affect the expression of glycolysis rate-limiting enzyme HK1 (Figure 7F).

High m7G risk is related to adverse genetic alterations

Somatic mutations were common in ACC sam- ples. Missense mutation was the most fre- quent mutational form (Figure 8A), and single nucleotide polymorphism (SNP) was also the dominant variant type (Figure 8B). Meanwhile, C>T (n = 3,758) and C>A (n = 3,220) substitu- tions were the major types of SNPs (Figure 8C). The mean variant of each ACC sample was as high as 21.5 (Figure 8D). Moreover, the mutations of TTN, TP53, MUC4, MUC16, and CTNNB1 frequently occurred in ACC samples (Figure 8E). Different m7G risk levels displayed different mutational characteristics (Figure 8F, 8G). The total mutation frequency in the high-risk group was up to 83.33%, whereas that in the low-risk group was only 38.74%. Moreover, the frequencies of characteristic mutated genes in the high-risk group were substantially higher than those in the low-risk group, such as TP53, CTNNB1, and MUC4. These findings indicated that high m7G risk was associated with adverse genetic mutations of ACC. Nonetheless, the somatic mutations of m7G signature genes were rarely visible in ACC samples. METTL1 exhibited the highest muta- tion frequency at 7% (Figure 8H).

m7G risk scores may serve as biomarkers of the efficacy of ICBs and mitotane treatments

No difference in m7G risk scores between patients responding and not responding to radiotherapy was observed (Figure 9A). With respect to the mitotane treatment, the most commonly used adjuvant option for ACC, m7G risk score in drug-sensitive patients was signifi- cantly higher than that in drug-resistant patients (Figure 9B).

We investigated potential linkages between m7G risk scores and ICB efficacy. TMB was markedly higher in the high-risk than in the low- risk group (Figure 9C). The TMB value was also positively correlated with m7G risk score (r = 0.561, P < 0.001; Figure 9D). High TMB is com-

monly accompanied by high production of tumor neoantigens, suggesting a good res- ponse for ICBs [17]. TIDE scores were lower in the high-risk than in the low-risk group (Figure 9E). Patients with low m7G risk were more sus- ceptible to suffer from T cell dysfunction (Figure 9E). These findings also suggested that high m7G risk may indicate good responses to ICBs. The high-risk group showed higher expres- sion of CD274 (PD-L1) than the low-risk group (Figure 9F). LAG3 expression was positively correlated with m7G risk score (r = 0.272, P = 0.015; Figure 9K). However, the expression of other ICs was not associated with m7G risk score (Figure 9G-J, 9L). Further, the IMvigor 210 cohort revealed that m7G risk scores were higher in CR/PR than in SD/PD patients (Figure 9M). The ORR in the high-risk group was 31.9%, which was also significantly higher than that in the low-risk group (15.5%; Figure 9N). Thus, a high m7G risk may indicate ICB treat- ment response.

METTL1 can promote the proliferation, migra- tion, and invasion of ACC cells

Considering that METTL1 was the most critical regulatory gene in the m7G process (Table 3), we further explored its roles in ACC progres- sion. Using 10 pairs of clinical samples, we confirmed that METTL1 was significantly upreg- ulated in ACC tissues compared to adjacent normal tissues (Figure 10A). The qPCR tests revealed that sh-METTL1 and OE-METTL1 could effectively manipulate METTL1 expres- sion (Figure 10B, 10C). Colony formation assays showed that METTL1 overexpression promoted, whereas METTL1 silencing inhibit- ed the proliferation of ACC cells (Figure 10D, 10E). The colony numbers in the OE-METTL1 group were significantly higher than that in other experimental groups; by contrast, the least colonies were observed in the sh-METTL1 group (Figure 10F, 10G). Regarding migrative and invasive abilities, upregulation METTL1 enhanced, whereas METTL1 deletion sup- pressed the migration ability of ACC cells (Figure 11A). Likewise, METTL1 stimulated invasion by ACC cells (Figure 11B). The results of quantitative analyses also were in line with these results (Figure 11C-F). Collectively, METTL1 exhibited oncogenic potential in ACC progression.

m7G and METTL1 in ACC

A

B

C

D

Variants per sample Median: 21.5

E

Variant Classification

Variant Type

SNV Class

Top 10 mutated genes

Missense_Mutation

TOG

1776

250

Nonsense_Mutation

SNP

TTN

12%

T>A

672

TP53

18%

Frame_Shift_Del

MUC4

14%

Splice_Site

1184

T>C

989

MUC16

14%

Frame_Shift_Ins

INS

CTNNB1

15%

In_Frame_Del

C>T

3758

PKHD1

9%

Nonstop_Mutation

592-

NF1

9%

C>G

1124

CNTNAPS

9%

In_Frame_Ins

DEL

PCDH15

8%

Translation_Start_Site

C>A

3220

ASXL3

8%

0

0

1000

2000

3000

4000

5000

6000

0

1000

2000

3000

4000

5000

6000

7000

0.00

0.25

0.50

0.75

1.00

0

7

=

&

F

Altered in 25 (83.33%) of 30 samples.

G

Altered in 19 (38.78%) of 49 samples.

943

290

TMB

TMB

0

No. of samples

9

0

6

0

0

No. of samples

TP53

30%

TP53

8%

CTNNB1

27%

CTNNB1

10%

MUC16

17%

MUC16

12%

MUC4

30%

MUC4

6%

TTN

20%

TTN

6%

PKHD1

17%

PKHD1

4%

CNTNAP5

17%

CNTNAP5

4%

ASXL3

13%

ASXL3

2%

DST

13%

DST

2%

HMCN1

10%

HMCN1

8%

Risk

Risk

Missense_Mutation

Frame_Shift_Del

Risk

Nonsense_Mutation

Risk

Frame_Shift_Ins

Nonsense_Mutation

Missense_Mutation

In_Frame_Del

· Frame_Shift_Del

· Multi_Hit

high

· Multi_Hit

high

low

low

H

Study of origin

Mutation spectrum

METTL1

7%*

NCBP1

3%*

NUDT1

3%*

NUDT5

3%*

Genetic Alteration

Inframe Mutation (unknown significance)

Amplification

Deep Deletion

No alterations

Not profiled

Study of origin

Adrenocortical Carcinoma (TCGA, Firehose Legacy)

Adrenocortical Carcinoma (TCGA, PanCancer Atlas)

Mutation spectrum

C>A

C>G

C>T

T>A

T>C

T>G

No data

m7G and METTL1 in ACC

Figure 8. Somatic mutation information of ACC and m7G signature genes. A. Variant classification of ACC samples in the TCGA cohort. B. Variant type of ACC samples in the TCGA cohort. C. SNV class of ACC samples in the TCGA cohort. D. Somatic variants of each ACC sample in the TCGA cohort. E. The top 10 genes with the highest mutational frequency. F, G. Mutational waterfall plots of high- and low-m7G risk groups. The top column shows the TMB of each ACC sample. The characteristic mutation genes are listed on the left side of the plots and their mutation frequencies are listed on the right side. H. The somatic mutation frequency of m7G signature genes based on cBioPortal database. TMB, tumor mutation burden; SNV, single nucleotide variation; SNP, single nucleotide polymorphism; INS, insertion; DEL, deletion.

A

Radiation Therapy

B

Mitotane therapy

C

Tumor mutation burden

D

ns

10

10

25

25

Spearman

r = 0.561

9

9

20

20

P < 0.001

m7G risk score

m7G risk score

8

8

TMB value

15

15

TMB

7

7

10

10

6

5

6

5

0

5

5

0

-5

Non-response

Response

Non-response

Response

High-Risk

Low-Risk

6

8

10

m7G risk score

E

TIDE analysis

F

ns

15

ns

1.0

Expression

10

score

0.5

High

High

Low

Low

5

0.0

ns

ns

I

ns

0

I

I

I

±

I

-0.5

TIDE

Dysfunction

Exclusion

CD274 CTLA4 BTLA HAVCR2 TIGIT LAG3

Figure 9. Therapeutic correlations of m7G risk scores. A. Differences in m7G risk scores between ACC patients with response and non-response to radiotherapy based on the TCGA cohort. B. Difference in m7G risk scores between mitotane-response and -nonresponse patients based on the TCGA cohort. C. Differences in TMB values between high- and low-m7G risk groups. D. Correlation between m7G risk scores and TMB. E. TIDE scoring results. F. Expression differences of six ICs between high- and low-m7G risk groups. G-L. Expression correlations between ICs and m7G risk score. M. Differences in m7G risk score between ICB-response and -nonresponse patients based on the IMvigor 210 cohort. N. The ORR of each m7G risk group based on the IMvigor 210 cohort. TMB, tumor mutation burden; ICS, immune checkpoints; ORR, objective response rate; CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease; * P < 0.05, ** P < 0.01, *** P < 0.001.

G

H

I

J

12

10

10

10

10

9

9

9

m7G risk score

m7G risk score

m7G risk score

m7G risk score

8

8

8

8

7

7

7

6

Spearman

6

r = - 0.157

Spearman

6

Spearman

6

Spearman

P = 0.166

5

r =- 0.155

P = 0.174

5

r = 0.110

P = 0.334

5

r = 0.041

P = 0.717

4

0.0

2.5

5.0

7.5

10.0

0.0

0.5

1.0

1.5

2.0

0.0

0.1

0.2

0.3

0

5

10

CD274 expression

CTLA4 expression

BTLA expression

HAVCR2 expression

M

N

K

L

IMvigor 210 cohort

IMvigor 210 cohort

1.0

11

6.0

31.9%

15.5%

12

10

0.8

m7G risk score

5.5

9

Proportion

m7G risk score

10

m7G risk score

0.6

CR/PR

8

5.0

SD/PD

8

0.4

7

P=0.027

4.5

6

Spearman

0.2

r = 0.272

6

Spearman

P = 0.015

5

r =- 0.176

4.0

0.0

P = 0.121

4

CR/PR

SD/PD

High-risk

Low-risk

0

5

10

15

0.0

0.5

1.0

1.5

2.0

LAG3 expression

TIGIT expression

Table 3. Biofunctions of m7G signature genes in various cancers
m7G Signature geneStudyCancer typeMain function
METTL1PMID: 34371184LC, ICC, HCCPromote cancer progression
PMID: 34352206
PMID: 34898034
NCBP1PMID: 31448526LUADPromote cancer progression
NUDT1PMID: 29075149GC, LUADPromote cancer progression
PMID: 21289483
NUDT5PMID: 33096144 PMID: 35247377NSCLC, GCPromote cancer progression

m7G, N7-methylguanosine; LC, lung cancer; ICC, intrahepatic cholangiocarcinoma; HCC, hepatocellular carcinoma; LUAD, lung adenocarcinoma; GC, gastric cancer; NSCLC, non-small-cell lung cancer.

Figure 10. Overexpression of METTL1 promotes the proliferation of ACC cancer cells. A. 10 pairs of clinical samples confirmed the expressive differences of METTL1 between adjacent normal and ACC tissues through qPCR. B, C. Tests of transfection efficiency of sh-METTL1 and OE-METTL1. D-G. Colony formation assays of each experimental group. sh-METTL1, the short hairpin RNA target METTL1; OE-METTL1, overexpression of METTL1; * P < 0.05, ** P < 0.01, *** P < 0.001.

A

B

C

6


The relative expression of METTL1

5

H295R

4

SW13

4

Relative METTL1

4


3

expression

Relative METTL1 expression

3

2

2

2

1

1

**

0

0

0

Adjacent normal Tumor

Blank

OE-Vector

OE-METTL1

sh-Vector

sh-METTL1

Blank

OE-Vector

OE-METTL1

sh-Vector

sh-METTL1

D

Blank

OE-vector

OE-METTL1

sh-vector

sh-METTL1

H295R

E

SW13

F

600

H295R

G

800

SW13

Colony numbers

Colony numbers

400

600


400

200


200


0

0

Blank

OE-vector

OE-METTLI

sh-vector

sh-METTL1

Blank

OE-vector

OE-METTLI

sh-vector

sh-METTL1

Figure 11. Overexpression of METTL1 promotes the migration and invasion of ACC cancer cells. A. The transwell- migration results of each experimental group. B. Transwell-invasion results of each experimental group. C, D. Quanti- tative statistics of migrative cells. E, F. Quantitative statistics of invasive cells. * P < 0.05, ** P < 0.01, *** P < 0.001.

A

Blank

OE-vector

OE-METTL1

sh-vector

sh-METTL1

20%

20x

20元

20x

H295R

Migration

50pm

Dum

50um

50um

50pm

20

20x

20x

SW13

50pm

50pm

B

Blank

OE-vector

OE-METTL1

sh-vector

sh-METTL1

20

20%

20

H295R

Invasion

50jim

4

50pm 2

4

Dum

50

20-

20x

20x

20x

20x

SW13

50pm

50pm

50pm

50um

50um

C

400

H295R

D

300

SW13

E

250

H295R

=

F

250


SW13

=

200

200-

Migrative cell numbers

300-

Migrative cell

numbers

200

Invasive cell

Invasive cell

T

numbers

150

numbers

150

200-

T

100-

100

100-

100-

=


50

50

0

0

0

0

Blank

OE-vector

OE-METTLI

sh-vector

sh-METTLI

Blank

OE-vector

OE-METTLI

sh-vector

sh-METTLI

Blank

OE-vector

OE-METTLI

sh-vector

sh-METTL1

Blank

OE-vector

OE-METTLI

sh-vector

sh-METTL1

METTL1 expression affects the infiltration lev- els of CD8+ T cells and macrophages in ACC tissues

ssGSEA results revealed that as the core mem- ber of the m7G risk signature, METTL1 expres- sion was negatively correlated with the infiltrat- ing levels of CD8+ T cells (Figure 12A), and it was positively correlated with that of macro- phages (Figure 12B). Further, the somatic copy number alteration (SCNA) of METTL1 was also associated with the infiltration levels of CD8+ T cells and macrophages (Figure 12C). Arm-level deletion of METTL1 was accompa- nied by the increased abundance of CD8+ T cells and decreased abundance of macro- phages (Figure 12C).

We then confirmed the effects of METTL1 on infiltration levels of immune cells through

immunofluorescence assays. ACC sample with high-expression METTL1 showed very low fluo- rescence intensity of CD8+ T cells (red), but conspicuous that of METTL1 (green). In con- trast, ACC samples with low-expression METTL1 presented high fluorescence intensity of CD8+ T cells (red; Figure 12D). Regarding macrophages, the clinical samples exhibited the opposite trend to the fluorescence inten- sity of CD8+ T cells. The fluorescence intensi- ties of macrophages (CD163, red) and METTL1 (green) were both strong in ACC samples with high METTL1 expression, whereas in ACC samples with low METTL1 expression, they were weak (Figure 12E). Hence, METTL1 did not only stimulate the malignant behaviors of ACC cells but also affected the infiltration lev- els of CD8+ T cells and macrophages in ACC tissues.

m7G and METTL1 in ACC

A

B

C

ACC

0.62

0.62

0.5

**

*


Enrichment of CD8+ T cells

Enrichment of Macrophages

0.60

0.60

Infiltration Level

0.4

Copy Number

0.58

0.58

Arm-level Deletion

0.3

Diploid/Normal

Arm-level Gain

0.56

0.56

0.2

High Amplication

0.54

Spearman

0.54

Spearman

r =- 0.264

r = 0.286

0.1

P = 0.019

P = 0.011

0.52

0.52

4.0

4.5

5.0

5.5

6.0

6.5

2

4

6

B Cell

CD8+ T Cell

CD4+ T Cell

Macrophage

Neutrophil

Dendritic Cell

The relative expression of METTL1

The relative expression of METTL1

D

DAPI

CD8

METTL1

Merge

40x

40x

40x

40x

High

50pm

50um

50um

50um

40x

40x

40x

40x

Low

50um

50um

50um

50um

Figure 12. The immune effects of METTL1. A, B. Associations of METTL1 expression with immune abundances of CD8+ T cells and macrophages based on the TCGC cohort. C. The relationships between METTL1 SCNA and the infiltration levels of six core immune cells based on the TIMER database. D. Immunofluorescence tests on two ACC clinical samples showed the infiltration levels of CD8+ T cells under different METTL1 expression levels. E. Immunofluorescence tests on two ACC clini- cal samples showed the infiltration levels of macrophages (CD163) under different METTL1 expression levels. SCNA, somatic copy number alteration; * P < 0.05, ** P < 0.01, *** P < 0.001.

E

DAPI

CD163

METTL1

Merge

40x

40x

40x

40x

High

50um

50um

50μm

50um

40x

40x

40x

40x

Low

50um

50um

50μm

50um

Figure 13. Effects of METTL1 expression on xenograft tumors. A-C. Silencing METTL1 suppresses ACC tumor growth in a xenograft model. D. Difference of tumor weight between sh-vector and sh-METTL1 groups. E. Difference of tu- mor volume between sh-vector and sh-METTL1 groups. F, G. qPCR assays showed the differences of METTL1/P53 expression between sh-vector and sh-METTL1 groups. * P < 0.05, ** P < 0.01, *** P < 0.001.

A

D

800

P=1.5e-07

sh-vector

Tumor weight (mg)

600

400

200

0

sh-vector

sh-METTL1

E

800


1

2

3

4

5

6

Tumor volume(mm3)

600

sh-vector

sh-METTLI

B

400

**

sh-METTL1

200

**

0

0

5

10

15

20

Days

F

Relative METTL1 expression

1.5


1.0

1

2

3

4

5

6

C

0.5

sh-vector

0.0

sh-vector sh-METTL1

G

4


sh-METTL1

Relative P53 expression

3

1

2

3

4

5

6

2

5

7

8

9

10

11

12

13

14

15

1

0

T

sh-vector

sh-METTL1

Silencing METTL1 suppresses tumor growth in a xenograft model

Visually, the tumor burden of nude mice in the sh-METTL1 group was lower than that in the negative control group (Figure 13A, 13B). After the mice were sacrificed, we confirmed that silencing METTL1 indeed suppressed xeno-

graft tumor growth (Figure 13C). Tumor weight in the sh-METTL1 group was significantly lower than that in the sh-vector group (Figure 13D), and tumor volume exhibited the same trend (Figure 13E). qPCR revealed that METTL1 expressions of xenograft tumors in the sh-MET- TL1 group were significantly lower than that in the sh-vector group (Figure 13F). However, the

mRNA levels of P53, a classical tumor suppres- sor gene, were substantially higher in the sh- METTL1 than in the sh-vector group (Figure 13G). These results highlighted that METTL1 deletion decelerated ACC growth and increased P53 expression.

Potential regulatory mechanisms of METTL1 in ACC progression

Using TargetScanHuman, miRDB and ENCORI databases, we predicted potential upstream miRNAs of METTL1. The intersection part of three databases was obtained through a Venn diagram, miR-1277-3p and miR-885-5p were screened as the candidate regulators (Figure 14A). Next, the minimum free energy (MFE) of these miRNAs was quantified via RNAhybrid database. MEF of miR-1277-3p and miR-885- 5p was -15.5 and -27.4 kcal/mol, indicating the latter was more accessible to bind to METTL1 (Figure 14B, 14C). The binding site between miR-885-5p and METTL1 was also predicted with the aid of TargetScanHuman database. As shown in Figure 14D, miR-885- 5p may target the 3’-UTR region of METTL1 namely 5’-GUAAUGGA-3’.

Furthermore, the potential transcription factor (TF) of METTL1 was also investigated. CEBPB was speculated as the upstream TF of METTL1 based the intersection of hTFtarget, Gene- cards and ALGGEN databases (Figure 14E). The motif sequence of CEBPB exhibited the specificity and conservation of its binding site (Figure 14F). Theoretically, the promoter sequence with the highest binding probability with CEBPB was 5’-TATTGCACAAT-3’ (Figure 14F). Using JASPAR database, we predicted the binding site between CEBPB and METTL1 (Figure 14G). the most probable site was locat- ed between the 479th and 489th bases upstream of the METT1 transcription starting site (TSS), and the sequence was 5’-CGTTT- CACCAT-3’ (Figure 14H). Collectively, miR-885- 5p and CEBPB may participate ACC progres- sion through regulating METTL1.

Discussion

ACC is a rare urological carcinoma with an inci- dence of 0.7-2.0/million [37]. Due to the high degree of malignancy and early metastases, the five-year OSR of ACC patients is commonly less than 20%. To maximize the patients’ sur-

vival outcomes, multiple therapies such as molecular target treatment (MTT) and ICBs were explored for use in ACC treatment. Nonetheless, MTT produces only negligible results [38], and selecting suitable cases for ICBs is a persistent problem. By contrast to m6A, m7G has not received sufficient attention although it may exert important functions dur- ing cancer regulation and treatment [39-41].

Prognostic evaluation is vital for deciding on therapeutic strategies, however, traditional clinicopathological indicators of ACC do not allow for accurate prognosis. TNM-staging, clin- ical stage, Ki-67, and histological grade can be used for stratifying patient survival outcomes, however, up to 25% of patients experience a different outcome than predicted [42]. Hence, better indicators are required to compensate for the deficiency of current prognostic meth- ods. In the present study, we established a novel m7G risk signature, and the m7G risk score showed outstanding predictive accuracy for OSR and was identified as the only indepen- dent ACC-prognostic factor. More importantly, m7G risk score could improve the predictive accuracy and making-decision benefit of tradi- tional AJCC-Stage prognostic system (Figure 4C, 4E). These findings assured us that m7G risk score possessed highly prognostic value.

RNA methylation profoundly affects the anti- cancer immune response and the efficacy of immunotherapy [9, 43]. For instance, methyl- transferases METTL3/14 can enhance the response to anti-PD-1 treatment in colorectal cancer (CRC) and melanoma [44]. The activa- tion of retinoic acid-inducible gene-I (RIG-I) which is an innate immune receptor and is responsible for triggering type-I IFN response, relies on m7G recognition [45]. Moreover, METTL1/WDR4-mediated tRNA m7G can af- fect the immune landscape of head and neck squamous cell carcinoma (HNSCC) by altering the proportion of CD8+ T cells, NK cells, and CD4+ T cells [33]. In the current study, we also confirmed that the m7G risk score was strongly correlated to the immune microenvironment of ACC. High m7G risk significantly suppressed the immune enrichment of CD8+ T cells and NK cells, but it stimulated that of macrophages and TFH. Regarding the most potent anti-tumor immune cells, the functions of CD8+ T cells and NK cells in eradicating tumor cells did not need

A

B

TargetScanHuman

RNAhybrid

Target: METTL1(NM005371.6)

127

Length: 1378

5’

MiRNA: miR-1277-3p

37

6

2

Length: 22

MFE: - 15.5kcal/mol

6

1

23

C

5’

miRDB

ENCORI

Target: METTL1(NM005371.6)

Length: 1378

MiRNA: miR-885-5p

miR-1277-3p miR-885-5p

Length: 22

MFE: - 27.4kcal/mol

D

E

hTFtarget

F

2.0

ATT&CASAAz

44

1.5

20

0

1

0.5

239

2

20

0.0

1

2

3

4

5

6

7

8

9

10

11

H

GeneCards

ALGGEN

5’-CGTTTCACCAT-3’

Predictive Binding Site

CEBPB

-2000

-489 ~ - 479

+1

a JASPAR

METTL1 TSS

0 TargetScanHuman Prediction of microRNA forgetsPredicted consequential pairing of target region (top) and miRNA (bottom)Site typeContext++ scoreContext++ score percentileWeighted context++ scoreConserved branch lengthPCTPredicted relative Kp
Position 84-91 of METTL1 3' UTR hsa-miR-885-5p5' GGAAGAAAGUUCUACGUAAUGGA. 3' UCUCCGUCCCAUCACAUUACCU8mer-0.4099-0.400.072N/A-3.286

G

NameScoreRelative score YSequence IDStartEndStrand APredicted sequence 4
MA0466.1.CEBPB10.4993430.9299949439183495hg38_knownGene_ENST00000324871.12479489-CGTTTCACCAT
MA0466.2.CEBPB5.91766930.8649728020946857hg38_knownGene_ENST00000324871.12479488+ATGGTGAAAC
MA0466.2.CEBPB5.6514920.861368825631875hg38_knownGene_ENST00000324871.124453+TTTGTGAAAC
MA0466.2.CEBPB5.44399550.8585593698214778hg38_knownGene_ENST00000324871.12525534-CTTGCGCTAC
MA0466.2.CEBPB5.00903840.8526701538325228hg38_knownGene_ENST00000324871.12479488-GTTTCACCAT

Figure 14. Potential regulatory mechanism of METTL1 in ACC. A. The predicting intersection of three miRNA da- tabases. B, C. The minimum free energy of miR-1277-3p and miR-885-5p based on RNAhybrid database. D. The predicted binding site between miR-885-5p and 3’-UTR region of METTL1 based on TargetScanHuman database. E. The predicting intersection of three TF databases. F. The motif sequence of CEBPB. G. Five predicted binding sites between CEBPB and promoter region of METTL1 based on JASPAR database. H. CEBPB binding site with highest predictive score. MFE, minimum free energy; TSS, transcription starting site.

to further elaborate [46]. Macrophages medi- ate immunotolerance and immune evasion through releasing CCL2, CCL5, and VEGF cytokines [47, 48]. However, the roles of TFH cells in immune regulation is more complex. As TFH cells can produce CXCL13 which exerts immune-protective functions, they are strongly associated with long survival time of patients with breast cancer [49]. Nevertheless, TFH cells and cytotoxic transcriptional programs are functionally exclusive, thus TFH cells may be detrimental to anti-tumor immune and ICB ther- apy [50]. In view of these facts, the effects of m7G risk on the immune microenvironment of ACC are complicated and multifaceted.

ICBs represent a revolutionary change in can- cer treatment, however, identifying suitable cases is challenging. Currently, several bio- markers and methods have been tested to predict the efficacy of ICIs, such as TMB [51], IC expression [52], and TIDE scoring [19]. Surprisingly, the m7G risk score was associated with all the above predictive markers, which demonstrated its potential for predicting ICIs therapeutic response. Although PD-L1 expres- sion and TMB may each inform on the use of ICIs in most cancers [52], considerable con- troversy on these biomarkers remains. For instance, low TMB does not preclude respons- es to ICIs, especially in patients with Kaposi sarcoma [53] and Merkel cell carcinoma [54]. Moreover, experimental determination of TMB requires whole exome sequencing, which is technically complicated and highly expensive, thus limiting its clinical applicability [51]. Therefore, the m7G risk score sheds new light on ICIs prediction.

The catalytic function of METTL1 is a prerequi- site of the m7G process [8]. As expected, METTL1 was a member of the m7G risk signa- ture, which was consistent with the core identi- ty of METTL1 in m7G. Recent studies confirmed its pivotal roles in cancer onset and develop- ment. For example, METTL1/WDR4-mediated m7G can promote the progression of lung can- cer [36], and METTL1 drives oncogenic trans- formation through accelerating the m7G mo-

dification of Arg-TCT tRNA [55]. The METTL1- m7G-EGFR axis facilitates the progression of bladder cancer [56]. In the present study, we confirmed the oncogenic potential of METTL1 in ACC for the first time. METTL1 overexpres- sion substantially enhanced the proliferation, migration, and invasion abilities of ACC cells. Moreover, METTL1 deletion significantly sup- pressed xenograft tumor growth. These find- ings thus confirmed its potential as a tumor therapeutic target.

Glycolysis termed ‘Warburg effect’ can satisfy the metabolic requirements of cell proliferation and regulate cancer metastasis, thus acting as a pivotal hallmark of solid cancers [57]. In the present study, we found that glycolysis enrichment was concomitant with high-m7G risk scores, indicating m7G modification may drive glycolysis metabolism. However, only HK1 was affected by METTL1 among four glycolysis rate-limiting enzymes. Tissue-specific expres- sions of these glycolytic enzymes were the pos- sible reasons. Duan K et al. have confirmed that although HK1 and PKM were both upregu- lated in ACC compared to normal adrenal corti- cal tissue and adrenal cortical adenoma (ACA), PKM expression was overall low in ACC [58]. Moreover, HK2 mainly expressed in insulin-sen- sitive tissues, such as colon and fat, but not adrenal [59]. PFKFB3 mainly expressed the cancers of brain, skeletal muscle and liver [60]. In light of these facts, PKM, HK2 and PFKFB3 may rarely express in ACC compared to HK1, leading their expressions not to be affected by METTL1.

There are some limitations to this study. First, the m7G risk signature remains to be tested in a real clinical cohort. Second, since the detec- tion of m7G modifications have a certain degree of difficulty, it is intractable to determine the specific relationships between m7G risk score and m7G modification level. Third, due to the fact that genetic mutation analysis was reliant on the whole-exome next-generation sequenc- ing, we were unable to validate mutational fea- tures caused by m7G risk score in our current experimental condition. Fourth, the specific

m7G and METTL1 in ACC

mechanisms of METTL1 in ACC progression remain experimental validation. Hence, further research is necessary.

Conclusions

Although m7G is one of the most frequent RNA modifications, its roles in ACC remain obscure. Herein, we developed a novel m7G risk signa- ture for ACC clinical assessments. m7G risk score acted a biomarker for assessing progno- sis, anti-tumor immune response, glycolysis metabolism and predicting the efficacy of ICBs and mitotane treatments. It greatly con- tributes to acquire the disease state of ACC patients, in turn advancing individualized thera- py. Moreover, we confirmed the pro-oncogenic roles of METTL1 in ACC progression, which highlighted its great potentials as a novel anti- cancer target. In conclusion, m7G, an unre- solved epigenetic aspect, is expected to advance the paradigm of ACC treatment and clinical assessment.

Acknowledgements

All authors would like to thank Second Affiliated Hospital of Xi’an Jiaotong University for its sup- port. Dr. Fangshi Xu would like to thank his wife Dr. Danrui Cai for her encourage. All authors also thank Bullet Edits Limited for the linguistic editing of the manuscript. The animal study was reviewed and approved by the Ethics Committee of the Second Affiliated Hospital of Xi’an Jiaotong University.

For METTL1 testing on tumor samples, all patients provided written informed consent.

Disclosure of conflict of interest

None.

Address correspondence to: Xueyi Li and Bin- cheng Ren, Department of Rheumatology and Immunology, Second Affiliated Hospital of Xi’an Jiaotong University, No. 157, West Five Road, Xi’an 710004, Shaanxi, China. Tel: +86-29-876793- 23; E-mail: 13992891987@139.com (XYL); ren- bincheng7@163.com (BCR)

References

[1] Ettaieb M, Kerkhofs T, van Engeland M and Haak H. Past, present and future of epigenetics

in adrenocortical carcinoma. Cancers (Basel) 2020; 12: 1218.

[2] Pittaway JFH and Guasti L. Pathobiology and genetics of adrenocortical carcinoma. J Mol Endocrinol 2019; 62: R105-R119.

[3] Shariq OA and Mckenzie TJ. Adrenocortical carcinoma: current state of the art, ongoing controversies, and future directions in diagno- sis and treatment. Ther Adv Chronic Dis 2021; 12: 20406223211033103.

[4] Fassnacht M, Terzolo M, Allolio B, Baudin E, Haak H, Berruti A, Welin S, Schade-Brittinger C, Lacroix A, Jarzab B, Sorbye H, Torpy DJ, Stepan V, Schteingart DE, Arlt W, Kroiss M, Leboulleux S, Sperone P, Sundin A, Hermsen I, Hahner S, Willenberg HS, Tabarin A, Quinkler M, de la Fouchardière C, Schlumberger M, Mantero F, Weismann D, Beuschlein F, Gelderblom H, Wilmink H, Sender M, Edgerly M, Kenn W, Fojo T, Müller HH and Skogseid B. Combination che- motherapy in advanced adrenocortical carci- noma. N Engl J Med 2012; 366: 2189-2197.

[5] Xu F, Guan Y, Ma Y, Xue L, Zhang P, Yang X and Chong T. Bioinformatic analyses and experi- mental validation of the role of m6A RNA meth- ylation regulators in progression and prognosis of adrenocortical carcinoma. Aging (Albany NY) 2021; 13: 11919-11941.

[6] Shen C, Liu J, Yang X, Jiao W and Wang Y. De- velopment and validation of an m6A RNA methylation regulators-based signature for predicting the prognosis of adrenocortical car- cinoma. Front Endocrinol (Lausanne) 2021; 12: 568397.

[7] Jin Y, Wang Z, He D, Zhu Y, Hu X, Gong L, Xiao M, Chen X, Cheng Y and Cao K. Analysis of m6A-related signatures in the tumor immune microenvironment and identification of clinical prognostic regulators in adrenocortical carci- noma. Front Immunol 2021; 12: 637933.

[8] Tomikawa C. 7-methylguanosine modifications in transfer RNA (tRNA). Int J Mol Sci 2018; 19: 4080.

[9] Zhang M, Song J, Yuan W, Zhang W and Sun Z. Roles of RNA methylation on tumor immunity and clinical implications. Front Immunol 2021; 12: 641507.

[10] Ma J, Han H, Huang Y, Yang C, Zheng S, Cai T, Bi J, Huang X, Liu R, Huang L, Luo Y, Li W and Lin S. METTL1/WDR4-mediated m(7)G tRNA modifications and m(7)G codon usage pro- mote mRNA translation and lung cancer pro- gression. Mol Ther 2021; 29: 3422-3435.

[11] Zhou Y, Zhou B, Pache L, Chang M, Khoda- bakhshi AH, Tanaseichuk O, Benner C and Chanda SK. Metascape provides a biologist- oriented resource for the analysis of systems- level datasets. Nat Commun 2019; 10: 1523.

[12] Budczies J, Klauschen F, Sinn BV, Győrffy B, Schmitt WD, Darb-Esfahani S and Denkert C. Cutoff Finder: a comprehensive and straight- forward Web application enabling rapid bio- marker cutoff optimization. PLoS One 2012; 7: e51862.

[13] Chen B, Khodadoust MS, Liu CL, Newman AM and Alizadeh AA. Profiling tumor infiltrating im- mune cells with CIBERSORT. Methods Mol Biol 2018; 1711: 243-259.

[14] Jin Y, Wang Z, He D, Zhu Y, Chen X and Cao K. Identification of novel subtypes based on ssG- SEA in immune-related prognostic signature for tongue squamous cell carcinoma. Cancer Med 2021; 10: 8693-8707.

[15] Luo J, Xie Y, Zheng Y, Wang C, Qi F, Hu J and Xu Y. Comprehensive insights on pivotal prognos- tic signature involved in clear cell renal cell carcinoma microenvironment using the ESTI- MATE algorithm. Cancer Med 2020; 9: 4310- 4323.

[16] Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B and Liu XS. TIMER: a web server for compre- hensive analysis of tumor-infiltrating immune cells. Cancer Res 2017; 77: e108-e110.

[17] Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A and Peters S. Devel- opment of tumor mutation burden as an im- munotherapy biomarker: utility for the oncolo- gy clinic. Ann Oncol 2019; 30: 44-56.

[18] Shi Y, Lei Y, Liu L, Zhang S, Wang W, Zhao J, Zhao S, Dong X, Yao M, Wang K and Zhou Q. Integration of comprehensive genomic profil- ing, tumor mutational burden, and PD-L1 ex- pression to identify novel biomarkers of immu- notherapy in non-small cell lung cancer. Cancer Med 2021; 10: 2216-2231.

[19] Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B, Liu J, Freeman GJ, Brown MA, Wucherpfennig KW and Liu XS. Signatures of T cell dysfunction and exclusion predict can- cer immunotherapy response. Nat Med 2018; 24: 1550-1558.

[20] Brueckl WM, Ficker JH and Zeitler G. Clinically relevant prognostic and predictive markers for immune-checkpoint-inhibitor (ICI) therapy in non-small cell lung cancer (NSCLC). BMC Can- cer 2020; 20: 1185.

[21] Balar AV, Galsky MD, Rosenberg JE, Powles T, Petrylak DP, Bellmunt J, Loriot Y, Necchi A, Hoffman-Censits J, Perez-Gracia JL, Dawson NA, van der Heijden MS, Dreicer R, Srinivas S, Retz MM, Joseph RW, Drakaki A, Vaishampay- an UN, Sridhar SS, Quinn DI, Durán I, Shaffer DR, Eigl BJ, Grivas PD, Yu EY, Li S, Kadel EE 3rd, Boyd Z, Bourgon R, Hegde PS, Maria- thasan S, Thåström A, Abidoye OO, Fine GD and Bajorin DF. Atezolizumab as first-line treat- ment in cisplatin-ineligible patients with locally

advanced and metastatic urothelial carcino- ma: a single-arm, multicentre, phase 2 trial. Lancet 2017; 389: 67-76.

[22] Chen Y and Wang X. miRDB: an online data- base for prediction of functional microRNA tar- gets. Nucleic Acids Res 2020; 48: D127-D131.

[23] McGeary SE, Lin KS, Shi CY, Pham TM, Bisaria N, Kelley GM and Bartel DP. The biochemical basis of microRNA targeting efficacy. Science 2019; 366: eaav1741.

[24] Li JH, Liu S, Zhou H, Qu LH and Yang JH. star- Base v2.0: decoding miRNA-ceRNA, miRNA- ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res 2014; 42: D92-97.

[25] Rehmsmeier M, Steffen P, Hochsmann M and Giegerich R. Fast and effective prediction of microRNA/target duplexes. RNA 2004; 10: 1507-1517.

[26] Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, Stein TI, Nudel R, Lieder I, Mazor Y, Kaplan S, Dahary D, Warshawsky D, Guan-Golan Y, Kohn A, Rappaport N, Safran M and Lancet D. The GeneCards Suite: from gene data mining to disease genome sequence analyses. Curr Protoc Bioinformatics 2016; 54: 1.30.31-31.30.33.

[27] Messeguer X, Escudero R, Farré D, Núñez O, Martínez J and Alba MM. PROMO: detection of known transcription regulatory elements us- ing species-tailored searches. Bioinformatics 2002; 18: 333-334.

[28] Zhang Q, Liu W, Zhang HM, Xie GY, Miao YR, Xia M and Guo AY. hTFtarget: a comprehensive da- tabase for regulations of human transcription factors and their targets. Genomics Proteomics Bioinformatics 2020; 18: 120-128.

[29] Karolchik D, Hinrichs AS, Furey TS, Roskin KM, Sugnet CW, Haussler D and Kent WJ. The UCSC Table Browser data retrieval tool. Nucleic Acids Res 2004; 32: D493-496.

[30] Castro-Mondragon JA, Riudavets-Puig R, Rau- luseviciute I, Lemma RB, Turchi L, Blanc-Ma- thieu R, Lucas J, Boddie P, Khan A, Manosalva Pérez N, Fornes O, Leung TY, Aguirre A, Ham- mal F, Schmelter D, Baranasic D, Ballester B, Sandelin A, Lenhard B, Vandepoele K, Wasser- man WW, Parcy F and Mathelier A. JASPAR 2022: the 9th release of the open-access da- tabase of transcription factor binding profiles. Nucleic Acids Res 2022; 50: D165-D173.

[31] Xu F, Wang H, Pei H, Zhang Z, Liu L, Tang L, Wang S and Ren BC. SLC1A5 prefers to play as an accomplice rather than an opponent in pan- creatic adenocarcinoma. Front Cell Dev Biol 2022; 10: 800925.

[32] Jing J, Sun J, Wu Y, Zhang N, Liu C, Chen S, Li W, Hong C, Xu B and Chen M. AQP9 is a prog- nostic factor for kidney cancer and a promising

indicator for M2 TAM polarization and CD8+ T- cell recruitment. Front Oncol 2021; 11: 770565.

[33]

Chen J, Li K, Chen J, Wang X, Ling R, Cheng M, Chen Z, Chen F, He Q, Li S, Zhang C, Jiang Y, Chen Q, Wang A and Chen D. Aberrant transla- tion regulated by METTL1/WDR4-mediated tRNA N7-methylguanosine modification drives head and neck squamous cell carcinoma pro- gression. Cancer Commun (Lond) 2022; 42: 223-244.

[34] Xia P, Zhang H, Xu K, Jiang X, Gao M, Wang G, Liu Y, Yao Y, Chen X, Ma W, Zhang Z and Yuan Y. MYC-targeted WDR4 promotes proliferation, metastasis, and sorafenib resistance by induc- ing CCNB1 translation in hepatocellular carci- noma. Cell Death Dis 2021; 12: 691.

[35] Ganapathy-Kanniappan S and Geschwind JF. Tumor glycolysis as a target for cancer therapy: progress and prospects. Mol Cancer 2013; 12: 152.

[36] Ma J, Zhong M, Xiong Y, Gao Z, Wu Z, Liu Y and Hong X. Emerging roles of nucleotide metabo- lism in cancer development: progress and prospect. Aging (Albany NY) 2021; 13: 13349- 13358.

[37] Else T, Kim AC, Sabolch A, Raymond VM, Kan- dathil A, Caoili EM, Jolly S, Miller BS, Giordano TJ and Hammer GD. Adrenocortical carcinoma. Endocr Rev 2014; 35: 282-326.

[38] Jasim S and Habra MA. Management of adre- nocortical carcinoma. Curr Oncol Rep 2019; 21: 20.

[39] Han H, Yang C, Ma J, Zhang S, Zheng S, Ling R, Sun K, Guo S, Huang B, Liang Y, Wang L, Chen S, Wang Z, Wei W, Huang Y, Peng H, Jiang YZ, Choe J and Lin S. N(7)-methylguanosine tRNA modification promotes esophageal squamous cell carcinoma tumorigenesis via the RPTOR/ ULK1/autophagy axis. Nat Commun 2022; 13: 1478.

[40] Chen B, Jiang W, Huang Y, Zhang J, Yu P, Wu L and Peng H. N(7)-methylguanosine tRNA modi- fication promotes tumorigenesis and chemore- sistance through WNT/ß-catenin pathway in nasopharyngeal carcinoma. Oncogene 2022; 41: 2239-2253.

[41] Ouyang W, Jiang Y, Bu S, Tang T, Huang L, Chen M, Tan Y, Ou Q, Mao L, Mai Y, Yao H, Yu Y and Lin X. A prognostic risk score based on hypox- ia-, immunity-, and epithelialto-mesenchymal transition-related genes for the prognosis and immunotherapy response of lung adenocarci- noma. Front Cell Dev Biol 2021; 9: 758777.

[42] Clay MR, Pinto EM, Fishbein L, Else T and Kisel- jak-Vassiliades K. Pathological and genetic stratification for management of adrenocorti- cal carcinoma. J Clin Endocrinol Metab 2022; 107: 1159-1169.

[43] Yang B, Wang JQ, Tan Y, Yuan R, Chen ZS and Zou C. RNA methylation and cancer treatment. Pharmacol Res 2021; 174: 105937.

[44] Wang L, Hui H, Agrawal K, Kang Y, Li N, Tang R, Yuan J and Rana TM. m(6)A RNA methyltrans- ferases METTL3/14 regulate immune respons- es to anti-PD-1 therapy. EMBO J 2020; 39: e104514.

[45] Devarkar SC, Wang C, Miller MT, Ramanathan A, Jiang F, Khan AG, Patel SS and Marcotrigia- no J. Structural basis for m7G recognition and 2’-O-methyl discrimination in capped RNAs by the innate immune receptor RIG-I. Proc Natl Acad Sci U S A 2016; 113: 596-601.

[46] Farhood B, Najafi M and Mortezaee K. CD8(+) cytotoxic T lymphocytes in cancer immunother- apy: a review. J Cell Physiol 2019; 234: 8509- 8521.

[47] Cassetta L and Pollard JW. Targeting macro- phages: therapeutic approaches in cancer. Nat Rev Drug Discov 2018; 17: 887-904.

[48] DeNardo DG and Ruffell B. Macrophages as regulators of tumour immunity and immuno- therapy. Nat Rev Immunol 2019; 19: 369-382.

[49] Bindea G, Mlecnik B, Tosolini M, Kirilovsky A, Waldner M, Obenauf AC, Angell H, Fredriksen T, Lafontaine L, Berger A, Bruneval P, Fridman WH, Becker C, Pagès F, Speicher MR, Trajanos- ki Z and Galon J. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 2013; 39: 782-795.

[50] Crotty S. T follicular helper cell biology: a de- cade of discovery and diseases. Immunity 2019; 50: 1132-1148.

[51] Jardim DL, Goodman A, de Melo Gagliato D and Kurzrock R. The challenges of tumor muta- tional burden as an immunotherapy biomark- er. Cancer Cell 2021; 39: 154-173.

[52] Yarchoan M, Albacker LA, Hopkins AC, Monte- sion M, Murugesan K, Vithayathil TT, Zaidi N, Azad NS, Laheru DA, Frampton GM and Jaffee EM. PD-L1 expression and tumor mutational burden are independent biomarkers in most cancers. JCI Insight 2019; 4: e126908.

[53] Galanina N, Goodman AM, Cohen PR, Framp- ton GM and Kurzrock R. Successful treatment of HIV-associated Kaposi sarcoma with im- mune checkpoint blockade. Cancer Immunol Res 2018; 6: 1129-1135.

[54] Nghiem PT, Bhatia S, Lipson EJ, Kudchadkar RR, Miller NJ, Annamalai L, Berry S, Chartash EK, Daud A, Fling SP, Friedlander PA, Kluger HM, Kohrt HE, Lundgren L, Margolin K, Mitch- ell A, Olencki T, Pardoll DM, Reddy SA, Shantha EM, Sharfman WH, Sharon E, Shemanski LR, Shinohara MM, Sunshine JC, Taube JM, Thompson JA, Townson SM, Yearley JH, Topa- lian SL and Cheever MA. PD-1 blockade with

[55] pembrolizumab in advanced merkel-cell carci- noma. N Engl J Med 2016; 374: 2542-2552.

Orellana EA, Liu Q, Yankova E, Pirouz M, De Braekeleer E, Zhang W, Lim J, Aspris D, Send- inc E, Garyfallos DA, Gu M, Ali R, Gutierrez A, Mikutis S, Bernardes GJL, Fischer ES, Bradley A, Vassiliou GS, Slack FJ, Tzelepis K and Greg- ory RI. METTL1-mediated m(7)G modification of Arg-TCT tRNA drives oncogenic transforma- tion. Mol Cell 2021; 81: 3323-3338, e3314.

[56] Ying X, Liu B, Yuan Z, Huang Y, Chen C, Jiang X, Zhang H, Qi D, Yang S, Lin S, Luo J and Ji W. METTL1-m(7) G-EGFR/EFEMP1 axis promotes the bladder cancer development. Clin Transl Med 2021; 11: e675.

[57] Chang X, Liu X, Wang H, Yang X and Gu Y. Gly- colysis in the progression of pancreatic cancer. Am J Cancer Res 2022; 12: 861-872.

[58] Duan K, Gucer H, Kefeli M, Asa SL, Winer DA and Mete O. Immunohistochemical analysis of the metabolic phenotype of adrenal cortical carcinoma. Endocr Pathol 2020; 31: 231-238.

[59] Ciscato F, Ferrone L, Masgras I, Laquatra C and Rasola A. Hexokinase 2 in cancer: a prima Donna playing multiple characters. Int J Mol Sci 2021; 22: 4716.

[60] Shi L, Pan H, Liu Z, Xie J and Han W. Roles of PFKFB3 in cancer. Signal Transduct Target Ther 2017; 2: 17044.

m7G and METTL1 in ACC

Supplementary Table 1. The clinical characteristics of TCGA and GSE19750 cohorts
ItemsTCGA-ACCGSE19750
Sample size7922
Survival status
Dead2718
Alive524
AgeNA
<60/16
≥60/6
Clinical stage
Stage I91
Stage II377
Stage III161
Stage IV154
Unknown29
TNA
T19/
T242/
T38/
T418/
Unknown2/
MNA
M062/
M115/
Unknown2/
NNA
N068/
N19/
Unknown2/

NA, not available.

Supplementary Table 2. The m7G-related genes

Gene symbol

n=34

METTL1

WDR4

NSUN2

AGO2

CYFIP1

CYFIP2

DCPS

EIF3D

EIF4A1

EIF4E

EIF4E1B

EIF4E2

EIF4E3

EIF4G3

GEMIN5

IFIT5

LARP1

LSM1

NCBP1

NCBP2

NCBP2L

NCBP3

SNUPN

DCP2

NUDT1

NUDT10

NUDT11

NUDT16

NUDT16L1

NUDT3

NUDT4

NUDT4B

NUDT5

NUDT7

Supplementary Table 3. The specific sequences of sh-METTL1 and OE-METLL1
GeneSequence (5'-+3')
sh-METTL1 OE-METLL1CCGGGATGACCCAAAGGATAAGAAACTCGAGTTTCTTATCCTTTGGGTCATCTTTTTG
METTL1-XbaI-F: GCTCTAGAATGGCAGCCGAGACTCGGAACGTGGCCGG
METTL1-BamHI-R: CGGGATCCTCAGTGACCAGGCAGGCTGGTTTGGG

OE, over expression.

Supplementary Table 4. The detailed description of the gene sets used in metabolic analyses
NamesGene countsDescription
GO glycolytic process106Fermentation that includes the anaerobic conversion of glucose to pyruvate via the glycolytic pathway.
Hallmark Glycolysis200Genes encoding proteins involved in glycolysis and gluconeogenesis.
Biosynthetic process470The energy-requiring part of metabolism in which simpler substances are transformed into more complex ones, as in growth and other biosynthetic processes.
WP Nucleotide Metabolism19Nucleotide metabolism.
KEGG DNA Replication36DNA replication.
Amino acid and derivative Metabolic Process101The chemical reactions and pathways involving amino acids, organic acids containing one or more amino sub- stituents, and compounds derived from amino acids.
Hallmark Fatty acid Metabolism158Genes encoding proteins involved in metabolism of fatty acids.
Supplementary Table 5. The primer lists
GenePrimerSequence (5'-+3')
METLL1Forward5'-AGCTATACCCAGAGTTCTTCGCTCCAC-3'
Reverse5'-ACAGCCTATGTCTGCAAACTCCACT-3
TP53Forward5'-TAACAGTTCCTGCATGGGCGGC-3'
Reverse5'-AGGACAGGCACAAACACGCACC-3'
GAPDHForward5'-GTCGCCAGCCGAGCCACATC-3
Reverse5'-CCAGGCGCCCAATACGACCA-3'
Supplementary Figure 1. The risk plots of m7G risk signature. The top figure shows the m7G score of each ACC sample in TCGA cohort, which is arranged from low to high. The middle figure shows the shows the distribution of survival outcomes of each ACC sample. The bottom figure shows the expressive trends of four m7G signature genes in high- and low-m7G risk groups.

10

Risk score

8

6

12.5

Survival time

10.0

Status

7.5

· Alive

5.0

· Dead

2.5

Group

4

METTL1

2

NCBP1

0

NUDT1

-2

NUDT5

100%

Relative Percent

80%

60%

40%

20%

0%

Supplementary Figure 2. The infiltrating levels of 21 immune cells in each TCGA-ACC sample.

TCGA-OR-ASJA-01A-11R-A298-0

TCGA-OR-ASL4-01A-11R-A298-02

TCGA-OR-ASJS-01A-11R-A298-07

TCGA-OR-ASKW-01A-11R-A298-07

TCGA-OR-ASJE-01A-11R-A298-07

TCGA-OR-ASKX-01A-11R-A298-07

TOGA-OR-ASIT-01A-11R-A298-07

TCGA-P6-ASOF-01A-11R-A298-07

TCGA-OR-ASL9-01A-11R-A298-07

TCGA-OR-ASJL-01A-11R-A298-07

TCGA-OR-ASK9-01A-11R-A295-07

m7G and METTL1 in ACC

TCGA-OR-ASJV-01A-11R-A29S-07

TCGA-OR-ASKV-01A-11R-A298-07

TCGA-OR-ASJW-01A-11R-A29S-07

TCGA-OR-ASLD-01A-11R-A298-07

TCGA-OR-ASKT-01A-11R-A298-07

TCGA-PK-ASHA-01A-11R-A298-07

TCGA-OR-ASKO-01A-11R-A298-07

TCGA-OR-A5J7-01A-11R-A298-07

TCGA-PK-ASH8-01A-11R-A298-07

TCGA-OR-ASJP-01A-11R-A298-07

TCGA-OR-ASJC-01A-11R-A298-07

TCGA-OR-ASLJ-01A-11R-A298-07

TCGA-OR-ASLP-01A-11R-A298-07

TCGA-P6-ASOG-01A-22R-A298-07

TCGA-OR-ASLA-01A-11R-A298-07

TCGA-OR-ASKO-01A-11R-A298-07

TCGA-OR-ASLN-01A-11R-A298-07

TCGA-OR-ASLM-01A-11R-A298-07

TCGA-OR-ASK5-01A-11R-A298-07

TCGA-OR-ASK6-01A-11R-A298-07

TOGA-OR-ASK3-01A-11R-A298-07

TCGA-OR-ABLE-01A-11R-A298-07

TCGA-OR-A5LL-01A-11R-A298-07

TCGA-OR-A512-01A-11R-A298-07

TOGA-OR-ASJJ-01A-11R-A298-07

TCGA-OR-ASK8-01A-11R-A298-07

TCGA-PA-ASYG-01A-11R-A298-07

TOGA-OR-ASKU-01A-11R-A298-07

TCGA-OR-ASJO-01A-11R-A29S-07

TCGA-OR-ASLB-01A-11R-A298-07

TCGA-OR-ASK2-01A-11R-A298-07

TCGA-OR-A5J3-01A-11R-A298-07

TCGA-OR-ASL6-01A-11R-A298-07

TCGA-OR-ASJS-01A-11R-A298-07

TCGA-OR-ASJR-01A-11R-A298-07

TOGA-OR-ASJZ-01A-11R-A298-07

TCGA-OR-ASJO-01A-11R-A298-07

TCGA-OR-ASJO-01A-11R-A298-07

TCGA-OR-ASKY-01A-11R-A29S-07

TCGA-OR-ASKZ-01A-11R-A298-07

TCGA-OR-A5JK-01A-11R-A298-07

TCGA-OR-ASJY-01A-31R-A298-07

TCGA-OU-ASPI-01A-12R-A298-07

TCGA-OR-ASL8-01A-11R-A298-07

TOGA-OR-AS/1-01A-11R-A298-07

TCGA-OR-ASL3-01A-11R-A298-07

TOGA-OR-ASJM-01A-11R-A298-07

TCGA-OR-A5JX-01A-11R-A298-07

TOGA-OR-ASLS-01A-11R-A298-07

TOGA-PK-ASH9-01A-11R-A298-07

TCGA-OR-ASJI-01A-11R-A298-07

TCGA-PK-ASHB-01A-11R-A298-07

TCGA-OR-ASLG-01A-11R-A298-07

TCGA-OR-ASLT-01A-11R-A298-07

TOGA-OR-ASJG-01A-11R-A298-07

TCGA-OR-ASK4-01A-11R-A298-07

TCGA-OR-ASK1-01A-11R-A29S-07

Mast cells activated

Mast cells resting

Dendritic cells resting Dendritic cells activated

Macrophages M2

Macrophages M1

Macrophages MO

Monocytes

NK cells activated

NK cells resting

= T cells gamma delta

= T cells follicular helper T cells regulatory (Tregs)

= T cells CD4 memory resting T cells CD4 memory activated

T cells CD4 naive

T cells CD8

Plasma cells

B cells naive B cells memory

TCGA-OR-ASLO-01A-11R-A298-07

TCGA-OR-ASJB-01A-11R-A298-07

TCGA-OR-ASLK-01A-11R-A298-07

TCGA-OR-ASJ9-01A-11R-A298-07

TOGA-OR-ASLR-01A-11R-A298-07

TCGA-OR-ASJ6-01A-31R-A298-07

TCGA-OR-ASJF-01A-11R-A298-07

TCGA-OR-A5J8-01A-11R-A298-07

TCGA-OR-ASLC-01A-11R-A298-07

Neutrophils

Eosinophils

TCGA-OR-ASL5-01A-11R-A298-07

TCGA-OR-ASLH-01A-11R-A298-07

m7G and METTL1 in ACC

Supplementary Figure 3. Immune heatmap of each m7G-risk group. The bar on the right of figure shows the differ- ent immune analysis items. Different colors indicate the level of immune score.

TumorPurity

TumorPurity

ESTIMATEScore

3

1

ImmuneScore

StromalScore

2

Subtype

0.5

aDCs

1

Type-II IFN Response

ESTIMATEScore

iDCs

0

3000

NK cells

Mast-cells

-1

Tfh

CD8+ T cells

-2

-2000

Inflammation-promoting

Cytolytic activity TIL

-3

ImmuneScore

2000

Check-point

T cell co-inhibition

T cell co-stimulation Th1 cells

-1000

APC co-inhibition

StromalScore

Treg

1500

CCR

Parainflammation

B cells

-1500

Neutrophils

pDCs

Subtype

T helper cells

Low

APC co-stimulation

High

HLA

Macrophages

Th2 cells

Type-I IFN Response

DCs

MHC-class I