ERBB3 mediates the PI3K/AKT/mTOR pathway to alter the epithelial‑mesenchymal transition in cervical cancer and predict immunity filtration outcome

Cervical cancer is the fourth most common cancer among women worldwide, and the prognosis of advanced/recurrent cervical cancer remains poor. Metastasis and invasion of this type of cancer are closely associated with the tumor microenvironment. Studying the complex interactions between tumor progression and immune cells or stromal cells can provide new insights into treatment for patients with aggressive tumor, recurrence and drug resistance. In the present study, a bioinformatics method (Gene Expression Profiling Interactive Analysis, differentially expressed genes, Gene Ontology, Kyoto Encyclopedia of Genes and Genomes, protein-protein interactions and survival analysis) was used to explore the mRNA and protein level discrepancy gene signature of ERBB3 via the PI3K/AKT/mTOR pathway from the speculation in immuno-oncology and experimental verification of different cervical cancer cell lines. The high expression of ERBB3 in cervical cancer tissues (especially HPV-positive and adenocarcinoma-related) promoted the activation of the PI3K/AKT/mTOR pathway. The increased expression of MMP9 changed the macrophage infiltration in the tumor microenvironment and affected prognosis of patients with cervical cancer. In conclusion, the present study identified 14 EMT-related genes and 30 genes involved in the PI3K/AKT/mTOR pathway in cervical cancer, and they might provide novel clues for future treatment. The MMP family may be a notable factor associated with tumor cells and immune cells.


Introduction
Cervical cancer is the most common malignant tumor of the female genital tract worldwide (1). Persistent high-risk human papillomavirus (HPV) oncoproteins E6 and E7 are the main pathogenic factors of cervical cancer (2). At the molecular level, the HPV E6 and E7 proteins directly activate Akt, and this pathway is further stimulated in cervical cancer cells by amplifications and mutations of the PI3K genes. As it has a key role in the control of HPV gene expression and development of cervical cancer, the PI3K/AKT/mammalian target of rapamycin (mTOR) pathway may have potential as a therapeutic target for cervical cancer (3)(4)(5).
The tumor microenvironment (TME) contains stromal cells and immune cells that shape cancer development and impact responses to tumor therapy (6). Cancer cell proliferation, angiogenesis and metastasis also contribute to the establishment of an immunosuppressive environment. These factors are associated with tumor progression and poor clinical outcomes (7,8). However, factors that contribute to immunosuppression in the TME are poorly defined.
Epithelial-mesenchymal transition (EMT) plays an important role in tumor development from initiation to metastasis. EMT contributes to the majority of the hallmarks of cancer and continues to be an attractive target for cancer therapy (9). Classical EMT is characterized by the phenotype change of epithelial cells to cells with mesenchymal properties, but EMT is also associated with multiple other molecular processes, including tumor immune evasion (10). Immunosuppression occurs as a direct consequence of the EMT program or develops through some additional, still-uncharacterized signaling channels (11).
The PI3K/AKT/mTOR pathway can promote migration and induce EMT in numerous types of tumors, including cervical cancer (8,12). EMT-related changes in the expression of various receptor tyrosine kinases (RTKs) (13) have been reported, although the proliferation and survival dependence of specific RTKs under different conditions remains to be fully elucidated. The expression of an EGFR family member-ERBB3 (14) is associated with the epithelial phenotype of the cell line and the sensitivity to EGFR inhibition. ERBB3 heterodimerizes (15) with additional EGFR family members after stimulation with ERBB3 mediates the PI3K/AKT/mTOR pathway to alter the epithelial-mesenchymal transition in cervical cancer and predict immunity filtration outcome various ligands, including neuregulin (16). ERBB3 contains multiple binding sites for p85, which is the regulatory subunit of PI3K (17). This allows direct recruitment and activation of PI3K signals via ERBB3 (18). Although changes in ERBB3 expression have been observed, the functional consequences of these changes and the relationship with downstream signals after EMT (19) have not been fully described. These targets involved in cervical cancer are not functionally exclusive; rather, they are intertwined and reciprocal, and together they form intricate TME networks to meet context-specific needs for cellular function. To improve understanding of the correlation between TME and prognosis of cervical cancer, the present study assessed cervical cancer cell lines and tissues to reveal the roles of ERBB3 in the EMT induction of TME harboring immunosuppression, migration and invasion of cervical cancer, and to explore whether PI3K/AKT/mTOR signaling is involved in this process.
The current study intended to explore how ERBB3 mediates the PI3K/AKT/mTOR pathway and changes the tumor immune microenvironment to affect the EMT status of cervical cancer, which may provide further understanding of MMPs involved in immunotherapy.

