ERBB3 methylation and immune infiltration in tumor microenvironment of cervical cancer

ERBB3, a member of the ERBB family of receptor tyrosine kinases, plays an important role in cancer, despite its lack of intrinsic carcinogenic mechanism of cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC). Research on bioinformatics methods through multi-omics, this work proves that ERBB3 gene mutation, methylation modification have extensive regulatory mechanisms on the CESC microenvironment. We found that ERBB3 is involved in carcinogenesis of cervical cancer and is not associated with its prognosis. The carcinogenic mechanism is mainly related to the suppression of the immune system between tumor infiltrating lymphocytes (TILs) and the methylation of the RNA level. Our study indicated ERBB3 is more likely to be a carcinogenic factor than a key prognostic factor for cervical cancer. Methylation of ERBB3 may work as a checkpoint immunotherapy target in CESC, DNA methylation modification of the 4480 base pair downstream of ERBB3 transcription initiation site was the highest.

The Human Protein Atlas Gene set enrichment analysis database. The "Human Proteome" chapters provide a knowledge-based analysis and entry into the Human Protein Atlas (https:// www. prote inatl as. org/) from different defined transactions of the human tissue proteome. Analysis of four genes RNA-seq expressed in 291 cervical cancer samples, The RNA-seq data is reported as average FPKM (number Fragments Per Kilobase of exon per Million reads).
GO and KEGG analysis. Using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) 35 (https:// david. ncifc rf. gov/), Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis based on co-expressed genes (https:// www. kegg. jp/), p < 0.05. Survival prognosis analysis. We used the "survival map" module of gepia2 to obtain the Overall Survival (OS) and Disease-Free Survival (DFS) significance map data of four genes in ERBB family of CESC in all TCGA. To assess correlations between gene expression level and OS/DFS, patients were classified into low and high mRNA expression subgroups using median expression (50%) as the cut-off. The hypothesis was tested by log rank test, and the survival map was obtained through the "survival analysis" module of gepia2. The correlation between overall survival time and m6A regulators expression is measured using Kaplan-Meier Plotter Database (https:// kmplot. com/ analy sis/) 36 .

Correlation between molecular and methylation degree of DNA methylation sites. Data from
TCGA RNA-seq data in Level 3 HTSeq-FPKM format in CESC (cervical squamous cell carcinoma and adenocarcinoma) project and methylation data from illumina human methylation 450. RNAseq data in FPKM (fragments per kilobase per million) format is converted into TPM (transcripts per million reads) format and log2 conversion is performed.
Statistical data. Statistical analysis on the TCGA data was performed with R software (version 3.6.0). Differentially expressed genes (DEGs) were identified using R package of mainly ggplot2 version 3.3.3. The Cluster-Profiler package was used for enrichment analysis and visualization. In the previous analysis, except specifically mentioned, p value less than 0.05 was considered statistically significant.

Results
Gene differential expression and clinical significance of ERBB3 in CESC. We first compared the expression pattern of EGFR, ERBB2, ERBB3, ERBB4 in CESC tumor and normal tissues. ERBB3 obtain a significant difference in CESC (Fig. 1A). The mRNA levels of different primary cervical cancer cell lines indicate that ERBB3 is highly expressed in cervical malignant cell lines dominated by adenocarcinoma: AV3, Hela (Fig. 1B). Further research on pathological types and HPV typing, we found that in the three types of cells with higher expression of ERBB3, there was no statistical difference between SiHa and HeLa, while the statistical difference between SiHa and C33A significantly (p < 0.0001). The difference between the two cell lines is mainly HPV infection. The low expression of ERBB3 in HPV(-)-C33A inspired us to study the significance of ERBB3 and HPV in the carcinogenicity of cervical cancer in the future.
At the same time, GEPIA2 provided the structural map of the isomeric protein domain 37 based on Pfam prediction to show the structural differences among the isomers (Fig. 1C). Figure 1D-F analyze the clinical significance of the ERBB family. The DFS Prognostic (Fig. 1D) and OS rate (Fig. 1E) of EGFR, ERBB2, ERBB3 and ERBB4 have no significance (p > 0.05). Figure 1G shows the accordion diagram of each gene in cervical cancer staging. We found that the difference in the expression of the four genes between different stages of cervical cancer were not significant. Figure 1F and Table 1 list the RNA-seq data of ERBB Family in 291 cervical cancer samples. RNA-seq data of EGFR, ERBB2, ERBB3, and ERBB4 in cervical and ectocervix tissue type reported as average FPKM (number Fragments Per Kilobase of exon per Million reads), generated by the The Cancer Genome Atlas (TCGA).
According to the TCGA database in Fig. 2B, the alteration frequency of ERBB3 in adenocarcinoma is higher than squamous cell carcinoma. This is consistent with ERBB3 mRNA expression in different cell lines (Fig. 1B).
We showed the mutation site of V104ML with the case ID list as TCGA-C5-A1MI-01, TCGA-C5-A7CO-01, and TCGA-EA-A5FO-01 in receptor domain, which has the highest change frequency in the three-dimensional structure of ERBB3 (Fig. 2C,D). The potential correlation between mutation status and progression-free, overall, disease-specific and disease-free survival of CESC is not significant (Fig. 2E). m6A regulators is differently expressed in CESC cancer. We analyzed the expression of 22 major m6A RNA methylation modulators in 607 CESC patients from the TCGA dataset. This study showed that mutations in the m6A regulator of the human CESC genome were associated with pathogenesis. There were nine m6A regulatory factors significantly increased in CESCs: HNRNPC, YTHDF1, HNRNPA2B1, IGF2BP1, VIRMA,  IGF2BP2, YTHDF2, RBM15, and IGF2BP3. There were seven m6A regulators with significantly reduced expression: ZCCHC4, METTL3, ZC3H13, YTHDC1, YTHDC2, METTL16, and FTO. This study focuses on the role of m6A in CESC. Detection of genetic variation of m6A regulators in 607 patients using cBioPortal database (Fig. 3), among which IGF2BP2 displayed the highest incidence rate (17%). www.nature.com/scientificreports/ By comparing the expression level of 22 m6A regulator genes in tumors with adjacent normal tissues, these seven genes were down-expression in tumor tissues: ZCCH4 (Fig. 4A) 4V) were up-regulated in tumor tissues. Yellow background indicated the most significant (p < 0.0001) genes up-regulated in tumor tissues, including HNRNPC, YTHDF1, IGF2BP1, YTHDF2, RBM15, and IGF2BP3. Blue background indicated the most significant genes down-regulated in tumor tissues, including YTHDC1, YTHDC2, METTL16, and FTO. Through further analysis, the absolute expression of HNRNPC in tumor tissues was the highest, and the absolute expression of YTHDC2 was the lowest.
In Fig. 5, we obtain the prognostic value of 22 m6A regulators in CESC through the Kaplan-Meier plot database. We revealed higher levels of ZC3H13, WTAP, HNRNPC, YTHDF3, and VIRMA were significantly associated with worse outcomes in CESC (HR > 1, blue background), while YTHDC1, YTHDF1 were protective factors of cervical cancer (HR < 1, orange background). It indicates these m6A regulators had key roles in CESC prognosis.
Relationship between ERBB3 and m6A regulators. The upper part of Fig. 6 is the expression of ERBB3 gene, the lower part is the expression of m6A regulator genes after Z-score transformation. The thermal maps show ERBB3 expression has the most significant relation with YTHDC1, METTL14, RBM15, RBM15B, CBLL1, ZC3H13, METTL3, YTHDC2, and ZCCHC4. ERBB3 had the highest correlation with YTHDC1.
For the sake of investigating the downstream pathways of hub m6A regulators in CESC, we performed GO and KEGG analysis using co-expression genes of 6 m6A regulators. The results showed that YTHDC1 (Fig. 7A) and YTHDC2 (Fig. 7B) belong to RNA binding protein families. Both YTHDC1 and YTHDC2 can recognize and bind RNA containing N6 methyladenosine (m6A); The YTH domain mediates this binding [38][39][40] . M6A is a modifier existing in mRNA and some non-coding RNA internal sites, and plays a role in mRNA splicing, processing and stability regulation. YTHDC1 (also known as splicing factor yt521) can be used as a key regulator of exon addition or exon skipping to regulate selective splicing. YTHDC1 can promote exon increase by recruiting srsf3 into the region containing m6A, and inhibit exon skipping by blocking srsf10 binding to these same regions 20 . RBM15 (Fig. 7C,D) plays a major role in RNA modification and mRNA metabolic regulation.  (Fig. 8H). Beta value is a data format for estimating the degree of methylation. It represents the ratio of signal intensity between methylated and unmethylated bases. Beta values are usually between 0 and 1, 0 for no methylation and 1 for complete methylation. The probe of cg11835619 (TSS + 4480) has the highest relation with ERBB3 DNA methylation (r = 0.62, p < 0.001). Correlation analysis only shows that a methylation site regulates the correlation degree of ERBB3 gene, but further experiments (artificially regulating this methylation site to detect whether the expression of this gene has changed) must be done to clarify the causality.
ERBB3 methylation andtumor immunity in the microenvironment (TIME) classification. There are interactive regulatory mechanisms in the local immune microenvironment of cervical cancer tumors. Cells, related cytokines, and immune effector cells have different distribution patterns and infiltration densities. We studied the relationship between ERBB3 DNA methylation modification and the infiltration of various immune components in the tumor immune microenvironment (TILS, Immunomodulator and Chemokine), and explored the possible mechanism of ERBB3 shaping the immune environment of cervical cancer based on DNA methylation.

