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).

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/).

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% CO₂ 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.

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).

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/).

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% CO₂ conditions. Subsequently, the morphology of cells undergoing pyroptosis was observed after 4 h using an inverted fluorescence Axio Observer A1 microscope (Zeiss AG).

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).

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.

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 log₂(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).

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.

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).

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.

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.

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).

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).

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).

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).

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.

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).