This article is mentioned in:

Introduction

Esophageal carcinoma (ESCA) is the seventh most commonly encountered cancer worldwide. According to the latest global cancer statistics, 604,000 estimated new cases of esophageal cancer and 544,000 related deaths have been reported, with the mortality rate ranking sixth among all types of cancer (1). Adenocarcinoma (AC) and squamous cell carcinoma (SCC) are the most frequent histological subtypes of ESCA, of which SCC cases comprise 90% of all ESCA cases; China, Central Asia, East Africa and South Africa have the highest rates of ESCA (2). A large body of evidence suggests that smoking, alcoholism, obesity, poor nutrition and viral infection may increase the risk of developing ESCA (3). At present, the treatment for ESCA is relatively simple and mainly involves esophagectomy; however, this may lead to severe post-operative complications and can markedly increase the risk of mortality (4). Hence, detailed insight into the mechanisms of the development and progression of ESCA are crucial for the identification of novel therapeutic targets.

Human papillomavirus (HPV) is a double-stranded DNA virus that easily infects the mucosal epithelial tissues of the anogenital and upper aerodigestive tracts (5). High-risk HPV infections are the leading cause of cervical carcinomas, and almost all cervical cancers are the result of HPV infection (6). There is increasing evidence to indicate that HPV is also associated with the development of vaginal, penile, anal, and head and neck cancers (5). Globally, 570,000 female and 60,000 male cancer-related cases are attributable to HPV each year, demonstrating an increasing trend (7). Although the association between HPV and ESCA is generally controversial, studies have found that ESCA samples are accompanied by HPV infection, and HPV infection is a negative prognostic factor for patients with ESCA (8,9). It has been previously demonstrated by the authors' laboratory that HPV18 E6E7 induces the immortalization and malignant transformation of human embryonic esophageal epithelial cells (10,11). Therefore, this evidence may suggest that HPV infection is indeed one of the key oncogenic factors for ESCA.

Cell death can be categorized as accidental cell death and regulated cell death (RCD). RCD includes apoptosis, necrosis, autophagy, pyroptosis, ferroptosis, entotic cell death, netotic cell death, parthanatos, lysosome-dependent cell death, alkaliptosis and oxeiptosis (12,13). Several studies have revealed that HPV inhibits cell death to promote the malignancy of cancer cells. For example, HPV E7 inhibits the necrosis of SiHa cells, while E7 knockdown induces cell necrosis (14). The tumor suppressor p53 has been demonstrated to induce cell cycle arrest and apoptosis, and p53 is a recognized target protein of HPV E6. HPV E6 degrades p53 to promote cell proliferation and inhibit apoptosis (15,16). In autophagy, HPV16 E5 downregulates the expression of key autophagy genes, including Beclin 1, autophagy-related (ATG)5, microtubule-associated protein 1 light chain 3 alpha (LC3), Unc-51 like autophagy activating kinase 1 (ULK1), Unc-51 like autophagy activating kinase 2 (ULK2), autophagy related 4A cysteine peptidase (ATG4A) and ATG7, while E6E7 inhibits autophagosome-lysosome fusion (17). Although there are limited studies available on whether HPV affects pyroptosis, it has also been confirmed that HPV can inhibit the occurrence of pyroptosis (18). HPV is one of the key carcinogenic factors of ESCA, and the majority ESCA samples are also accompanied by HPV infection (8,9). However, the mechanisms through which HPV is involved in ESCA occurrence and particularly the mechanisms through it affects cell death pathways in ESCA, warrant further investigation.

Transcriptomics is the study of the expression and regulation of genes overall. Using transcriptomics to study tumor pathogenesis, explore pathway functions and identify biomarkers for tumor diagnosis and prevention has always been a hotspot in tumor research. In the present study, differentially expressed genes were obtained through the original transcriptomic data of HPV18-positive and HPV-negative ESCC samples. A bioinformatics analysis of differentially expressed genes, combined with experimental verification, was conducted in order to examine the effects of HPV18 on the death pathway of ESCC at the transcriptional level.

Materials and methods

RNA-Seq raw data acquisition and analysis