Materials and methods
Data collection and preprocessing. RNAseq data in the transcripts per million (TPM) format from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/) and The Genotype-Tissue Expression (GTEx) (https://gtexportal. org/) were uniformly processed using the Toil process (20). Through extraction of the cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) data in TCGA and corresponding normal tissue data in GTEx, the 306 cervical tumor samples were classified as the malignant group and the three samples adjacent to cancer from TCGA together with the 10 normal cervix tissues from GTEx were classified as the non-malignant group. The RNAseq data in TPM format were log2 transformed for expression comparison between samples.
Related hub genes selection. A total of 14 EMT-related genes were selected based on the article 'Guidelines and definitions for research on epithelial-mesenchymal transition' written by the EMT International Association (TEMTIA) in 2020 (21). When searching for pathway related genes, GeneCards (https://www.genecards.org/) was used; the key words 'PI3K/AKT' and 'PI3K/AKT/mTOR' were searched for and genes with a relevance score >4.0 were selected.
Gene Expression Profiling Interactive Analysis (GEPIA). GEPIA (http://gepia.cancer-pku.cn/index.html) is a user-friendly web portal for gene expression analysis based on TCGA and GTEx data. In the current study, expression analysis of ERBB3 was evaluated using the project ID of TCGA-CESC. In the module 'Expression DIY' of GEPIA, the expression of ERBB3 between pan-cancer and normal tissue samples was studied with the option of matching normal TCGA data and GTEx data (Fig. 1A).
Differentially-expressed genes (DEGs). Batch correction, normalization and difference analysis of RNA-seq data from GSE63514, GSE9750, and GSE44001 were performed to screen for DEGs in CESC samples. GSE63514 (22), GSE9750 (23) and GSE44001 (24,25) were obtained from the NCBI Gene Expression Omnibus database (ncbi.nlm.nih.gov/geo). The GSE63514 dataset used the GPL570 [HG-U133 Plus 2] Affymetrix Human Genome U133 Plus 2.0 Array platform, which contained 28 cervical cancer tissue samples and 24 normal samples. The GSE9750 dataset used the GPL96 Affymetrix Human Genome U133A Array platform and included 33 cervical cancer tissue samples that were primarily marked by HPV16 or HPV18, and 21 normal cervical samples. The GSE44001 dataset used the GPL14951 Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip, which contained 300 cervical cancer tissue samples. GSE63514 and GSE9750 were set as the reference group and GSE44001 as the test group, and the R software limma package was used to identify DEGs between the groups (26). A total of 13,473 DEGs, including 6,514 downregulated and 6,595 upregulated genes, were identified in cervical cancer. The results were visualized using R software (version 3.6.3) (statistical analysis and visualization) with the R package ggplot2 [version 3.3.3] (27) to generate a volcano plot (Fig. 1C), which identified important genes.

Genomic alteration types and alteration frequency analysis.
Genomic alteration types (missense mutation with putative driver or unknown significance, amplification and no alterations) and alteration frequency of 14 EMT-associated genes and 30 PI3K/AKT/mTOR pathway-associated genes were obtained from the cBioPortal for Cancer Genomics (http://www. cbioportal.org), using the 'OncoPrint' module and 'Cancer Types Summary' module for visualization.
Immune cell infiltration estimation. For the immune infiltration analysis, transcriptome or other omics data was used to calculate the score of immune cells in the tissue through algorithms, and inferred the infiltration of immune cells in the tissue. Single-sample Gene Set Enrichment Analysis (GSEA) in immune infiltration, which uses the markers of each type of immune cells (28), was used as the gene set to calculate the enrichment of each type of immune cells in each sample.

Gene Ontology (GO) Term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and
GSEA. GO and KEGG (29,30) analyses were applied to explore the biological functions of target genes in CESC. GO analysis is a powerful bioinformatics tool to determine the biological processes (BPs), cellular components (CCs) and molecular functions (MFs) related to ERBB3. GSEA was utilized to investigate the potential mechanisms of ERBB3. GO, KEGG and GSEA were performed using the R (version 3.6.4) package ClusterProfiler (31). P<0.1 and q<0.2 were selected as the cut-off level.
Protein interactions and biological processes. The direct and indirect relationship between ERBB3 and 30 hub genes in the PI3K/AKT/mTOR signaling pathway were analyzed using the online tool STRING (https://string-db.org).
Risk survival analysis. Kaplan-Meier curves (33) can describe the survival status of each group of patients or the survival status of each group of experimental animals. The present study analyzed RNA-seq data from TCGA-CESC cohort (n=304) and selected the median as the cutoff value. Hazard ratio (HR) is defined as the ratio of the two risk rates. When HR is >1, the research object is a risk factor; when HR is <1, the research object is a protective factor; when HR=1, the research object has no effect on survival time. As an outcome index, overall survival (OS) refers to the time to death. The prognostic data come from a Cell article (34). Data filtering: control/normal (not all items have control/normal) were removed and clinical information was retained. For the. nomogram chart, on the basis of multifactor regression analysis, the ruler score was set to characterize the situation of each variable in the multifactor regression model, and finally the total score was calculated to predict the probability of event occurrence (35).
Statistical analysis. Software R (version 3.6.3) (37) was used for statistical analysis and visualization. For differential analysis of single gene expression, the R package of ggplot2 (version 3.3.3) (38) was used for visualization. For multigene association analysis, we used the R package of igraph (version 1.2.6) (39) and ggraph package (version 2.0.5) (40). For GO-KEGG analysis and GSEA, the R package of ggplot2 and cluster Profiler package was used (33). Visualization of Kaplan-Meier OS analysis was based on the use of the R package of survminer (0.4.9 version) (41) and for statistical analysis of survival data the survival package (3.2-10 version) was used. Wilcoxon rank sum test was used to assess differences in gene expression. Spearman's rank correlation coefficient was used to assess the significance of correlations. The qPCR data are presented as the mean ± standard deviation of three experiments, and qPCR and RNA-seq data were analyzed using one-way ANOVA followed by Tukey's multiple comparison test. P<0.05 was considered to indicate a statistically significant difference.
The total number of gene identifications (IDs) after removing the null value was 35,905. Among them, 1,706 IDs met the |log2(FC)|>2 and P<0.05 threshold. Under this threshold, there were 1,221 IDs with high expression (positive logFC) and 485 with low expression (negative logFC). The genes that met this threshold and were significantly associated with EMT included MMP3 and SNAI2 (downregulated genes), and KRT12 (upregulated gene) (Fig. 1C). The expression of seven hub genes in the PI3K/AKT/mTOR pathway were increased in cervical cancer: EIF4EBP1, GSK3B, HRAS, KRAS, NRAS, PIK3CB and PIK3R2 (Fig. 2).
Figs. 3 and 4 show the somatic variation pattern in cervical cancer. These schematics represented the distribution of the number of protein-altering somatic mutations and copy number variations in 607 samples of cervical cancer. The highest frequencies of mutations among the EMT-related genes were revealed in MMP3 (11%, including amplification, deep deletion and missense mutation); PIK3CA (37%, including amplification and missense mutation); and PTEN (11%, including deep deletion and truncating mutation). PIK3CA mutation status assesses the hotspot mutations of the PIK3CA gene (42). . Therefore, it is of certain significance to study this pathway in relation to carcinogenesis and prognosis of cervical cancer (Fig. 2B-D).

Selection of EMT-related genes with GO and KEGG analyses.
It was revealed that in cervical cancer the ERBB3 gene was enriched in the PI3K/AKT signaling pathway (NES=-1.707, P=0.045, FDR=0.034; Fig. 5).
TWIST1, tight junction protein 1 (TJP1), MMP9, MMP3 and vimentin (VIM) (Fig. 6) with all annotated functional molecules were compared using hypergeometric distribution tests to determine which functional roles were involved in that stack. The function of genes was divided into three categories: BPs, CCs and MFs.
Proteomics. By comparing with the RNA level in normal cervical epithelium (Fig. 9A), it was revealed that the expression level of MMP9 in cervical cancer was significantly different. Through IHC analysis of the mRNA-protein expression of the EMT-related genes in cervical cancer tissues, it was revealed that the expression level of MMP family was relatively higher compared with those of other EMT-related genes (Fig. 9B). EMT of cervical cancer is associated with upregulation of MMP9. According to our previous findings of ERBB3 sequencing of cervical cancer cell lines (43), the mRNA levels of different primary cervical cancer cell lines indicated that The upright triangles represent the highly expressed genes, the inverted triangles represent the poorly expressed genes and the circles represent genes that have no significantly different expression in cervical cancer. (C) Volcano map of differentially-expressed genes, where red represents the upregulated genes and light blue represents the downregulated genes. Seq, sequencing; ERBB3; TPM, transcripts per million reads; ns, not significant; CESC, cervical squamous cell carcinoma and adenocarcinoma.
ERBB3 is highly expressed in cervical malignant cell lines dominated by SiHa and HeLa (Fig 9C). After further research on pathological types and HPV typing, it was revealed that in the three types of cells with higher expression of ERBB3, there was no significant difference between SiHa and HeLa, while there was a significant difference between SiHa and C33A (P<0.0001). These data show that ERBB3 is closely associated with adenocarcinoma and HPV-positive cervical carcinoma.     Among the abovementioned eight types of immune cell infiltration, the prognosis of the group with a higher ERBB3 mRNA level was decreased (Fig. 11).

EMT-related genes change the immune cell infiltration.
If the absolute value of R is below 0.3, there is no straight phase off; R of ≥0.3 denotes linear correlations; R values between 0.3 and 0.5 indicate low-degree correlations; R values between 0.5 and 0.8 refer to significant correlations; and R values of 0.8 and above are high-degree correlations. As shown in Table II, correlation between the EMT-related factors and immune cell infiltration was not high. Among the EMT-related factors, MMP9, MMP2, and ZEB1 were closely associated with the immune system. The EMT status (Fig. 12) may be related to MMP9 changing the tumor immune microenvironment through dendritic cells and macrophages (10).
Risk score was constructed using eight selected genes through the multifactor analysis of the disease prognosis model, and three prognostic types of OS, disease-specific survival and progression-free interval (PFI) were analyzed. ERBB3, CD47, MMP9, TWIST1, CDH2, PTEN, VIM and ZEB1 were included. Multivariate analysis revealed that CDH2, MMP9, and VIM were significant factors in the assessment of PFI. Gene signatures of cervical cancer cell immune-oncology microenvironment positively correlated with the patients' survival. These analyses (Fig. 14) indicated that the selected genes and constructed risk prognostic models had good prognostic value.

Discussion
High-risk HPV16 DNA is integrated into the host cell genome (HPV16: q21-q31 of chromosome 2.13; HPV18: chromosome 24.8) (44), disrupts the open reading frame, and causes overexpression of E6 and E7 genes (45). It has been demonstrated that E6 and E7 exert carcinogenic effects by combining with cell cycle regulators, such as p53 (a transcription factor related to the PI3K pathway as shown in Table I) and retinoblastoma (46). E6 can interact with E6-related protein E6AP to form a complex and bind to p53 (47), hydrolyze p53, and cause the loss of negative regulation of cell proliferation induced by p53, thereby leading to uncontrolled cell proliferation and malignant transformation. The present study also revealed that ERBB3 was highly expressed in HPV-infected cell lines and was associated with adenocarcinoma. Therefore, it can be speculated that HPV-positive cervical cancer cells and adenocarcinoma are the carcinogenic factors or prognostic factors of cervical cancer.
To investigate the association between the actual activation of the PI3K pathway and immune infiltration, the DEGs of ERBB3 in CESC was assessed (GSE63514, GSE9750 and GSE44001 datasets from the Gene Expression Omnibus database were analyzed) for hub genes of the PI3K/AKT/mTOR pathway and cancer progression. Pathway enrichment, protein-protein interaction and pathway crosstalk analyses were performed to identify key genes and pathways. The current study illustrated that cancer-immune interactions might differ depending on specific alterations in the PI3K pathway, demonstrating that genetic aberrations in malignant cells influence the immune landscape of tumors.
The diversity of EMT creates a wide range of heterogeneity in tumors, and may provide tumor cells with increased adaptability and resistance, enabling them to survive and proliferate in a complex TME, and metastasize and invade lymph and blood vessels. The present study demonstrated that MMPs, especially MMP9 as a prominent representative, are highly relevant for TME and immune cells. MMPs are a family of zinc-dependent endopeptidases (48,49). The biological function of MMPs is to degrade various molecules used for cell adhesion and regulate the interaction between cells and the extracellular matrix. Recent studies (50)(51)(52) have shown that MMPs are highly associated with the microenvironment of tumors and immune cells, and targeting MMPs may overcome the barriers of immunosuppression. However, the present study revealed that the expression of MMP9 was not a significant predictor of OS in patients with cervical cancer. MMP9 was significantly associated with ERBB3.
In the complex TME, the same anti-infection immune cells can be destroyed by tumor cells (53). As a result, the antitumor immune cells not only do not destroy the transformed cells, but they even change to immune cells that promote tumor growth and metastasis (54,55). These immune cells secrete factors that promote survival, promote migration, and resist detection. Hence, the present study discussed the mechanism of accelerating cervical cancer tumor progression from the perspective of EMT-associated immune evasion.
Cellular immunity is necessary for clearing HPV-infected and HPV-transformed tumor cells. HPV-specific CD8 cytotoxic T lymphocytes (CTLs) are needed for the immune defense against cervical cancer. However, the function of CTLs may be blunted by systemic and local immunosuppressive environments associated with tumor growth (56,57). A series of clinical trials (58)(59)(60) have shown that the immune system is unable to completely eradicate the tumor despite the presence of HPV-specific T cells in HPV-associated neoplastic tissue, which suggests the possible existence of systemic immunosuppression and an immunosuppressive TME that significantly influence the efficacy of therapeutic vaccines. Tumor-associated macrophages (TAMs) may cause the disruptive change of antitumor immunity in TME and promote tumor growth and metastasis (61). TAMs are a heterogeneous population of cells that display a range of phenotypes depending on the type of tumor and their location in TME (62,63). TAMs are commonly the most abundant infiltrating leukocytes in most tumors and are predominantly thought to have pro-tumor effects. These include both immunosuppressive effects in The translational expression level (mRNA and protein) of the eight EMT-related genes was positively correlated with disease status as they were upregulated in cervical squamous cell carcinoma and adenocarcinoma samples. ANOVA followed by Tukey's multiple-comparison tests; ** P<0.01 and **** P<0.0001. EMT, epithelial-mesenchymal transition; CDH1, E-cadherin; KRT12, cytokeratin 12; TJP1, tight junction protein 1; CDH2, N-cadherin; VIM, vimentin; ZEB1, zinc finger E-box-binding homeobox 1; ILK, integrin-linked protein kinase; FPKM, fragments per kilobase per million; ERBB3, Erb-B2 receptor tyrosine kinase 3; TPM, transcripts per million. addition to pro-antigenic and metastatic effects. TAMs also promote tumor immune evasion through expression of signal regulatory protein α (SIRPα) (64,65). SIRPα is a receptor for CD47 (66), a cell surface protein that typically protects normal cells from phagocytosis by macrophages or dendritic cells. CD47 is frequently overexpressed on tumor cells and plays a key role in tumor escape by binding to SIRPα and sending macrophages a 'don't eat me' signal (67,68). Blockade of the CD47-SIRPα signal has been shown to stimulate phagocytosis, leading to tumor cell elimination (69).
In the present study, TAMs were revealed to drive tumor angiogenesis and progression in a spontaneous model of cervical cancer through the production of MMP-9. Previous studies (60,70) showed that MMP9 alone did not significantly affect the survival rate of cervical cancer; however, in the present study, when TME macrophages decreased, there was an impact on OS, and the OS with high ERBB3 was significantly reduced. The GO and KEGG enrichment of MMP9 demonstrated that it participated in the biological process of the IL-17 signaling pathway. In patients with bullous pemphigoid, monocyte-derived macrophages, but not lymphocytes, respond to CXCL10 (upregulated by IL-17) by increasing MMP-9 release, potentially creating an inflammatory loop associated with disease outcome (71).
Late after pattern recognition receptors stimulation, bone-marrow-derived dendritic cells (BMDCs) induce the glucose transporter GLUT1 (72,73) and commit to aerobic glycolysis via the mTOR-HIF-1α/iNOS axis (74), which generates NO, inhibiting the electron transport chain. This process might decrease expression of MHC and co-stimulatory molecules by activated BMDCs (75).
Intratumoral immune cells are more often present in tumors with increasing PI3K downstream phosphorylation (76)(77)(78). This is most pronounced for MMP9-positive cells. Research data shows that disturbances in the PI3K pathway may help immune escape (79). A prospective trial in cervical cancer suggested that PI3K pathway alterations may be associated with the composition of TME (80,81).
Previous studies (82)(83)(84) have suggested that when cells undergo EMT and shift progressively from an epithelial to a mesenchymal state, genetic alterations include decreased expression levels of CDH1, cytokeratin 12 and TJP1, increased It can be inferred from the present study that tumor infiltration by CD8-positive lymphocytes was associated with PIK3CA mutations and worse clinical outcome. At the molecular level, EMT transcriptional factors, including SNAl, ZEB1 and TWIST1, regulated immunosuppressive cells or enhanced the expression of immunosuppressive checkpoint molecules through the production of chemokines, thereby resulting in immunosuppression of TME. Immunosuppressive factors can also induce EMT in tumor cells. This mutual feedback between EMT and immunosuppression promotes tumor progression (85).
In conclusion, integrating the characteristics of biomarkers in multiple dimensions can ensure the most efficient management choice for each patient with cancer (Fig. 15). The EMT status cannot be assessed based on one or a few molecular markers, but should be assessed in conjunction with changes in cell characteristics to assess the current ability of cell metastasis and distant invasion. Immuno-oncology research can generate the discriminating power and richness of data   The significance of correlations between immune cell infiltration in cervical cancer tissue and the expression of different genes was determined using Spearman's rank correlation coefficient. a     required for these features. In particular, simple MMP changes are not a prognostic factor in cervical cancer. When ERBB3 activates the PI3K pathway to change immune cell infiltration, the cervical cancer prognosis model is meaningful.

Acknowledgements
Not applicable.

Availability of data and materials
The datasets generated and analyzed during the current study are available in the cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) data in TCGA (https://portal.gdc.cancer.gov/). The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions
XY performed the data analysis work and aided in writing the manuscript. WZ designed the study and edited the manuscript. XY and WZ confirm the authenticity of all the raw data. All authors have read and approved the final manuscript.

Ethics approval and consent to participate
Not applicable.

Patient consent for publication
Not applicable.