Title
Author
DOI
Article Type
Special Issue
Volume
Issue
Functional, microarray and bioinformatics analyses of human papillomavirus E7-specific pathways in cervical carcinogenesis
1Doctoral Program in Obstetrics and Gynecology, Graduate School of Comprehensive Human Sciences, University of Tsukuba, 305-8577 Tsukuba, Japan
2Department of Obstetrics and Gynecology, Institute of Medicine, University of Tsukuba, 305-8575 Tsukuba, Japan
*Corresponding Author(s):minaguchit@md.tsukuba.ac.jp (Takeo Minaguchi)
| History | Submitted: 24 November 2024 | Accepted: 27 December 2024 | Published: 15 May 2025 |
| Copyright: | ©2025 The Author(s). Published by MRE Press. |
Background: Human papillomavirus (HPV) vaccines are still ineffective for already infected patients. The development of new prophylactics and therapeutics for cervical cancer is hence urgently required particularly for contries and regions where medical resources are lacking. The mRNAs encoding the HPV oncogenes E6 and E7 are bicistronic and generated from the same pre-mRNA. Alternative splicing produces different mRNA variants, with the resultant balance of E6/E7 levels impacting downstream functions. Accumulating evidence suggests that E7 may contribute more closely to cervical carcinogenesis than E6. The aim of this study was to explore the E7-specific carcinogenic pathways. Methods: Small interfering RNAs knocking down either E6 or E6/E7 were transfected into CaSki and HeLa cells, and the analyses of cellular effects, microarray and bioinformatics were conducted. The survival of patients with cervical cancer based on gene expressions was analyzed utilizing a web database tool. Results: E7 knockdown induced G1 cell cycle arrest and inhibited cellular proliferation, clonogenicity and transformation. Microarray analyses identified 15 E7-specific differentially expressed genes (DEGs), and their functional annotations included “epithelial-to-mesenchymal transition”, “chromatin remodeling”, “focal adhesion: phosphatidylinositol 3-kinase (PI3K)-AKT serine/threonine kinase (Akt)-mammalian target of rapamycin (mTOR)-signaling pathway” and “DNA-binding transcription factor activity”. Pathway interaction analysis revealed that G1 cell cycle arrest was the most significant and pivotal pathway. Among the E7-specific DEGs, aldo-keto reductase family 1 member B10 (AKR1B10), family with sequence similarity 78 member A (FAM78A) and AHNAK nucleoprotein 2 (AHNAK2) exhibited tumor-suppressive effects, whereas forkhead box A1 (FOXA1), SMAD family member 9 (SMAD9) and plastin 1 (PLS1) showed oncogenic effects in survival analysis, being consistent with the expression fold differences by the microarray. Conclusions: The identified DEGs were suggested to be involved in E7-specific cervical carcinogenesis by cooperating via multiple mechanisms, providing potential novel targets and biomarkers for the more efficient prevention and treatment of cervical cancer.
Cite this article
Nan Qi, Takeo Minaguchi, Kaoru Fujieda, Asami Suto, Hiroya Itagaki, Yuri Tenjimbayashi, Ayumi Shikama, Nobutaka Tasaka, Azusa Akiyama, Sari Nakao, Toyomi Satoh. Functional, microarray and bioinformatics analyses of human papillomavirus E7-specific pathways in cervical carcinogenesis. European Journal of Gynaecological Oncology. 2025; 46(5): 8-23. doi: 10.22514/ejgo.2025.061
1. Introduction
A total of 348,874 women died from cervical cancer and the age-standardized mortality rate was 7.1 in 2022 worldwide (https://gco.iarc.fr/today/home). At present, cervical cancer is the fourth most common cancer and the fourth leading cause of cancer death in women (https://gco.iarc.fr/today/home). The major causative factor for the development of cervical cancer is human papillomavirus (HPV) infection [1]. Although prophylactic HPV vaccines have been shown to reduce cervical cancer risk, they are still ineffective for already infected patients. It is urgently necessary to develop new prophylactics and therapeutics particularly for countries and regions where medical resources are lacking.
HPV has a round virion containing 8-Kb circular double-stranded DNA. The HPV genome contains the six early genes, E1, E2, E4, E5, E6 and E7, as well as the two late genes, L1 and L2, which encode capsid proteins [2]. E6 and E7 encode the oncoproteins, which serve pivotal roles in cervical carcinogenesis by degrading and inactivating the tumor suppressors p53 and Rb, respectively [2]. The major downstream functions of p53 include apoptosis [3], while those of Rb include G1 cell cycle arrest [4]. In addition to these, E6 and E7 mediate a diverse range of other oncogenic pathways [5, 6]. To date, >200 types of HPV have been identified based on the sequence of the L1 gene. High-risk HPVs (HR-HPVs) include types 16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58, 59, 68, 73 and 82 [7], and are responsible for >99.7% of cervical cancers [8]. Among HR-HPVs, HPV16 and 18 are the most prevalent types; 55% of cervical cancers are HPV16-positive and 15% are HPV18-positive [8].
The mRNAs encoding E6 and E7 are bicistronic and are generated from the same pre-mRNA [9]. E6 is mainly translated from unspliced mRNA, and E7 is translated from spliced mRNAs [9]. Alternative splicing produces different mRNA variants, and the resultant balance of E6/E7 protein levels impacts downstream cellular effects [9]. Accumulating evidence [10, 11, 12, 13], including our previous study [14], suggests that E7 may contribute more closely to cervical oncogenesis than E6. Therefore, targeting the E7-specific carcinogenic pathway is thought to provide more efficient prophylactic and therapeutic potentials for cervical cancer. However, the underlying molecular mechanisms are not yet fully clarified. The present study explored the E7-specific pathways of cervical carcinogenesis by utilizing microarray and bioinformatics analyses in HeLa and CaSki human cervical cancer cells, in which E6 or E6/E7 was knocked down. The present findings have significant implications for understanding detailed molecular mechanisms of cervical carcinogenesis, as well as for developing more efficient novel prophylactics and therapeutics for cervical cancer.
2. Materials and methods
2.1 Cell lines and culture
CaSki (HPV16-positive) and HeLa (HPV18-positive) cell lines were purchased from American Type Culture Collection. Cells were maintained in Dulbecco’s Modified Eagle Medium (Nissui Pharmaceutical Co., Ltd., Tokyo, Japan, 05915) containing 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc., Grand Island, NY, USA, 10437-028) and 1 U/mL penicillin-streptomycin-amphotericin B suspension (FUJIFILM Wako Pure Chemical Corporation, Osaka, Japan, 161-23181) at 37 °C with 5% CO2.
2.2 Small interfering RNA (siRNA)
A total of 3 × 105 cells were seeded in 60-mm dishes in antibiotics-free medium and cultured overnight. Cells were transfected with siRNAs at a concentration of 10 nM (HeLa cells) or 40 nM (CaSki cells) using Lipofectamine RNAiMAX (Thermo Fisher Scientific, Inc., Tokyo, Japan, 13778-150) following the manufacturer’s instructions. The siRNAs (16E6E7, 16E6, 18E6 and 18E6E7) were synthesized by GE Healthcare Dharmacon, Inc. (Tokyo, Japan), and their targets and sequences are described in Table 1 (Ref. [15]). The siRNA sequences were subjected to National Center for Biotechnology Information (NCBI)-Nucleotide Basic Local Alignment Search Tool (BLASTN) for similarity against the human Reference Sequence (RefSeq) RNA database with no significant matches identified (https://blast.ncbi.nlm.nih.gov/Blast.cgi). In principle, the full-length mRNA encodes the E6 protein, and the spliced short mRNAs encode the E7 protein [9]. siRNAs 16E6 and 18E6 are designed to target the spliced-out region of mRNA, and siRNAs 16E6E7 and 18E6E7 are designed at the non-spliced-out region of mRNA. Accordingly, 16E6 and 18E6 inhibit only the full-length mRNA, thereby knocking down only the E6 protein, and 16E6E7 and 18E6E7 inhibit both the full-length and the spliced mRNAs, thereby knocking down both the E6 and E7 proteins [15] (Table 1). Non-specific (NS) siRNA was purchased from GE Healthcare Dharmacon, Inc (Tokyo, Japan, D-001210-03-50; Target sequence: AUGUAUUGGCCUGUAUUAG). Cells were analyzed after 24 or 48 h of incubation.
| siRNA | Target | Namea | Sequencea |
| 16E6E7 | HPV16 E6 & E7 | 198 | 5′-GCACACACGUAGACAUUCGdTdT-3′ |
| 16E6 | HPV16 E6 | 209 | 5′-UCCAUAUGCUGUAUGUGAUdTdT-3′ |
| 18E6 | HPV18 E6 | 219 | 5′-CUCUGUGUAUGGAGACACAdTdT-3′ |
| 18E6E7 | HPV18 E6 & E7 | 220 | 5′-UGGAGUUAAUCAUCAACAUdTdT-3′ |
| aThe names and sequences of the siRNAs are based on a previous study [15]. HPV: human papillomavirus; siRNA: small interfering RNA. |
2.3 Western blotting
Proteins were extracted from cells using M-PER mammalian protein extraction reagent (Thermo Fisher Scientific, Inc., Tokyo, Japan, 78503) supplemented with 1X cOmplete Mini protease inhibitor cocktail (Roche Diagnostics, Tokyo, Japan, 04693124001) and 1X PhosSTOP phosphatase inhibitor cocktail (Roche Diagnostics, Tokyo, Japan, 04906845001), separated by 10% sodium dodecyl sulfate (SDS)-polyacrylamide gel electrophoresis (PAGE) and transferred to nitrocellulose membranes. Membranes were blocked with 5% non-fat dry milk and incubated at 4 °C overnight with primary antibody. The following primary antibodies were used: HPV type 16 E6 polyclonal antibody (1:1000; PA5-117355; Thermo Fisher Scientific, Inc., Tokyo, Japan), HPV16/18 E6 (1:200; sc-460; Santa Cruz Biotechnology, Inc., Santa Cruz, CA, USA), HPV16 E7 (1:200; sc-51951; Santa Cruz Biotechnology, Inc., Santa Cruz, CA, USA), HPV18 E7 (1:100; sc-365035; Santa Cruz Biotechnology, Inc., Santa Cruz, CA, USA), p53 (1:1000; sc-126; Santa Cruz Biotechnology, Inc., Santa Cruz, CA, USA), anti-human Retinoblastoma protein (1:1000; 554136; Becton, Dickinson and Company, Franklin Lakes, NJ, USA) and anti-β-tubulin (1:1000; 556321; Becton, Dickinson and Company, Franklin Lakes, NJ, USA). After incubation with goat anti-mouse/rabbit IgG (H + L) HRP-conjugated secondary antibody (1:5000; W4021/W4011; Promega Corporation, Madison, WI, USA), proteins were detected by enhanced chemiluminescence using ECL Select Western Blotting Detection Reagent (GE HealthCare Japan, Tokyo, Japan, RPN2235). The intensity of the target signal was detected using a LAS 500 image analyzer (Cytiva, Tokyo, Japan) and quantified using ImageJ 1.53k software (https://imagej.net/ij/).
2.4 Cell proliferation assay
A total of 2.6 × 105 cells were seeded in 100-mm dishes, incubated overnight and transfected with the siRNAs. After 24 h of incubation, cells were trypsinized, collected and counted using an auto cell counter TC10 (Bio-Rad Laboratories, Inc., Hercules, CA, USA).
2.5 Colony formation assay
A total of 2000 cells were seeded in 6-well plates in triplicate and incubated overnight. Cells were transfected with the siRNAs and cultured in an incubator for 3 weeks. Cells were fixed with 100% methanol and stained with Giemsa, and the colonies consisting of >50 cells in each well were counted manually under a microscope.
2.6 Flow cytometry
A total of 7.8 × 105 cells were seeded in 100-mm dishes, incubated overnight and transfected with the siRNAs. After 24 h of incubation, cells were trypsinized, collected, fixed in ice-cold 70% ethanol and stained with ice-cold propidium iodide (PI)/phosphate-buffered saline (PBS)-Triton/ribonuclease (RNase) A. The cell cycle distribution was analyzed using a FACSVerse flow cytometer (Becton, Dickinson and Company, Franklin Lakes, NJ, USA).
2.7 TdT-mediated dUTP nick-end labeling (TUNEL) assay
A total of 1.2 × 105 cells were seeded on glass coverslips in 6-well plates, incubated overnight and transfected with the siRNAs. After 48 h of incubation, apoptotic cells were analyzed under a fluorescence microscope (Zeiss Axio Imager A1, Zeiss, Oberkochen, BW, Germany) using a DeadEnd Fluorometric TUNEL System (Promega Corporation) according to the manufacturer’s instructions.
2.8 Wound healing assay
A total of 2.4 × 105 cells were seeded in 6-well plates and incubated until 90–100% confluent. The cell monolayer was scraped in a straight line with a sterile p1000 pipet tip, and the first images of the scratch at three random points were acquired under a phase-contrast microscope. The cells were transfected with the siRNAs and incubated for 24 h in serum-free medium and the second images were acquired. Distances (μm) between one side of the scratch and the other were measured using ImageJ 1.53k software.
2.9 Transformation assay
Cell transformation was analyzed using a CytoSelect 96-well Cell Transformation Assay Kit (CBA-135; Cell Biolabs, Inc., San Diego, CA, USA) according to the manufacturer’s instructions. A total of 3 × 105 cells were seeded in 60-mm dishes, incubated overnight and transfected with the siRNAs. After 24 h of incubation, cells were trypsinized and collected, and cell suspension/agar matrix was prepared. Cell suspension/agar matrix including 2.5 × 103 transfected cells was added into each well of a 96-well sterile microplate containing 50 μL each of agar matrix layer. After 7 days of incubation, anchorage-independent growth was measured based on the absorbance at 570 nm using a plate reader (Sunrise; Tecan Group, Ltd., Männedorf, ZH, Switzerland).
2.10 Microarray and bioinformatics analyses
Total RNAs were extracted using TRIzol reagent (Thermo Fisher Scientific, Inc., Tokyo, Japan, 15596026) from CaSki and HeLa cells transfected with the siRNAs and submitted to Macrogen, Inc., and gene expression profiling was conducted using a SurePrint G3 Human Gene Expression 8 × 60k v3 Microarray (Agilent Technologies, Inc.). Differentially expressed genes (DEGs) were identified using Biometrics Research Branch (BRB)-Array tools v4.6.2 Beta 1 (https://brb.nci.nih.gov/BRB-ArrayTools/) using the class comparison analysis with pairing samples by cell line and setting p < 0.05 and |log fold change| ≥ 2 as the cutoff criteria. E7-specific DEGs were extracted using a Venn diagram illustrated by FunRich 3.1.3 software (http://funrich.org/index.html). Gene Ontology (GO) and pathways for the DEGs were identified by functional annotation analyses using the Database for Annotation, Visualization and Integrated Discovery v2022 (DAVID; https://david.ncifcrf.gov/home.jsp; p < 0.05; gene count ≥2; Expression Analysis Systematic Explorer (EASE) score ≤0.1) [16] and Metascape v3.5 (https://metascape.org; p < 0.01; minimum overlap, 3; minimum enrichment, 1.5) [17]. Protein-Protein Interactions (PPIs) were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins v11.5 (STRING; https://string-db.org/; confidence score >0.4) [18] and hub genes, modules and pathway interactions were analyzed using Cytoscape 3.9.1 software (https://cytoscape.org/) [19] using cytoHubba (maximal clique centrality method), Molecular Complex Detection (MCODE) (degree cutoff, 2; node score cutoff, 0.2; k-core, 2; max. depth, 100) and EnrichmentMap (connectivity cutoff of Jaccard similarity, 0.4) tools, respectively. Recurrence-free survival (RFS) in patients with cervical squamous cell carcinoma (n = 174) was compared between groups of patients based on the mRNA expression levels of the DEGs by Kaplan-Meier plotter (https://kmplot.com/analysis/) which utilized gene expression data and survival information from Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/), European Genome-phenome Archive (EGA; https://ega-archive.org/), and The Cancer Genome Atlas (TCGA; https://www.cancer.gov/ccg/research/genome-sequencing/tcga), and p-values were calculated by the log-rank test [20].
2.11 Statistical analysis
All data are presented as the mean ± standard deviation (SD) of at least three independent experiments. Differences were compared by one-way analysis of variance (ANOVA) with Tukey’s post-hoc test or unpaired Student’s t-test using R version 4.0.5 (https://www.r-project.org/). p < 0.05 was considered to indicate a statistically significant difference.
3. Results
3.1 E6/E7 knockdown in CaSki and HeLa cells
In HPV16-positive CaSki cells, both siRNAs 16E6E7 and 16E6 decreased E6 expression and increased p53 expression compared with that in the mock and siRNA NS groups, and siRNA 16E6E7 decreased E7 expression and increased Rb expression compared with that in the mock, siRNA NS and 16E6 groups (Fig. 1A; Table 1). In HPV18-positive HeLa cells, both 18E6 and 18E6E7 decreased E6 expression and increased p53 expression compared with that in the mock and NS groups, and 18E6E7 decreased E7 expression and increased Rb expression compared with that in the mock, NS and 18E6 groups (Fig. 1A; Table 1). Accordingly, comparing the cellular effects of 16E6 with those of mock and NS should indicate mainly the results of HPV16 E6 inhibition, and comparing the cellular effects of 16E6E7 with those of 16E6 should indicate mainly the results of HPV16 E7 inhibition (Table 1). Comparing the cellular effects of 18E6 with those of mock and NS should indicate mainly the results of HPV18 E6 inhibition, and comparing the cellular effects of 18E6E7 with those of 18E6 should indicate mainly the results of HPV18 E7 inhibition (Table 1).
3.2 Effect of knockdown of E6/E7 on cellular proliferation and clonogenicity
Cellular proliferation was significantly inhibited in the 16E6E7 group compared with the 16E6 group in CaSki cells (p < 0.05; Fig. 1B,C; Table 2). Colony formation was significantly inhibited in the 16E6E7 group compared with the 16E6 group, and in the 16E6 and 18E6 groups compared with the mock group in CaSki and HeLa cells, respectively (p < 0.05, p < 0.01 and p < 0.05, respectively; Fig. 1D–F; Table 2).
| Cellular functions | CaSki | HeLa | ||||
| 16E6 vs. Mock | 16E6 vs. NS | 16E6E7 vs. 16E6 | 18E6 vs. Mock | 18E6 vs. NS | 18E6E7 vs. 18E6 | |
| Cell number | 54% | 120% | 31%a | 35% | 50% | 58% |
| Colony number | 44%b | 102% | 12%a | 20%a | 29% | 58% |
| G0/G1 population | 122%b | 103% | 120%b | 132%b | 118%a | 131%a |
| Sub-G1 population | 72% | 67% | 90% | 152% | 193% | 87% |
| TUNEL-positive cells | 298%a | 134% | 186% | 419%b | 235%a | 111% |
| Cell migration | 89% | 102% | 85% | 57%a | 73% | 84% |
| Cell transformation | 64% | 76% | 89% | 57%a | 61%a | 62%b |
| ap < 0.05 and bp < 0.01, one-way ANOVA with Tukey’s post-hoc test or unpaired Student’s t-test. NS: Non-specific; TUNEL: TdT-mediated dUTP nick-end labeling. |
Fig. 1.Knockdown of E6/E7 by transfection with mock, NS, 16E6E7, 16E6, 18E6 or 18E6E7 siRNAs in CaSki and HeLa cells (Table 1). (A) Western blotting of E6, E7, p53, Rb and Tubulin, and the relative intensity to Tubulin. (B) Cell numbers following knockdown of E6 based on the cell proliferation assay. (C) Cell numbers following knockdown of E7 based on the cell proliferation assay. (D) Representative well images of the colony formation assay. (E) Colony numbers following knockdown of E6 based on the colony formation assay. (F) Colony numbers following knockdown of E7 based on the colony formation assay. Data are presented as the mean ± SD. Differences were compared by one-way ANOVA with Tukey’s post-hoc test or unpaired Student’s t-test. *p < 0.05, **p < 0.01. NS: Non-specific; siRNA: small interfering RNA.
3.3 Effect of knockdown of E6/E7 on the cell cycle and apoptosis
Flow cytometry revealed that the G0/G1 cell cycle population was significantly increased in the 16E6E7 group compared with the 16E6 group, and in the 16E6 group compared with the mock group in CaSki cells (p < 0.01 and p < 0.01, respectively; Fig. 2A–C; Table 2), as well as in the 18E6 group compared with the mock and NS groups, and in the 18E6E7 group compared with the 18E6 group in HeLa cells (p < 0.01, p < 0.05 and p < 0.05, respectively; Fig. 2A–C; Table 2). The sub-G1 population was not significantly changed in either CaSki or HeLa cells (Fig. 2A,D,E; Table 2). The TUNEL assay showed that 16E6 significantly increased the proportion of apoptotic cells compared with that in the mock group in CaSki cells (p < 0.05; Fig. 2F,G; Table 2), and that 18E6 significantly increased the proportion of apoptotic cells compared with that in the mock and NS groups in HeLa cells (p < 0.01 and p < 0.05, respectively; Fig. 2F–H; Table 2).
Fig. 2.Analyses of the cell cycle and apoptosis in CaSki and HeLa cells transfected with mock, NS, 16E6E7, 16E6, 18E6 or 18E6E7 small interfering RNAs (Table 1). (A) Cell cycle distribution. (B) Percentages of cells in the G0/G1 phase following knockdown of E6. (C) Percentages of cells in the G0/G1 phase following knockdown of E7. (D) Percentages of cells in the sub-G1 phase following knockdown of E6. (E) Percentages of cells in the sub-G1 phase following knockdown of E7. (F) Representative images of the TUNEL assay (magnification, ×200). (G) Percentages of TUNEL-positive apoptotic cells following knockdown of E6. (H) Percentages of TUNEL-positive apoptotic cells following knockdown of E7. Data are presented as the mean ± SD. Differences were compared by one-way ANOVA with Tukey’s post-hoc test or unpaired Student’s t-test. *p < 0.05, **p < 0.01. NS: Non-specific; TUNEL: TdT-mediated dUTP nick-end labeling; PI-A: propidium iodide-area; DAPI: 4′,6-diamidino-2-phenylindole.
3.4 Effect of knockdown of E6/E7 on cell migration and transformation
The wound healing assay revealed that 18E6 significantly inhibited wound healing compared with that in the mock group in HeLa cells (p < 0.05; Fig. 3A–C; Table 2). The transformation assay showed that 18E6 significantly inhibited anchorage-independent growth compared with the mock and NS groups in HeLa cells (p < 0.05 and p < 0.05, respectively; Fig. 3D; Table 2), and that 18E6E7 also significantly inhibited anchorage-independent growth compared with the 18E6 group in HeLa cells (p < 0.01; Fig. 3E; Table 2).
Fig. 3.Analyses of migration and transformation in CaSki and HeLa cells transfected with mock, NS, 16E6E7, 16E6, 18E6 or 18E6E7 small interfering RNAs (Table 1). (A) Representative images of the wound healing assay (magnification, ×10). (B) Percentages of cell migration following knockdown of E6 based on the wound healing assay. (C) Percentages of cell migration following knockdown of E7 based on the wound healing assay. (D) Anchorage-independent growth (absorbance at 570 nm) following knockdown of E6 based on the transformation assay. (E) Anchorage-independent growth following knockdown of E7 based on the transformation assay. Data are presented as the mean ± SD. Differences were compared by one-way ANOVA with Tukey’s post-hoc test or unpaired Student’s t-test. *p < 0.05, **p < 0.01. NS: Non-specific.
3.5 Microarray and bioinformatics analyses
Microarray analyses were subsequently conducted using RNAs extracted from CaSki and HeLa cells in which E6/E7 were knocked down using the same siRNAs as in the aforementioned experiments. A total of 13,980 genes were analyzed after normalization and filtration. A total of 31 genes were identified as E7-associated DEGs by comparing between 16E6- and 16E6E7-transfected CaSki cells, and between 18E6- and 18E6E7-transfected HeLa cells (Fig. 4A,B). A total of 150 genes were identified as E6-associated DEGs by comparing between NS- and 16E6-transfected CaSki cells, and between NS- and 18E6-transfected HeLa cells (Fig. 4A,B). A total of 222 genes were identified as E6, E7-associated DEGs by comparing between NS- and 16E6E7-transfected CaSki cells, and between NS- and 18E6E7-transfected HeLa cells (Fig. 4A,B). Using FunRich software, the present study further determined the DEGs identified in both the E7-associated group and the E6, E7-associated group (Fig. 4B). One E6-associated gene was excluded, and finally 15 genes were extracted as the E7-specific DEGs (Fig. 4B; Table 3). Functional annotation analyses of the 15 DEGs using Metascape identified “positive regulation of DNA-binding transcription factor activity”, “focal adhesion: phosphatidylinositol 3-kinase (PI3K)-AKT serine/threonine kinase (Akt)-mammalian target of rapamycin (mTOR)-signaling pathway” and “chromatin remodeling” (Fig. 4C,D; Table 4). Functional annotation analyses of the 15 DEGs using DAVID identified “negative regulation of epithelial to mesenchymal transition”, “HPV infection”, “PI3K-Akt signaling pathway” and “focal adhesion: PI3K-Akt-mTOR-signaling pathway” (Table 5).
| Gene symbol | Description | Fold difference of geometric mean gene expressions | Parametric p-valuea | Log-ratio of normalized gene expressions | |
| CaSki | HeLa | ||||
| SLC44A1 | solute carrier family 44 member 1 | 3.07 | <0.001 | 1.52 | 1.71 |
| AKR1B10 | aldo-keto reductase family 1 member B10 | 0.43 | <0.001 | −1.11 | −1.34 |
| RTKN2 | rhotekin 2 | 2.12 | <0.001 | 1.07 | 1.09 |
| NFE2 | nuclear factor, erythroid 2 | 0.50 | 0.001 | −1.14 | −0.88 |
| ARL5B | ADP ribosylation factor like GTPase 5B | 2.11 | 0.003 | 0.75 | 1.40 |
| FOXA1 | forkhead box A1 | 2.20 | 0.005 | 0.68 | 1.60 |
| PPP2R2B | protein phosphatase 2 regulatory subunit Bbeta | 2.54 | 0.010 | 0.62 | 2.07 |
| SMAD9 | SMAD family member 9 | 2.49 | 0.011 | 0.58 | 2.05 |
| FAM78A | family with sequence similarity 78 member A | 0.48 | 0.011 | −0.51 | −1.60 |
| PTEN | phosphatase and tensin homolog | 2.37 | 0.012 | 0.55 | 1.95 |
| GALNT1 | polypeptide N-acetylgalactosaminyltransferase 1 | 2.53 | 0.012 | 0.57 | 2.11 |
| H1-1 | H1.1 linker histone, cluster member | 2.14 | 0.013 | 0.41 | 1.78 |
| PLS1 | plastin 1 | 2.33 | 0.026 | 0.35 | 2.09 |
| CCNE1 | cyclin E1 | 2.06 | 0.027 | 0.32 | 1.78 |
| AHNAK2 | AHNAK nucleoprotein 2 | 0.49 | 0.030 | −0.28 | −1.75 |
| aPaired t-test with random variance model. |
Fig. 4.Microarray analysis and functional annotation using RNAs from small interfering RNA-transfected cells. (A) Volcano plots of E7-associated DEGs (comparing 16E6- vs. 16E6E7-transfected CaSki cells and 18E6- vs. 18E6E7-transfected HeLa cells), E6-associated DEGs (comparing NS- vs. 16E6-transfected CaSki cells and NS- vs. 18E6-transfected HeLa cells) and E6, E7-associated DEGs (comparing NS- vs. 16E6E7-transfected CaSki cells and NS- vs. 18E6E7-transfected HeLa cells). p-values (-log10) are plotted against FC (log2). Upregulated (log2FC ≥1 and p < 0.05) and downregulated (log2FC ≤−1 and p < 0.05) DEGs are denoted by blue nodes. (B) Venn diagram showing the overlap among E7-associated, E6-associated and E6, E7-associated DEGs. (C) Functional annotation using Metascape of the 15 E7-specific DEGs extracted from the Venn diagram. (D) Enriched ontology clusters based on Metascape analysis of the 15 E7-specific DEGs. FC: fold change. DEGs: differentially expressed genes; NS: Non-specific; GO: Gene Ontology; PI3K: phosphatidylinositol 3-kinase; Akt: AKT serine/threonine kinase; mTOR: mammalian target of rapamycin.
| Group ID | Category | Term | Description | Log10p-value | Symbols |
| 1_Summary | GO Biological Processes | GO:0051091 | Positive regulation of DNA-binding transcription factor activity | −3.55 | FOXA1, PTEN, RTKN2 |
| 1_Member | GO Biological Processes | GO:0051091 | Positive regulation of DNA-binding transcription factor activity | −3.55 | FOXA1, PTEN, RTKN2 |
| 1_Member | GO Biological Processes | GO:0051090 | Regulation of DNA-binding transcription factor activity | −2.87 | FOXA1, PTEN, RTKN2 |
| 2_Summary | WikiPathways | WP3932 | Focal adhesion: PI3K-Akt-mTOR-signaling pathway | −3.38 | FOXA1, PPP2R2B, PTEN, CCNE1 |
| 2_Member | WikiPathways | WP3932 | Focal adhesion: PI3K-Akt-mTOR-signaling pathway | −3.38 | FOXA1, PPP2R2B, PTEN |
| 2_Member | KEGG Pathway | hsa05165 | Human papillomavirus infection | −3.27 | CCNE1, PPP2R2B, PTEN |
| 2_Member | WikiPathways | WP4172 | PI3K-Akt signaling pathway | −3.24 | CCNE1, PPP2R2B, PTEN |
| 2_Member | KEGG Pathway | hsa04151 | PI3K-Akt signaling pathway | −3.18 | CCNE1, PPP2R2B, PTEN |
| 3_Summary | GO Biological Processes | GO:0006338 | Chromatin remodeling | −3.20 | H1-1, FOXA1, NFE2 |
| 3_Member | GO Biological Processes | GO:0006338 | Chromatin remodeling | −3.20 | H1-1, FOXA1, NFE2 |
| 3_Member | GO Biological Processes | GO:0006325 | Chromatin organization | −2.51 | H1-1, FOXA1, NFE2 |
| GO: Gene Ontology; FOXA1: forkhead box A1; PTEN: phosphatase and tensin homolog; RTKN2: rhotekin 2; PPP2R2B: protein phosphatase 2 regulatory subunit Bbeta; CCNE1: cyclin E1; H1-1: H1.1 linker histone, cluster member; NFE2: nuclear factor, erythroid 2; PI3K: phosphatidylinositol 3-kinase; Akt: AKT serine/threonine kinase; mTOR: mammalian target of rapamycin; ID: identifier. |
| Term | Description | Count | % | p-value | Genes | List total | Pop hits | Pop total | Fold enrichment | Bonferroni | Benjamini | FDR |
| GO:0010719 | Negative regulation of epithelial to mesenchymal transition | 2 | 13 | 0.025 | FOXA1, PTEN | 14 | 38 | 19308 | 72.6 | 1 | 1 | 1 |
| hsa05165 | Human papillomavirus infection | 3 | 20 | 0.030 | PPP2R2B, CCNE1, PTEN | 8 | 331 | 8164 | 9.25 | 0.81 | 0.88 | 0.88 |
| hsa04151 | PI3K-Akt signaling pathway | 3 | 20 | 0.034 | PPP2R2B, CCNE1, PTEN | 8 | 354 | 8164 | 8.65 | 0.85 | 0.88 | 0.88 |
| WP3932 | Focal adhesion: PI3K-Akt-mTOR-signaling pathway | 3 | 20 | 0.035 | FOXA1, PPP2R2B, PTEN | 9 | 302 | 7879 | 8.70 | 0.92 | 1 | 1 |
| WP4172 | PI3K-Akt signaling pathway | 3 | 20 | 0.043 | PPP2R2B, CCNE1, PTEN | 9 | 339 | 7879 | 7.75 | 0.96 | 1 | 1 |
| FDR, false discovery rate; GO: Gene Ontology; FOXA1: forkhead box A1; PTEN: phosphatase and tensin homolog; RTKN2: rhotekin 2; PPP2R2B: protein phosphatase 2 regulatory subunit Bbeta; CCNE1: cyclin E1; PI3K: phosphatidylinositol 3-kinase; Akt: AKT serine/threonine kinase; mTOR: mammalian target of rapamycin. |
The subsequent functional enrichment analysis of the 15 E7-specific DEGs using STRING identified the PPI network shown in Fig. 5A. The Cytoscape software identified one module (Fig. 5B) and the top 8 hub genes (Fig. 5C; Table 6) in the PPI network. The identified pathway interactions for the PPI network based on Cytoscape are shown in Fig. 5D, and the pathways included “G1 cell cycle arrest”, “focal adhesion: PI3K-Akt-mTOR-signaling”, “formation of senescence-associated heterochromatin foci (SAHF)” and “apoptotic signaling pathway”.
| Rank | Name | Gene | Score |
| 1 | 9606.ENSP00000269305 | TP53 | 56 |
| 2 | 9606.ENSP00000262643 | CCNE1 | 54 |
| 2 | 9606.ENSP00000266970 | CDK2 | 54 |
| 4 | 9606.ENSP00000361021 | PTEN | 51 |
| 5 | 9606.ENSP00000281708 | FBXW7 | 24 |
| 5 | 9606.ENSP00000228872 | CDKN1B | 24 |
| 7 | 9606.ENSP00000244573 | H1-1 | 6 |
| 8 | 9606.ENSP00000250448 | FOXA1 | 2 |
| FOXA1: forkhead box A1; PTEN: phosphatase and tensin homolog; CCNE1: cyclin E1; H1-1: H1.1 linker histone, cluster member; TP53: tumor protein p53; CDK2: cyclin dependent kinase 2; FBXW7: F-box and WD repeat domain containing 7; CDKN1B: cyclin dependent kinase inhibitor 1B. |
Fig. 5.Analyses of the PPI network, module, hub genes and pathway interactions of the 15 E7-specific DEGs by STRING and Cytoscape. (A) PPI network generated using STRING. (B) Module of the PPI network. (C) Top 8 hub genes of the PPI network based on the maximal clique centrality method. (D) Pathway interactions of the PPI network. PPI: Protein-Protein Interaction; DEGs: differentially expressed genes; STRING: Search Tool for the Retrieval of Interacting Genes/Proteins.
The prognosis of patients with cervical squamous cell carcinoma was further analyzed according to the mRNA expression levels of the E7-specific DEGs using the web-based survival analysis tool Kaplan-Meier plotter. High expression levels of aldo-keto reductase family 1 member B10 (AKR1B10), family with sequence similarity 78 member A (FAM78A) and AHNAK nucleoprotein 2 (AHNAK2) (p = 0.00095, p = 0.0034 and p = 0.014; Fig. 6) were significantly associated with improved RFS, whereas high expression levels of forkhead box A1 (FOXA1), SMAD family member 9 (SMAD9) and plastin 1 (PLS1) were significantly associated with worse RFS (p = 0.0014, p = 0.027 and p = 0.039; Fig. 6).
Fig. 6.Recurrence-free survival curves obtained utilizing Kaplan-Meier plotter in patients with cervical squamous cell carcinoma (n = 174) according to the mRNA expression levels of the E7-specific DEGs. HR: hazard ratio; AKR1B10: aldo-keto reductase family 1 member B10; FOXA1: forkhead box A1; SMAD9: SMAD family member 9; FAM78A: family with sequence similarity 78 member A; PLS1: plastin 1; AHNAK2: AHNAK nucleoprotein 2.
4. Discussion
Our previous study on separately analyzing E6/E7 mRNAs in liquid-based cytology samples indicated that the presence of E7 mRNAs correlated with progression from low-grade cervical intraepithelial neoplasia (CIN) to invasive cancer by contrast with the presence of E6 mRNA, and that the presence of both E6 and E7 mRNAs was associated with upgraded abnormal cytology in the followed-up patients with CIN1–2 in contrast to the presence of E6 mRNA [14]. These findings are in line with multiple reports, where the positive rate of E7 mRNA/protein expression increased with CIN grade [10, 11, 12], and E7 mRNA predicted progression to high-grade squamous intraepithelial lesion/CIN2+ [13]. Furthermore, in vitro studies have revealed that the expression of E7 by itself could immortalize human keratinocytes at a low frequency but E6 had no such activity, and that the combined expression of E6 and E7 could efficiently immortalize most types of primary cells [21, 22]. Furthermore, in a transgenic mouse model, E7 alone, but not E6 alone, was sufficient to induce high-grade CIN and invasive cancer, and E6 addition resulted in larger and more extensive tumors [23]. Taken together, accumulating evidence suggests that E7 may contribute more closely to cervical carcinogenesis than E6. Although numerous molecules and pathways downstream of E6 and E7 have been identified, the precise mechanism whereby E7 exhibits stronger involvement in cervical carcinogenesis than E6 is yet to be fully clarified, and this was explored in the present study.
The cellular effects of knocking down E6/E7 in HPV-positive CaSki and HeLa cells were evaluated, as summarized in Table 2. The cell proliferation and colony formation assays suggested that knockdown of E7 and knockdown of E6 both inhibited cellular proliferation and/or clonogenicity (Fig. 1B,D; Table 2). Flow cytometry suggested that knockdown of E7 and knockdown of E6 both induced G1 cell cycle arrest (Fig. 2A,B; Table 2). The TUNEL assay suggested that knockdown of E6 induced apoptosis (Fig. 2E; Table 2). The wound healing assay suggested that knockdown of E6 inhibited cell migration (Fig. 3B; Table 2). The transformation assay suggested that knockdown of E7 and knockdown of E6 both inhibited cellular transformation (Fig. 3C; Table 2). These observed differential cellular effects of knockdown of E6 or E7 were consistent with published findings [24].
Flow cytometry exhibited that knockdown of E7 induced G1 cell cycle arrest (Fig. 2B; Table 2), and the pathway interaction for the E7-specific PPI network using Cytoscape revealed that G1 cell cycle arrest was the most significant and pivotal pathway (Fig. 5D). These findings suggest that G1 cell cycle arrest may be the most important E7-specific carcinogenic pathway. Besides G1 cell cycle arrest, the functional annotations based on Metascape and DAVID analyses both included “focal adhesion: PI3K-Akt-mTOR-signaling pathway”, those of Metascape analysis included “chromatin remodeling”, and those of DAVID analysis included “negative regulation of epithelial to mesenchymal transition” (Tables 4 and 5). All of these pathways reportedly contribute to cellular transformation [25, 26, 27, 28]. The transformation assay demonstrated that knockdown of E7 inhibited cellular transformation (Fig. 3E; Table 2), suggesting that these three pathways may contribute to the E7-specifc cervical carcinogenesis via cellular transformation. Additionally, these pathways including G1 cell cycle arrest as well as PI3K-Akt-mTOR-signaling pathway are considered to provide the ideal molecular targets for more effective treatment of cervical cancer.
The survival analyses of groups of patients based on the expression levels of the E7-specific DEGs showed that AKR1B10, FAM78A and AHNAK2 were associated with improved prognosis, while FOXA1, SMAD9 and PLS1 were associated with worse prognosis in cervical cancer (Fig. 6). These oncogenic or tumor-suppressive prognostic significances were consistent with the fold differences of the gene expression levels in the microarray analyses (Table 3). These findings are further supported by published findings. AKR1B10, a member of the aldo/keto reductase superfamily, catalyzes the conversion of retinal to retinol, and thus, counteracts the formation of retinoic acid, which is involved in cell proliferation and differentiation [29]. High expression levels of FAM78A are associated with improved overall survival in pancreatic adenocarcinoma [30]. Knockdown of AHNAK2, a large nucleoprotein, has been reported to lead to radioresistance in esophageal squamous cell carcinoma cells [31]. Knockdown of FOXA1, the founding member of the FOX family of transcription factors [32], reverses chemoresistance by suppressing cell proliferation, migration and epithelial-mesenchymal transition, and high FOXA1 expression is associated with chemoresistant cells and worse patient survival in lung adenocarcinoma [33]. A single nucleotide polymorphism in SMAD9, a member of the SMAD family that transduces signals from the transforming growth factor (TGF)-β pathway, is associated with unfavorable survival in non-small cell lung cancer after radiotherapy [34]. PLS1, a member of the actin-binding protein family, has been implicated to promote metastasis of colorectal cancer [35]. Therefore, it is suggested that these DEGs may serve as useful biomarkers for the prognosis of cervical cancer.
Tumor protein p53 (TP53) is one of the most mutated human gene involved in numerous cellular functions [3], and is known to be involved even in Rb-mediated G1/S cell cycle checkpoint through cyclin dependent kinase inhibitor 1A/p21 [36]. Therefore, the finding that TP53 ranked as the most significant hub gene of the E7-specific PPI network should be reasonable (Table 6). Although Rb was not among the top 8 hub genes of the E7-specific PPI network, cyclin E1 and cyclin dependent kinase 2 (CDK2), ranked as the second and third hub genes (Table 6), are also known to be involved in the Rb-mediated G1/S cell cycle checkpoint [37]. Phosphatase and tensin homolog (PTEN) and cyclin dependent kinase inhibitor 1B/p27, ranked as the fourth and sixth hub genes (Table 6), are also known to be involved in the Rb-mediated G1/S checkpoint through the PI3K/Akt signaling pathway [38, 39]. Additionally, E7 reportedly interacts with numerous proteins other than Rb, and E6 interacts with numerous proteins other than p53 [24], supporting the present results on the hub genes in the PPI network of the E7-specific DEGs.
This study has some limitations. First, the siRNA experiments were based on transient transfection, not stable transfection. Second, the data on E7-specific DEGs were extracted by comparing results indirectly from siRNA experiments. Third, NS siRNA decreased E6 levels in both cell lines, and decreased E7 levels in HeLa cells, compared with those in the mock group. As possible consequences, NS siRNA significantly inhibited colony formation and increased the G0/G1 cell cycle population in CaSki cells compared with those in the mock group. Although other NS siRNAs were tested, the results were similar. Therefore, the significant cellular effects of knockdown of E6 were evaluated by comparing with mock as well as NS. Lastly, reverse transcription-polymerase chain reaction (PCR) experiments using the same siRNAs are necessary to validate the DEGs and hub genes to be E7-specific. Nevertheless, our reasonable findings in view of the previous multiple publications support the significance of the current study. Further exploring the E7-specific pathways in vivo as well as across other types of HR-HPVs will provide beneficial information for the management of cervical neoplastic diseases.
5. Conclusions
The present study explored the precise E7-specific pathways involved in cervical carcinogenesis. The current findings suggest that the identified DEGs contribute to the E7-specific cervical carcinogenesis by cooperating through the multiple pathways, providing significant implications for novel targets and biomarkers for more efficient prophylactic and therapeutic strategies for cervical cancer.
Availability of data and materials
The data generated in this study may be requested from the corresponding author.
Author contributions
NQ—conducted the experiments; acquired data and drafted the manuscript. TM—conceptualized the study; conducted statistical analysis and edited the manuscript. KF, AsS, HI, YT, AyS, NT, AA, SN and TS—interpreted data and critically reviewed the manuscript. AA—acquired funding. TS—supervised the study. NQ and TM—confirm the authenticity of all the raw data. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Not applicable.
Acknowledgment
Not applicable.
Funding
This study was partially supported by a Grant-in-Aid for Scientific Research (grant no. 21K09464) from the Ministry of Education, Culture, Sports, Science and Technology, Tokyo, Japan.
Conflict of interest
The authors declare no conflict of interest.