Discussion
ERBB receptor is a typical cell membrane receptor tyrosine kinase, which is activated by dimerization after binding to ligands. The ERBB receptor tyrosine kinase family contains four cell surface receptors: ERBB1/EGFR/ HER1, ERBB2/HER2, ERBB3/HER3 and ERBB4/HER4. HER3/ERBB3 is a member of the ERBB receptor protein tyrosine kinase family, but lacks tyrosine kinase activity. The tyrosine phosphorylation of ERBB3 depends on its binding to other ERBB tyrosine kinases. When binding to ligands, heterodimers were formed between ERBB3 and other ERBB proteins, and ERBB3 was phosphorylated by activated ERBB kinase on tyrosine residues 41,42 .
There are at least nine potential tyrosine phosphorylation sites in the tail region of carboxyl end of ERBB3. These sites are the common binding sites of signal transduction proteins (including Src family members, Grb2 and PI3 kinase p85 subunit), which can mediate downstream signal transduction of ERBB 43 . The Tyr1222 and Tyr1289 sites of ERBB3 are located in the YXXM motif and participate in PI3K signal transduction 44 . Researchers have found that ERBB3 is highly expressed in many cancer cells 45 . www.nature.com/scientificreports/ In this study, it was found that in ERBB family, the expression of ERBB3 in cervical cancer tissues was higher than that in normal tissues (Fig. 1A). It shows that ERBB3 may bes closely related to adenocarcinoma and HPV positive cervical carcinoma (Figs. 1B, 2B). Figure 1F shows the comparison of transcriptome expression of four ERBB families in cervical cancer tissues in Table 1. It is found that the high expression of ERBB3 in tumor tissues is consistent with the high expression of ERBB2. However, ERBB3 had no statistical significance in clinical disease staging and disease prognosis (Figs. 1D,E,G, 2E). Therefore, we speculate that ERBB3 is a pathogenic factor of cervical cancer rather than a prognostic factor. The hotspot mutation site of ERBB3 in cervical cancer is extracellular domain V104M/L ( Fig. 2A,C), the one has been shown to be a statistically significant mutation hotspot to promote oncogenic signalling 46 . Treatments of cells harboring V104M mutation in ERBB3 with ERBB antibodies and other inhibitors blocked oncogenic signaling 47 .
Abnormal methylation is a prominent feature of cancer. It is unclear how DNA methylation affects immune surveillance and tumor metastasis. N6-methyladenosine (m6A) is one of the most common and representative chemical modifications in eukaryotic RNA. It is a kind of dynamic and heritable information, which is widely present in a variety of organisms. M6A is mainly divided into its important roles in regulating gene expression, splicing, RNA editing, regulating RNA stability, and controlling mRNA degradation. It is a reversible epigenetic modification. Mainly divided into three categories: m6A methyltransferase (writer), m6A demethylase (eraser), m6A binding protein (reader).
Among them, we found higher expressions of HNRNPC, YTHDF1, IGF2BP1, YTHDF2, RBM15, and IGF2BP3, lower expression of YTHDC1, YTHDC2, METTL16, and FTO in CESC (Fig. 4). YTHDC1 and YTHDF1 have anti-cancer effects in cervical cancer, while ZC3H13, WTAP, HNRNPC, YTHDF3 and VIRMA have tumor-promoting effects in cervical cancer (Fig. 5). We compared the methylation regulatory factors with the highest correlation in Figs. 5 and 6, and obtained three factors in Fig. 7. Through further GO and KEGG clustering analysis of the functions of these three factors (YTHDC1, YTHDC2 and RBM15), we infer that ERBB3 methylation plays an important role in the occurrence and progression of cervical cancer.
According to the analysis of ERBB3 DNA methylation in Fig. 8, in cervical cancer, the correlation of DNA methylation modification of the 4480 base pair downstream of ERBB3 transcription initiation site was the highest, but whether the gene was finally regulated still needs further experiments to determine the causal relationship. www.nature.com/scientificreports/ GO and KEGG analysis confirm ERBB3 Methylation were involved in regulating RNA splicing, RNA stability, and cell proliferation create Tumor and immune system interaction, which integrates multiple heterogeneous data types (Fig. 7).
We researched the relationship between ERBB3 methylation and immune cell infiltration in cervical cancer microenvironment, finding that ( Table 2) the abundance of TH1, MDSC, Macrophage, effector memory CD8 T cell, activated CD8 T cell, immature B cell and regulatory T cell have the significant association with methylation of ERBB3 in cervical tumor immune microenvironment (R > 0.6). Relations between three kinds of immunomodulators and methylation of ERBB3 are as follows: TIGIT (r = 0.656), ICOS (r = 0.702), and HLA-E (r = 0.514). TIGIT (T cell Ig and ITIM domain) is a member of the poliovirus receptor (PVR)/Nectin family 48 . It is expressed in lymphocytes, especially effector and regulatory CD4 + T cells, follicular auxiliary CD4 + T cells, effector CD8 + T cells and natural killer (NK) cells. TIGIT plays an inhibitory role in multiple steps of the tumor immune cycle 49 . The immune regulation of ICOS is manifested in the following aspects: enhancing the pattern recognition receptor signal of dendritic cells, inducing CD4 + T cells to produce IL-10 and producing high affinity antibodies against specific antigens. CXCR5 (r = 0.671) is a G protein-coupled seven transmembrane receptor, belonging to the CXC chemokine receptor family 50 , and its ligand is the chemokine CXCL13 (r = 0.594). CCR5 (r = 0.656) is the receptor of intracellular β-chemokines (RANTES, MIP1α and MIP1β), which has the function of regulating the migration, proliferation and immunity of T cells and monocytes/macrophage cell lines 51 . It is mainly expressed in the memory quiescent T Lymphocytes, monocytes, immature dendritic cells, etc. on the cell membrane. When cancer cells spread in the body, secondary tumors called metastases can form. These secondary tumors cause approximately 90% of cancer patients' deaths. An important way to spread cancer cells is through the lymphatic system, which runs through the entire body like the vascular system and connects the lymph nodes to each other. When white blood cells migrate through the lymphatic system to coordinate defenses against pathogens, as a specific These results demonstrate that ERBB3 dynamically reshapes the composition and function of the immune macroenvironment in cervical cancer. The correlations between methylation and TIME shows potential correlation between methylation of ERBB3 and inhibitory immune checkpoints along with immunosuppressive Tregs, MDSCS and potentially macrophages indicates an immune permissive environment favoring tumors to grow. Chemokines CXCL19 and chemokine receptors CXCR5 mediate immune cell trafficking into the tumour micro environment based on the DNA methylation of ERBB3 in cervical cancer.

Conclusions
In summary, this work proves that ERBB3 gene mutation, methylation modification have extensive regulatory mechanisms on the CESC microenvironment. ERBB3 is more likely to be an important carcinogenic factor of cervical cancer, but it has no significant effect on the clinical stage and prognosis of the disease. The differences in methylation modification patterns lead to the heterogeneity and complexity of CESC tumor immune microenvironment. Our study found that m6A regulator was an important biomarker of CESC and was closely related to tumor immune infiltration in cervical cancer. A comprehensive assessment of the ERBB3 multi-omics will help to enhance our understanding of the characteristics of cell infiltration in CESC tumor microenvironment and guide more effective immunotherapy strategies. DNA methylation modification of the 4480 base pair downstream of ERBB3 transcription initiation site was the highest. Further experiments will verify the carcinogenic mechanism of this methylation site on ERBB3 in cervical cancer.