HPV18-positive and HPV-negative ESCC transcriptomic raw data were downloaded from the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/) through the Accession no. PRJNA530677 (SRX5628804, SRX5628801 and SRX5628804). In these data, HPV18-positive ESCC cells (EC9706 and EC109) and HPV-negative ESCC cells (HKESC01) were derived from a Chinese patient with ESCA (19). Subsequently, the data were processed using Hisat2 (http://www.ccb.jhu.edu/software/hisat/) + Stringtie (http://ccb.jhu.edu/software/stringtie). Read counts were adjusted by a scaling normalization factor using the edgeR package. Differential expression analysis was performed using edgeR (http://bioconductor.org).

Bioinformatics analysis

The Gene Ontology analysis (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed genes were performed using OmicsBean V1.0 (http://www.omicsbean.cn/). The GEPIA2 (http://gepia2.cancer-pku.cn/) website was used to perform the analysis of The Cancer Genome Atlas (TCGA) gene expression and Kaplan-Meier survival analysis associated with the differentially expressed genes in clinical samples of ESCA. The correlation of differentially expressed genes with tumor-associated fibroblast infiltration was performed using the TIMER2.0 database (http://timer.cistrome.org/). Protein-protein interactions were analyzed using STRING (https://string-db.org/).

Cells, cell culture and transfection

Two HPV18-positive ESCC cell lines (EC109 and EC9706) and two HPV-negative ESCC cell lines (KYSE150 and KYSE180) and 293T cells were cultured in DMEM medium (Gibco; Thermo Fisher Scientific, Inc.) with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.), 100 U/ml penicillin and 100 mg/ml streptomycin (HyClone; Cytiva), and placed in a 5% CO2 constant temperature cell culture incubator (Thermo Fisher Scientific, Inc.) at 37°C. The KYSE150 (TCHu236) and 293T (GNHu17) cells were purchased from The Cell Bank of Type Culture Collection of the Chinese Academy of Sciences. The EC109 and EC9706 cells were obtained from Chinese Center for Disease Control and Prevention. The KYSE180 cells was obtained from Shantou University Medical College. The KYSE150, KYSE180, EC109, EC9706 cells originate from the resected specimens of patients with ESCA (19,20). The HPV18 E6E7 gene was obtained from the NCBI gene database (https://www.ncbi.nlm.nih.gov/gene) and synthesized by Tsingke Biological Technology. The HPV18 E6E7 gene was assembled into a pLVX-puro vector, named pLVX-puro-E6E7. The pLVX-puro-E6E7 (8 µg), pMD2.G (1.25 µg) and psPAX2 (5.75 µg) plasmids were transfected into 293T cells using 20 µl jetPRIME transfection reagent (Polyplus-transfection SA) in 10-cm dish for lentiviral packaging. The culture medium was changed after 6 h of transfection and the supernatant was collected at 48 and 72 h, respectively. The supernatants of two times were mixed and concentrated with lentivirus concentration reagent (cat. no. BF06205; Biodragon-Immunotech. Inc.). Following KYSE150 and KYSE180 cell transfection with lentivirus and 8 µg/ml polybrene (cat. no. HY112735; MedChemExpress) for 37°C, 48 h, 8 µg/ml puromycin (cat. no. P8230; Beijing Solarbio Science & Technology Co., Ltd.) was added. The cells stably expressing HPV18 E6E7 were obtained after 10 days of continuous culture, named KYSE150-E6E7 and KYSE180-E6E7. The LC3-RFP-GFP (2 µg) plasmid was transfected into KYSE150 cells and KYSE150-E6E7 cells with 4 µl jetPRIME transfection reagent in a six-well plate for 48 h. The pLVX-puro plasmid was obtained from the Chinese Center for Disease Control and Prevention. pMD2.G and psPAX2 plasmids were obtained from the Institute of Biophysics of the Chinese Academy of Sciences. The LC3-RFP-GFP plasmid was purchased from Tsingke Biological Technology.

RNA extraction and reverse transcription-quantitative PCR (RT-qPCR)

In total, 11 genes were selected, and were enriched in pathways related to apoptosis, autophagy and pyroptosis to verify the NCBI BioProject database PRJNA530677 RNA-seq results acquired through bioinformatics analysis, as aforementioned. PCR primers can be found in PrimerBank (https://pga.mgh.harvard.edu/primerbank/) or were designed using the Oligo7 software according to the instructions (Table SI). Total RNA was extracted from the cells using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.). A PrimeScript RT kit with gDNA Eraser (cat. no. RR047A; Takara Biotechnology Co., Ltd.) was used to synthesize 1 µg of total RNA with a final volume of 20 µl, in order to synthesize the first-strand cDNA. In the process of reverse transcription, the conditions for removing genomic DNA are 2 min at 42°C and 0 sec at 4°C. The conditions for reverse transcription were 15 min at 37°C, 5 sec at 85°C and 0 sec at 4°C. SYBR Premix DimerEraserTM (cat. no. RR091A; Takara Biotechnology Co., Ltd.) was used to detect gene expression with a ViiA7 RT-qPCR instrument (Applied Biosystems; Thermo Fisher Scientific, Inc.). GAPDH and β-actin were used as the internal controls. The KYSE150 and KYSE180 cell lines were used as experimental controls. In total, four sets of each gene were repeated for a 10 µl PCR reaction. The pre-denaturation occurred at 95°C for 30 sec, and then for 40 cycles; each cycle included 5 sec at 95°C, 30 sec at 55°C and 30 sec at 72°C. The melting curve cycle included 0 sec at 95°C, 15 sec at 65°C and 0 sec at 95°C. The 2−ΔΔCq method was used to calculate relative gene expression levels (21).

Apoptosis induction and detection

The apoptosis of KYSE150 cells was induced by cisplatin (DDP; cat. no. HY17394; MedChemExpress). DDP is an antineoplastic chemotherapy agent by cross-linking with DNA and causing DNA damage in cancer cells, it also can activate extracellular-signal-regulated protein kinase (ERK) to cause tumor cell death (22). After treating the KYSE150 cells with 20, 40 and 80 µM DDP, and the KYSE150-E6E7 cells with 80 µM DDP for 24 h, the cells and supernatant were collected and centrifuged at 1,000 × g at room temperature for 5 min, and then were washed twice with PBS (Gibco; Thermo Fisher Scientific, Inc.) and centrifuged again at 1,000 × g for room temperature, 5 min each time. An Annexin V FITC Apoptosis detection kit (cat. no. AD10; Dojindo Laboratories, Inc.) was used to detect apoptosis, and 5 µl Annexin V FITC and 5 µl PI solution were added to each group. Apoptosis was analyzed using BD FACSCalibur flow cytometry (BD Biosciences) and FlowJo v10.6.1 software (https://www.flowjo.com/).

Pyroptosis induction

The metabolites α-ketoglutarate can activate pyroptosis through the GSDMC pathway (23); herein, 15 mM of α-ketoglutarate (cat. no. S30041, Shanghai Yuanye Bio-Technology Co., Ltd.) were added to the KYSE150 and KYSE150-E6E7 cells for 24 h, and the cells were cultured at 37°C under 5% CO2 conditions. Subsequently, the morphology of cells undergoing pyroptosis was observed after 4 h using an inverted fluorescence Axio Observer A1 microscope (Zeiss AG).

Autophagy induction and detection (using western blot anlaysis)

The LC3-RFP-GFP plasmid was transferred into cells. Afterwards, following an 24-h time period, the KYSE150 and KYSE150-E6E7 cells were induced with Hank's Balanced Salt Solution (HBSS; Gibco; Thermo Fisher Scientific, Inc.) for 3 h, and the cells were fixed with 4% paraformaldehyde (cat. no. 441244; MilliporeSigma) and observed and photographed using a laser confocal microscope LSM 5 (Zeiss AG). Following 3 h of HBSS induction, each group of cells was washed twice with ice-cold PBS and a RIPA lysis buffer (Cell Signaling Technology, Inc.) containing protease inhibitor mixture (Roche Diagnostics GmbH) was used for processing. The sample was centrifuged at a rate of 12,000 × g for 10 min at 4°C, and the supernatant was extracted. The BCA protein assay kit (Beijing Solarbio Science & Technology Co., Ltd.) was used to measure the total protein concentration of the sample, and the sample was diluted with PBS to 4 mg/ml. Subsequently, a 5X sodium sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) Sample Loading Buffer (Applygen Technologies, Inc.) was added and boiled in a water bath at 100°C for 5 min. The 40 µg/pore protein sample was added to 10% SDS-PAGE for processing; a semidry method was used to transfer it to a PVDF membrane (MilliporeSigma); 5% skimmed milk powder diluted with TBST was used at room temperature for 1 h for blocking. The membrane with was incubated with the diluted primary antibodies at 4°C overnight. The primary antibodies used included a rabbit anti-LC3B polyclonal antibody (1:1,000; cat. no. EPR18709; Abcam) and rabbit polyclonal antibody GAPDH (1:5,000; cat. no. 10494-1-AP; Proteintech Group, Inc.). Subsequently, the membrane was washed three times with TBST (cat. no. B1009; Applygen Technologies, Inc.) for 10 min each time, and then incubated with goat anti-rabbit IgG secondary antibody (cat. no. 072-06-15-06; KPL, Inc.) with marked Dyelight 680 at room temperature for 1 h. Finally, the Odyssey infrared imaging system (LI-COR Biosciences) was used for detection and quantification was performed using ImageJ 1.52a software (National Institutes of Health).

Statistical analysis and graphing

All the data graphs were constructed using GraphPad Prism 7.04 (https://www.graphpad.com/), Cytoscape V3.7.1 (https://cytoscape.org/) and TBtools v1.108 (24) software. The graphs of TCGA gene expression, Kaplan-Meier survival analysis and tumor-associated fibroblast infiltration were obtained bioinformatics analysis website. All statistical analyses were performed using GraphPad Prism 7.04 and SPSS 22.0 (IBM Corp.). The data results are expressed as the mean ± standard deviation. Comparisons between two groups were made using an unpaired Student's t-test. Comparisons between multiple groups were made using one-way ANOVA followed by Tukey's post hoc test. Kaplan-Meier analysis was performed with the log-rank test. Tumor-associated fibroblast infiltration data was made using Spearman's correlation analysis. P<0.05 was considered to indicate a statistically significant difference.

Results

RCD-related pathways are significantly dissimilar and closely related to classical cancer pathways

By comparing data for HPV18-positive ESCC cells (EC9706) with HPV-negative cells (HKESCC01), 8,962 distinct genes were obtained. Among those genes, 4,177 were upregulated genes and 4,785 downregulated in HPV18-positive ESCC cells (Fig. 1A and Table SII). In total, 9,892 distinct genes were obtained by comparing HPV18-positive ESCC cells, EC109, with the HPV-negative HKESCC01 cells, of which 3,516 were upregulated and 6,376 downregulated (Fig. 1B and Table SIII). Subsequently, 5,908 genes differentially expressed in the two sets of data concurrently, were selected (Fig. 1C). Subsequenlty, 5,393 genes with an absolute value of log2(Foldchange) >2 and padj <0.05 in the two groups of data were selected, of which 2,016 were found to be upregulated and 3,377 downregulated (Fig. 1D).

Figure 1. Gene expression analysis of the HPV18-positive and HPV-negative ESCC samples. (A) Differentially expressed genes between the EC9706 and HKESC01 cells. (B) Differentially expressed genes between the EC109 and HKESC01 cells. (C) Venn diagram of two groups of differentially expressed genes. (D) Expression of genes merged and screened from two groups of data, absolute value of log2(Foldchange) >2, padj <0.05. HPV, human papillomavirus; ESCC, esophageal squamous cell carcinoma.

Omicsbean was then used to analyze the combined and screened differentially expressed genes. KEGG and GO analyses were performed for the differential genes with log2 (fold change)>2 and padj <0.05. In the GO analysis, ‘multicellular organismal process’ (P=0.00e+00), ‘anatomical structure development’ (P=0.00e+00), ‘single-multicellular organism process’ (P=0.00e+00) were defined as the top three most significant pathway in ‘biological progress’ category. ‘Extracellular region part’ (P=1.75e-274), ‘extracellular region’ (P=1.28e-273), ‘extracellular space’ (P=8.20e-238) are the top three most significant pathways in the ‘cell component’ category. ‘Protein binding’ (P=0.00e+00), ‘binding’ (P=8.60e-165), ‘identical protein binding’ (P=4.22e-99) were the top three most significant pathways in the ‘molecular function’ category (Fig. 2A). In the KEGG analysis, the PI3K/Akt signaling pathway (P=4.47e-06), retinol metabolism (P=1.92e-05), cytokine-cytokine receptor interaction (P=2.03e-05), ECM-receptor interaction (P=3.80e-05), protein digestion and absorption (P=6.14e-05), purine metabolism (P=6.43e-05), cAMP signaling pathway (P=1.50e-04), drug metabolism-cytochrome P450 (P=2.08e-04), cell adhesion molecules (P=2.14e-04) and transcriptional misregulation in cancer (P=2.20e-04) were the top 10 most significant pathways (Fig. 2B and C). Of note, the development and progression of cancer involves a large number of significant pathways. A large number of RCD-related pathways were determined and these pathways were also closely related to the top 10 pathways in KEGG (Fig. 2D and E), which may provide evidence concerning the HPV inhibition of cell death, thereby promoting cancer progression. At the transcriptional level, it was suggested that HPV regulates multiple RCD pathways in ESCC.

Figure 2. Pathway analysis of differentially expressed genes. (A) GO analysis of the differentially expressed genes. (B) KEGG analysis of the differentially expressed genes. (C) The top 10 most significant pathways in KEGG analysis (D) Significant RCD-related pathways. (E) The association between RCD-related pathways and the top 10 pathways of KEGG; orange nodes are RCD-related pathways, blue nodes are the top 10 pathways of KEGG. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; RCD, regulated cell death.

HPV18 E6E7 significantly inhibits the apoptosis pathway

There were two significantly different apoptosis-related pathways, in which 18 genes were enriched (Table SIV). Among the enriched genes, three genes were selected [cyclin dependent kinase inhibitor 2A (CDKN2A), plakophilin 1 (PKP1) and desmoglein 3 (DSG3)] with a more significant difference and a clearly reported association with apoptosis (25–27). By validating their expression in HPV18-positive and HPV-negative ESCC cells, the results of RT-qPCR revealed that CDKN2A was highly expressed in HPV18-positive cells, DSG3 expression was comparably reduced in HPV18-positive cells and PKP1 was barely detected in HPV18-positive cells, although it was clearly or highly expressed in HPV-negative cells (Fig. 3).

Figure 3. Reverse transcription-quantitative PCR of apoptosis-related differentially expressed genes. (A) Relative mRNA expression of CDKN2A; data represent the mean ± SD of three independent experiments. (B) Relative mRNA expression of DSG3; data represent the mean ± SD of three independent experiments. (C) Relative mRNA expression of PKP1; data represent the mean ± SD of three independent experiments. *P<0.05 and **P<0.01. CDKN2A, cyclin dependent kinase inhibitor 2A; DSG3, desmoglein 3; PKP1, plakophilin 1.

An analysis of primary tumor samples from GEPIA2 website was also performed, which revealed CDKN2A to be overexpressed, and DSG3 and PKP1 underexpressed in ESCA (Fig. 4A). The correlation analysis of tumor-associated fibroblast infiltration demonstrated that the expression of CDKN2A negatively correlated, and the expression of DSG3 and PKP1 positively correlated with tumor-associated fibroblast infiltration (Fig. 4B). The differential expression of genes in tumors and the association between gene expression and tumor-associated fibroblast infiltration are important indicators for the evaluation the association between genes and the occurrence and development of tumor. These results suggested that the three genes selected were closely related to and may be important for the occurrence and development of ESCA. In the protein-protein interaction (PPI) network analysis, CDKN2A was one of the nodes with the most target proteins (Fig. 4C). This result indicated that CDKN2A may be a critical target of HPV18 oncoprotein and that HPV18 oncoprotein regulates the function of other proteins via CDKN2A.

Figure 4. Correlation analysis of the occurrence and development of ESCA and apoptosis-related differentially expressed genes. (A) The Cancer Genome Atlas differential expression analysis of apoptosis-related differentially expressed genes, data represent the mean ± SD, *P<0.05. (B) Correlation analysis of tumor-associated fibroblast infiltration of apoptosis-related differentially expressed genes; red indicates a positive correlation, blue indicates a negative correlation and gray indicates no significance. (C) Protein-protein interaction network analysis of apoptosis-related differentially expressed genes. ESCA, esophageal cancer.

The oncogenic ability of HPV is mediated by the HPV E6 and E7 oncogenes, which are regularly co-expressed (28). In order to investigate which genes of interest are regulated by HPV E6E7, the HPV18 E6E7 gene was inserted into a pLVX-puro lentiviral vector to stably obtain KYSE150-E6E7 and KYSE180-E6E7 cells expressing HPV18 E6E7 (Fig. S1). The mRNA expression of the CDKN2A, PKP1 and DSG3 genes was detected using RT-qPCR, and the presence of HPV18 E6E7 significantly promoted CDKN2A expression and inhibited the expression of DSG3, having no effect on the transcription levels of PKP1 (Fig. 5A and B). Since the KYSE150 cells are sensitive to DDP, the early and late apoptosis degree of KYSE150 cells was detected at three concentrations of DDP, namely 20, 40 and 80 µM. The results revealed that the early and late apoptosis level of KYSE150 cells was positively associated with DDP (Fig. 5C). DDP at a concentration of 80 µM induced the death of the KYSE150 cells, whereas HPV18 E6E7 significantly inhibited the DDP-induced apoptosis (early and late stage) of the KYSE150 cells (Fig. 5D). These results indicate that HPV18 may regulate apoptosis by affecting the expression of apoptosis-related genes.

Figure 5. HPV18 E6E7 inhibits the apoptosis of ESCC cells. (A) Reverse transcription-quantitative PCR results of CDKN2A, PKP1 and DSG3 in KYSE150 and KYSE150-E6E7; data represent the mean ± SD of three independent experiments. **P<0.01, as compared with the control group. (B) RT-qPCR results of CDKN2A, PKP1 and DSG3 in KYSE180 and KYSE180-E6E7; data represent the mean ± SD of three independent experiments. **P<0.01, as compared with the control group. (C) Detection of apoptosis in KYSE150 cell treated with 20, 40 and 80 µM cisplatin concentrations. (D) Comparison of the effect of cisplatin-induced apoptosis in KYSE150 cells and KYSE150-E6E7 cells. **P<0.01, as compared with the KYSE150 group. HPV, human papillomavirus; CDKN2A, cyclin dependent kinase inhibitor 2A; DSG3, desmoglein 3; PKP1, plakophilin 1.

HPV18 E6E7 significantly inhibits autophagy by inhibiting autophagosome-lysosomal fusion

There are three significantly different autophagy-related pathways, of which 104 genes were found to be enriched (Table SIV). Among the enriched genes, two genes with a more significant difference and a clear reported association with autophagy were selected [ULK1, adrenoceptor beta 2 (ADRB2)] (29,30). The RT-qPCR results demonstrated that the mRNA expression of the ULK1 and ADRB2 genes was reduced in the HPV18-positive cells (Fig. 6).

Figure 6. Reverse transcription-quantitative PCR of autophagy related differentially expressed genes. (A) Relative mRNA expression of ADRB2; data represent the mean ± SD of three independent experiments. (B) Relative mRNA expression of ULK1; data represent the mean ± SD of three independent experiments. ***P<0.001. ULK1, Unc-51 like autophagy activating kinase 1; ADRB2, adrenoceptor beta 2.

ADRB2 expression was also reduced in primary tumor samples of ESCA and the expression of ULK1 was not significantly associated with ESCA (Fig. 7A). The correlation analysis of tumor-associated fibroblast infiltration revealed that the expression of the ULK1 and ADRB2 genes positively correlated with tumor-associated fibroblast infiltration (Fig. 7B). These results suggested that ULK1 and ADRB2 may negatively regulate the occurrence and development of ESCA. In the PPI network analysis, ULK1 was demonstrated to interact with several known autophagy-related genes, including ATG4D, PRKAA2 and WDR45 (Fig. 7C).

Figure 7. Correlation analysis of occurrence and development of ESCA and autophagy related differentially expressed genes. (A) The Cancer Genome Atlas differential expression analysis of autophagy-related differentially expressed genes; data represent the mean ± SD, *P<0.05. (B) Correlation analysis of tumor-associated fibroblast infiltration of autophagy-related differentially expressed genes; red indicates a positive correlation, blue indicates a negative correlation and gray indicates no significance. (C) Protein-protein interaction network analysis of autophagy-related differentially expressed genes. ESCA, esophageal cancer.

The reduced expression of ULK1 and ADRB2 in KYSE150-E6E7 and KYSE180-E6E7 was also detected, revealing that HPV18 E6E7 decreases their mRNA expression levels (Fig. 8A and B). Subsequently, the LC3-RFP-GFP plasmid was transfected into KYSE150 and KYSE150 E6E7 cells and autophagy was induced using HBSS. LC3-RFP-GFP is an effective tool for visualizing the autophagy process (31,32). The simultaneous detection of the green and red fluorescent proteins was used for the detection of LC3 expression. Prior to the fusion of the autophagosome and lysosome, LC3 exhibits red and green signals, and merges into yellow in the autophagosome. Following the fusion of autophagosome with the lysosome, the acidic pH environment weakens GFP, and thus LC3 is presented as a single red dot. In the results of the present study, HPV18 E6E7 inhibited autophagosome-lysosome fusion in late autophagy without hindering LC3-I to LC3-II cleavage (Fig. 8C and D).

Figure 8. HPV18 E6E7 inhibits the autophagy of ESCC cells. (A) RT-qPCR results of ADRB2 and ULK1 in KYSE150 and KYSE150-E6E7; data represent the mean ± SD of three independent experiments. **P<0.01 and ***P<0.001, as compared with the control group. (B) RT-qPCR results of ADRB2 and ULK1 in KYSE180 and KYSE180-E6E7; data represent the mean ± SD of three independent experiments. **P<0.01 and ***P<0.001, as compared with the control group. (C) The autophagy signals of HBSS-induced KYSE150 cells and KYSE150-E6E7 cells were recorded by fluorescence microscopy; yellow arrows indicate autophagosomes and red arrows indicate autophagosome-lysosome fusion. (D) Western blot analysis of the ratio of LC3-II to LC3-I proteins. HPV, human papillomavirus; ESCC, esophageal squamous cell carcinoma; RT-qPCR, reverse transcription-quantitative PCR; ULK1, ADRB2, adrenoceptor beta 2; Unc-51 like autophagy activating kinase 1; LC3, microtubule associated protein 1 light chain 3 alpha.

HPV18 E6E7 inhibits the pyroptosis pathway and α-ketoglutarate-induced pyroptosis

One significantly different pyroptosis-related pathway was determined, enriching four genes through the GO analysis (Table SIV). Pyroptosis is a relatively new concept, leading to the enrichment of incomplete genes. Other pyroptosis-related genes were identified based on a detailed review of the published literature (23,33–37). In total, six genes were obtained [gasdermin A (GSDMA), gasdermin C (GSDMC), NLR family pyrin domain containing 3 (NLRP3), absent in melanoma 2 (AIM2), NLR family pyrin domain containing 1 (NLRP1) and Toll like receptor 1 (TLR1)] for validation through pathway enrichment and a literature survey. The RT-qPCR results revealed that the mRNA expression of the aforementioned six genes was reduced in HPV18-positive cells, with the mRNA expression of NLRP3 being barely detectable in the EC109 and EC9706 cells (Fig. 9).

Figure 9. RT-qPCR of pyroptosis-related differentially expressed genes. (A) Relative mRNA expression of GSDMA; data represent the mean ± SD of three independent experiments. (B) Relative mRNA expression of GSDMC; data represent the mean ± SD of three independent experiments. (C) Relative mRNA expression of NLRP3; data represent the mean ± SD of three independent experiments. (D) Relative mRNA expression of NLRP1; data represent the mean ± SD of three independent experiments. (E) Relative mRNA expression of TLR1; data represent the mean ± SD of three independent experiments. (F) Relative mRNA expression of AIM2; data represent the mean ± SD of three independent experiments. **P<0.01, ***P<0.001. RT-qPCR, reverse transcription-quantitative PCR; GSDMA, gasdermin A; GSDMC, gasdermin C; NLRP3, NLR family pyrin domain containing 3; NLRP1, NLR family pyrin domain containing 1; TLR1, Toll like receptor; 1 AIM2, absent in melanoma 2.

By applying GEPIA2 website analysis, it was defined that only GSDMC mRNA expression was reduced in ESCA tumors. However, there was no marked difference in the expression of GSDMA, NLRP3, AIM2, NLRP1 and TLR1 in ESCA tumors (Fig. 10A). The correlation analysis of tumor-associated fibroblast infiltration revealed that the expression of GSDMC, NLRP3, NLRP1, TLR1 and AIM2 positively correlated, whereas GSDMA expression was not correlated with tumor-associated fibroblast infiltration (Fig. 10B). The expression of TLR1 was also closely related to the survival rate of patients with ESCA, according to an additional a survival analysis (Fig. 10C). The results demonstrated that all genes, apart from GADMA, may be negatively associated with the occurrence and development of ESCA.

Figure 10. Correlation analysis of the occurrence and development of ESCA and pyroptosis related differentially expressed genes. (A) The Cancer Genome Atlas differential expression analysis of pyroptosis-related differentially expressed genes, *P<0.05. (B) Correlation analysis of tumor-associated fibroblast infiltration of autophagy-related differentially expressed genes; red indicates a positive correlation, blue indicates a negative correlation and gray indicates no significance. (C) The Kaplan-Meier survival analysis of TLR1 in ESCC. ESCA, esophageal cancer; TLR1, Toll like receptor 1.

GSDMA and GSDMC are members of the gasdermin superfamily and important executive molecules of pyroptosis (23,33). It was hypothesized that HPV may inhibit the expression of key executive molecules for the inhibition of pyroptosis. The mRNA expression of GSDMA, GSDMC, NLRP3, NLRP1, TLR1 and AIM2 genes was detected in KYSE150-E6E7 cells. The presence of HPV18 E6E7 significantly inhibited the expression of GSDMC and TLR1, having no marked effect on the transcription levels of GSDMA, NLRP1, NLRP3 and AIM2 (Fig. 11A and B). α-ketoglutarate has been demonstrated to induce pyroptosis through the GSDMC pathway (23). Subsequently, 4 h after the induction of pyroptosis by α-ketoglutarate, the HPV-negative group clearly exhibited pyroptosis, whereas in the HPV18 E6E7 group, pyroptosis was limited (Fig. 11C).

Figure 11. HPV18 E6E7 inhibits the pyroptosis of ESCC cells. (A) RT-qPCR results of GSDMA, GSDMC, NLRP3, NLRP1, TLR1 and AIM2 in KYSE150 and KYSE150-E6E7; data represent the mean ± SD of three independent experiments. ***P<0.001 and ****P<0.0001, compared with the control group. (B) RT-qPCR results of GSDMA, GSDMC, NLRP3, NLRP1, TLR1 and AIM2 in KYSE180 and KYSE180-E6E7; data represent the mean ± SD of three independent experiments. **P<0.01 and ***P<0.001, as compared with the control group. (C) Comparison of α-ketoglutarate-induced pyroptosis in KYSE150 cells and KYSE150-E6E7 cells, data represent the mean ± SD of three random visual fields. **P<0.01, as compared with the KYSE150 group. HPV, human papillomavirus; ESCC, esophageal squamous cell carcinoma; RT-qPCR, reverse transcription-quantitative PCR; GSDMA, gasdermin A; GSDMC, gasdermin C; NLRP3, NLR family pyrin domain containing 3; NLRP1, NLR family pyrin domain containing 1; TLR1, Toll like receptor; 1 AIM2, absent in melanoma 2.

Discussion

ESCA is a very aggressive type of tumor with an increasing incidence rate, poor prognosis and Single surgical treatment mode (38). HPV infection is one of the carcinogenic factors of ESCA, and it has also been shown to be negatively associated with the prognosis of ESCA. According to the findings of the present study, HPV may inhibit ESCC cell death and decrease the cell sensitivity to therapeutics, thereby weakening the therapeutic effects, and promoting the development and progression of ESCC.

Using a GO analysis, three types of RCD pathways were obtained. In the apoptosis-related pathway, three genes (CDKN2A, PKP1 and DSG3) were verified. CDKN2A has been alternatively termed in various ways, including ARF, MLM, P14, P16, P19, CMM2, INK4, MTS1, TP16, CDK4I, CDKN2, INK4A, MTS-1, P14ARF, P19ARF, P16INK4, P16INK4A, P16-INK4A. The p16 protein encoded by CDKN2A is usually highly expressed in HPV-related tumors; thus, p16 is an important indicator of HPV infection (39). It has been confirmed that the overexpression of p16 is caused by the deletion of pRB. The function of p16 also requires the synergistic effect of pRB, which is a classical downstream factor of HPV E7. The E7 protein may act on pRB to inactivate it, thereby exerting the carcinogenic effect of HPV (40,41). PKP1 expression is less widespread in lung cancer, with PKP1 having been reported to inhibit cell proliferation, colony formation, migration and invasion and promote apoptosis in lung cancer (27). Concurrently, several studies have demonstrated that PKP1 can inhibit the expression of SPOCK1 (42). Previous studies have confirmed that SPOCK1 can promote the occurrence and development of esophageal squamous cell carcinoma (43). DSG3 is a member of the desmoglein family and cadherin cell adhesion molecule superfamily, which is predominantly expressed in stratified epithelia (44). DSG3 has been reported to exhibit an oncogenic function in head and neck cancer, with DSG3 overexpression promoting tumor growth and migration (45). Among autophagy-related pathways, the present study focused on two genes, ULK1 and ADRB2, which also positively correlated with tumor-associated fibroblast infiltration. In autophagy regulation, ULK1 initiates autophagosome formation, whereas ADRB2 negatively regulates autophagy by disrupting the Akt-dependent Beclin1/VPS34/Atg14 complex (46,47). Members of the gasdermin family play a crucial role in pyroptosis. In the present study, the expression of GSDMA and GSDMC, two genes that are also direct executors of pyroptosis, was significantly reduced. At present, there are limited studies available on GSDMA; however, there have been extensive studies on the induction of pyroptosis by GSDMC. For example, PD-L1 can mediate GSDMC to convert apoptosis into pyroptosis, and α-ketoglutarate can induce pyroptosis through the DDR6-caspase-8-GSDMC pathway (23,48). NLRP1, NLRP3, AIM2 and TLR1 are upstream factors in different pyroptotic pathways (49–51). Through TCGA analysis, it was observed that the genes selected were closely related to the occurrence and development of ESCA. For example, CDKN2A was highly expressed in ESCA and negatively correlated with tumor-associated fibroblast infiltration. DSG3, PKP1, ADRB2, ULK1, NLRP1, GSDMC, NLRP3, AIM2, and TLR1 positively correlated with tumor-associated fibroblast infiltration. ESCC cell lines were then used for further investigation. Through an RNA-Seq analysis and experimental verification, it was confirmed that CDKN2A was highly expressed and DSG3, PKP1, ULK1, ADRB2, GSDMA, GSDMC, NLRP3, NLRP1 and TLR1 are less strongly expressed in HPV18-positive ESCC cells. This result demonstrates that HPV18 regulates these cell death-related genes at the transcriptional level, further promoting the malignancy of ESCC.

The oncogenic ability of HPV is mediated by the HPV E6 and E7 oncogenes (52). Herein, to explore the association between HPV18 E6E7 and these differentially expressed genes, KYSE150-E6E7 cells were constructed and KYSE180-E6E7 stably expressing HPV18 E6E7 using HPV-negative KYSE150 cells and KYSE180 cells. EC109 and EC9706 are HPV18-positive cell lines and were used as negative controls (53,54). Since the KYSE150 cells and KYSE180 cells lack the expression of HPV all subtypes, these two cells are excellent cell models for studying the carcinogenic effect of HPV18 E6E7 on ESCC. Detected at the transcriptional level using RT-qPCR, the presence of HPV18 E6E7 significantly promoted the expression of CDKN2A and inhibited the expression of DSG3, ADRB2, ULK1, GSDMC and TLR1, which may be directly regulated by HPV18 E6E7. The transcription of GSDMA, NLRP3 and NLRP1 was not affected, and these genes may be regulated by other oncoproteins of HPV18 or other factors. In future research, the authors aim to construct EC109 cells and EC9706 cells in which HPV18 E6E7 is knocked out, in order to verify these results and further explore the mechanisms through which HPV18 E6E7 affects ESCC cell death. DDP is one of the commonly used drugs in the treatment of ESCC; however, the resistance of ESCC to DDP has been reported to lead to an unsatisfactory treatment effect (55). In the present study, by using DDP to induce KYSE150 and KYSE150-E6E7 cell death, it was observed that the presence of HPV18 E6E7 significantly inhibited the apoptosis-inducing effect of cisplatin, indicating that HPV18 could increase the resistance of esophageal squamous cell carcinoma cells to cisplatin. Previous studies have revealed that HPV inhibits the occurrence of autophagy in various ways. HPV E5 interferes with the transcriptional activation of autophagy-related genes in the early stage of autophagy, including Beclin 1, ATG5, LC3, ULK1, ULK2, ATG4A and ATG7, HPV E6E7 inhibits the autophagosome–lysosome fusion of cells, thereby inhibiting the process of autophagy (17,56). In the experiments in the present study, it was confirmed that the presence of HPV18 E6E7 inhibited autophagosome-lysosome fusion at the late stage of autophagy in ESCC, having no effect on the cleavage of LC3-I to LC3-II. These results confirmed that HPV18 E6E7 mainly acts in the later stages of autophagy in ESCC cells. HPV18 E6E7 significantly inhibited the expression of GSDMC, and α-ketoglutarate was used to induce the pyroptosis of ESCC cells. The existence of HPV18 E6E7 significantly inhibited the occurrence of pyroptosis of ESCC cells, which may be related to the inhibition of GSDMC expression by HPV18 E6E7. Finally, it was considered that the presence of HPV18 E6E7 significantly inhibits apoptosis, autophagy and pyroptosis in ESCC cells, which is one of the important mechanisms for promoting the occurrence and development of ESCC.

In conclusion, the differentially expressed genes between HPV18-positive and HPV-negative ESCC samples were explored through a transcriptomic analysis, and it was observed that several differentially expressed genes were enriched in various pathways related to cell death. By constructing HPV18-positive ESCC cells, it was confirmed that HPV18 E6E7 may increase the resistance of cisplatin and α-ketoglutarate, inhibits autophagosome-lysosome fusion, thus inhibiting ESCC apoptosis, autophagy and pyroptosis. Finally, it was suggested that HPV may inhibit multiple cell death pathways, which is one of the critical mechanisms for promoting the occurrence and development of ESCC.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was funded by the National Natural Science Foundation of China (grant no. 42077399).

Availability of data and materials

The datasets generated and analyzed during the current study are available in the NCBI BioProject database repository through the Accession no. PRJNA530677 (19). The datasets used and/or analyzed during the current study are also available from the corresponding author on reasonable request.

Authors' contributions

ZZ and DT designed the study and performed the experiments. ZL, GW, QW and BW performed the experiments and data analysis, and contributed significantly to the writing of the manuscript. XH, YZ, MY and CS were responsible for data collection and analysis. ZZ and DT confirm the authenticity of all the raw data. All authors edited the manuscript. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References