Bioinformatics analysis of laryngeal squamous cell carcinoma based on the high infection rate of HPV in Northwest China

Tumor tissue specimen

Paraffin specimens of patients with laryngeal squamous cell carcinoma with complete clinical data and pathological specimens admitted to our hospital from January 2022 to June 2024. The criteria for case inclusion and exclusion were as follows: (1) The postoperative pathological assessment confirmed the diagnosis of primary laryngeal squamous cell carcinoma, with the tumor localized to a single area; (2) complete clinical data for the patient was available; (3) the RNA content of the tissue sample met the required standards upon analysis. The collection procedures were carried out with the informed consent of all patients under a protocol approved by The First People’s Hospital of Yinchuan. Each tissue sample acquired from patients was accompanied by a definitive pathological diagnosis. The study was approved by the Ethics Committee of Yinchuan First People’s Hospital (KY-2024-017) and written informed consent was obtained from all subjects.

PCR-reverse dot blot hybridization

HPV genotyping was detected by PCR in vitro amplification and DNA reverse dot hybridization. By using the specific primers designed with the genetic characteristics of HPV, the target fragments of 23 HPV genotypes can be amplified, and the 23 HPV genotypes can be detected, including 17 high-risk ones: HPV16, 18, 31, 33, 35, 39, 45, 51, 53, 56, 58, 59, 66, 68, 73, 82; six low-risk types: HPV6, 11, 42, 43, 81, 83. Then, the amplified products were hybridized with the typing probes fixed on the membrane strip, including 17 high-risk types and six low-risk types, and the infection of these HPV genotypes was determined according to the presence or non-presence of hybridization signals. Results: According to the location of blue spots on the membrane strip, the genotype information marked at the corresponding location was read: If blue spots appeared at only one genotype, it was a single infection of the corresponding genotype; If blue spots appear at multiple genotypes, it is a mixed infection of the corresponding genotypes. Positive judgment value: 120 clinical HPV-negative specimens were detected, and the gray value (int) of the detection results was calculated, and the positive judgment value was 8.8 int.

Cell culture and treatment

The LC cells TU212 (iCell-h220) was cultivated in IMDM (iCell-0008) appended with 10% fetal bovine serum (FBS, Gibco, Waltham, MA, USA), penicillin/streptomycin (Gibco, Waltham, MA, USA), and humidity of 95% air and 5% CO₂ at 37 °C.

Generation of cells stably expressing HPV16E7

To gain a more targeted and in-depth understanding of the functions of E7, we decided to transfect only E7 in this experiment. TU212 cells were seeded in 6 cm culture dishes at a density of 1.5 × 10⁶ cells per dish and grown until they reached approximately 90% confluence. The recombinant plasmid pCDH-CMV-T2A-EGFP-EF1A-Puro carrying the HPV16E7 gene was then stably transfected into the cells using Lipofectamine 3000 (Invitrogen, Waltham, MA, USA). Twenty-four hours post-transfection, puromycin (2 μg/mL; Invitrogen, Waltham, MA, USA) was added to the culture medium to select for stable HPV16E7-expressing cells (TU212HPV) over a period of 2 weeks. The successful establishment of stable expression was confirmed by RT-qPCR and Western blot analysis.

Reverse transcription‑quantitative polymerase chain reaction (RT‑qPCR)

Total RNA was isolated from TU212HPV cells and TU212 cells containing the empty plasmid (TU212CON) using Trizol Reagent (Thermo Fisher Scientific, Waltham, MA, USA). The RNA sample is mixed with the DNase I buffer and incubated at 37 °C for 15 min. After that, the DNase reaction is terminated by using EDTA. The SuperScript™ VILO™ cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA) was applied to reverse the overall RNA into cDNA according to the instructions. Following that, the mRNA level of HPV16E7 was detected using RT-qPCR. The experiment employed PowerTrack™ SYBR Green Master Mix for qPCR (A46109; Thermo Fisher Scientific, Waltham, MA, USA) and was conducted on the 7900HT Fast Real-Time PCR System (Applied Biosystems, Waltham, MA, USA). ACTB (β-actin) served as the reference gene for data normalization. The changes in the HPV16E7 mRNA level were calculated using the 2^−∆∆Ct method. The primer sequences of HPV16E7and β-actin were offered:

β-actin

Forward Primer: 5′-GGCATCGTGATGGACTCCG-3′

Reverse Primer: 5′-GCTCGTCTTCACGGTTCCAT-3′

HPV16E7

Forward Primer: 5′-ATGACAGCTCAGAGGAGGAGGATG-3′

Reverse Primer: 5′-AACCGAAGCGTAGAGTCACACTTG-3′

Western blot (WB) analysis

The amount of 50 µg of protein was isolated by using SDS-PAGE gel electrophoresis (G2044; Servicebio, Wuhan, China) and shifted to polyvinylidene difluoride (PVDF) membranes (ISEQ00010; Millipore, Billerica, MA, USA). The protein-free rapid blocking buffer (P30500; NCM Biotech, Suzhou, China) was used for membrane sealing for 15 min. The membranes were then incubated with the antibody at 4 °C overnight: Anti-E7 (1:3000, ab20191; Abcam, Cambridge, UK) and anti-tubulin (1:2000, ab179503; Abcam, Cambridge, UK) and the secondary antibody at room temperature for 1 h. The protein band was observed under the HRP ECL system (P10100; NCM Biotech, Suzhou, China), and its grey value was calculated by ImageJ version1.4 (National Institutes of Health, Bethesda, MD, USA).

Transcriptional-level differential gene expression analysis

A total of 20,263 genes were detected. Next is the transcriptional-level differential gene expression analysis. The gene expression values of the transcripts were calculated by StringTie (version 1.3.3b). To identify the differentially expressed genes (DEGs) between two samples, DESeq2 (Love, Huber & Anders, 2014) (version 1.12.4) was employed. A gene was considered to be significantly differentially expressed under two conditions: the q-value was ≤0.001 and the absolute value of the FoldChange was ≥2. When the normalized expression of a gene was zero in both samples, its expression value was adjusted to 0.01. This adjustment was necessary because a value of 0 cannot be properly depicted on a log-scale graph. Additionally, if the normalized expression of a specific gene in both libraries was less than 1, this gene was not included in the subsequent differential expression analysis.

Proteomics quantitative analysis: TMT-based differential protein expression study

ProteomeDiscoverer (v.2.4) was used to search all of the raw data thoroughly against the Homo sapiens Uniprot database. Database search was performed with Trypsin digestion specificity. Alkylation on cysteine was considered as fixed modifications in the database searching. For protein quantification method, TMT/ITRAQ was selected. A global false discovery rate (FDR) was set to 0.01 and protein groups considered for quantification required at least 1 peptide. In this study, mass spectrometry was performed on six samples. A total of 6,971 proteins were identified, among which 6,858 proteins had quantitative information. The thresholds of fold change (>1.2 or <0.833) and P-value < 0.05 (derived from t tests performed in R with FDR correction via the Benjamini-Hochberg (BH) strategy) were used to identify differentially expressed proteins (DEPs).

Functional analysis of differentially expressed genes

Functional enrichment analyses using Gene Ontology (GO) and KEGG were carried out to find significantly enriched GO terms or metabolic pathways in differentially expressed genes (DEGs). GO, an international gene function classification system, annotates gene functions via Molecular Function, Cellular Component, and Biological Process. The KEGG database uses the hypergeometric test to identify enriched metabolic or signal transduction pathways in DEGs against a reference genome. GO terms and KEGG pathways with a false discovery rate (P-value) < 0.05 are considered significant. The DAVID (Sherman et al., 2022) online database was used to analyze the biological, cell-related, and molecular functions of DEGs through GO and KEGG enrichment pathways.

PPI network construction and module analysis

The list of identified proteins was evaluated using version 12.0 of the STRING database (https://cn.string-db.org) to predict the identified proteins. The interaction score was set to a minimum required interaction score: 0.400. Homo sapiens was selected as the organism. Cytoscape 3.10.3 is a freely available platform for network visualization and analysis (Cytoscape_v3.10.3). PPI network construction in Cytoscape requires each protein in the input file to have the same identifiers. Therefore, UniProt ID mapping with Cytoscape was applied to standardize the identifiers. The interactors in the network were carefully evaluated to identify potential interactions between the nodes.

The Cancer Genome Atlas (TCGA) database expression trend verification

To further validate the expression differences of key genes in human papillomavirus (HPV)-positive head and neck squamous cell carcinoma (HNSCC) tissues, HPV-negative HNSCC tissues, and normal tissues, and to conduct survival analysis, we utilized the UALCAN data analysis tool (http://ualcan.path.uab.edu/) to query The Cancer Genome Atlas (TCGA) database. Specifically, we obtained data from 80 HPV-positive HNSCC tissue samples, 434 HPV-negative HNSCC tissue samples, and 44 normal tissue samples. Through these data, we were able to systematically evaluate the expression patterns of key genes across different sample types and their impact on patient survival. Additionally, to more broadly validate the expression differences of these key genes between tumor tissues and normal tissues, we also employed the GEPIA2 online tool (http://gepia2.cancer-pku.cn/#index). GEPIA2 provides a convenient platform that enables researchers to directly investigate publicly available cancer transcriptome data, thereby offering additional support and validation for our findings. This approach not only helps confirm the expression characteristics of candidate genes under different pathological conditions but also provides valuable insights into how HPV status influences the biological behavior of HNSCC. By integrating resources from multiple public databases, our study aims to offer a more comprehensive understanding of the molecular mechanisms and potential therapeutic targets in HNSCC.

Statistical analysis

Statistical analyses were performed using GraphPad Prism 10.3.0 software. The groups’ comparison and multiple groups’ comparison was analyzed by t-test and one-way ANOVA following Tukey’s test, respectively. All experiments were performed in triplicates at least. When the value of P < 0.05, the difference was identified as significance.

HPV infection rate, type distribution and demographic characteristics of laryngeal cancer patients in Ningxia, China

Figure 1 presents the infection type distribution of HPV among 115 patients. The specific data indicates that the number of HPV16 infections is the highest, reaching 57 individuals, which accounts for 49.6% of the total. Subsequently, there are three patients (2.6%) with HPV11 infection, four patients (3.5%) with HPV33 infection, one patient each (0.9%) with HPV18 and HPV35 infection, and two patients (1.7%) with HPV52 infection. Among the 68 patients infected with HPV, one patient was co-infected with HPV16 and HPV18, two patients were co-infected with HPV16 and HPV33, and one patient was co-infected with HPV16 and HPV52. 51 people were not infected. It is evident from the data distribution that HPV16 is the most prevalent HPV type, while the incidences of other types of HPV are relatively low. The overall HPV infection rate of laryngeal cancer in Ningxia, China is 55.65%.

Among HPV-positive laryngeal cancer patients, ages ranged from 42 to 84 years, with median 61, first quartile 55, third quartile 67, and mean about 63.18. Females accounted for 6.25%. For age comparison between HPV-negative and -positive patients, an unpaired t test was done. Results: P-value 0.8232 (ns), two-tailed, t = 0.2239, df = 113, no significant difference as P < 0.05 not met. For gender difference analysis, a Fisher’s exact test was used. P-value 0.6916 (ns), two-sided, no statistical significance as P < 0.05 not satisfied. In conclusion, no significant differences in age or gender were detected between the two groups.

A laryngeal carcinoma cell harboring HPV16E7 was successfully transfected

A laryngeal cancer cell line TU212HPV was established, which stably expresses HPV16E7. As depicted in Fig. 2, the TU212HPV cells, upon successful construction, exhibited green fluorescence when visualized under a laser scanning microscope (Fig. 2A). Moreover, RT-qPCR and Western blot assays were employed to determine the expression levels of E7 mRNA and protein (Figs. 2B, 2C). TU212CON represents the TU212 laryngeal cancer cell line carrying an empty plasmid.

DEGs identification

DESeq2 (version 1.12.4) was used to determine differentially expressed genes (DEGs) between two samples. Genes were considered as significant differentially expressed if q-value ≤ 0.001 and |FoldChange| ≥ 2. Gene expression differences were visualized by volcano plot (Fig. 3A). After differential expression analysis, a total of 1,336 significantly differentially expressed genes (DEGs)were identified (including 797 up-regulated genes and 539 down-regulated genes.) in TU212HPV when compared to TU212.

Ontology analysis of differential genes

Gene IDs of 1,336 differential genes were input into DAVID database for gene ontology analysis. It turns out that over the course of biological process (BP), differentially expressed genes are mainly involved in developmental process, multicellular organismal process, locomotion, response to stimulus and biological process involved in interspecies interaction between organisms. In terms of molecular function (MF), differentially expressed genes are mainly involved in the functions of transcription regulator activity and transporter activity (Fig. 3B).

KEGG analysis of differential genes

Analysis using the DAVID database analysis reveals that the differentially expressed genes primarily participate in pathways such as neuroactive ligand-receptor interaction, cytoskeleton in muscle cells, MAPK signaling pathway, and cytokine-cytokine receptor interaction. Additionally, these genes are involved in pathways known to be associated with cancer development and progression, including transcriptional misregulation in cancer, cytokine-cytokine receptor interaction, MAPK signaling pathway, and JAK-STAT signaling pathway (Fig. 3C). The results of GO and KEGG analyses of the upregulated and downregulated differentially expressed genes have been placed in GO and KEGG Statistical Summary Table of the Supplemental Materials.

Differential proteins identification

In the quantitative proteomic analysis, the differential proteins were screened by specific criteria. Specifically, the screening was based on the quantitative protein ratios and the P-values of statistical tests. I chose the TMT method for protein quantification. The thresholds of fold change (>1.2 or <0.833) and P-value <0.05 were set to identify the differentially expressed proteins (DEPs). Eventually, a total of 236 differential proteins were detected, among which 124 were highly expressed and 112 were lowly expressed (Fig. 4A). The complete raw results of high-throughput sequencing and protein sequencing have been uploaded to the public repository Zenodo (DOI 10.5281/zenodo.15389257).

Ontology analysis of differential proteins

A total of 124 UniProt accession numbers of highly expressed proteins were input into the DAVID database, and 126 DAVID IDs were retrieved for further gene ontology analysis. The results indicate that in biological processes, the differentially expressed proteins are mainly involved in major functions such as the response to stimulus, the biological process involved in interspecies interaction between organisms, the immune system process, detoxification, biological regulation, and the viral process. In terms of molecular functions, the differentially expressed proteins are mainly engaged in catalytic activity functions (Fig. 4B).

By the same method, 112 differential proteins with low expression were analyzed. The results showed that in biological processes, these differentially expressed proteins were mainly involved in localization, cellular process, and developmental process. At the cell component level, no enrichment was observed when p < 0.05. In terms of molecular function, the differentially expressed proteins were mainly involved in cytoskeletal motor activity and structural molecule activity (Fig. 4C).

KEGG analysis of differential proteins

Utilizing the DAVID database to conduct an analysis of the differentially expressed proteins between the two groups, the resultant KEGG enrichment pathways are presented in Figs. 4D, 4E. Upon careful examination and interpretation of these figures, it becomes evident that the proteins exhibiting increased expression levels as a direct result of HPV infection possess a more profound and significant impact on the complex processes of oncogenesis and progression within laryngeal cancer cells. Their roles are likely intertwined with various molecular mechanisms and signaling pathways that contribute to the malignant transformation and advancement of the disease. Given the notable influence of these upregulated proteins, our research focus will be strategically centered and predominantly directed towards a more in-depth exploration and understanding of their functions, interactions, and potential as therapeutic targets.

PPI network analysis

In this study, we focused specifically on the differentially expressed genes with increased expression for further analysis. To deeply investigate the interaction relationships among these genes, 124 such genes were imported into the STRING database, and a PPI network was successfully constructed. This network consisted of 124 nodes (i.e., target proteins) and 244 edges (representing protein-protein interactions). Among them, the degree value was employed to quantify the number of interactions between a protein and other proteins. After calculation, the average degree value of all nodes in this network was 3.94. Subsequently, with the help of Cytoscape 3.10.3 software, nodes with a degree value less than double the average degree value (7.88) were removed. Then, the obtained PPI network was re-screened to eliminate nodes with values less than the average betweenness centrality (BC), average closeness centrality (CC), and average degree value. Eventually, a PPI network composed of 11 nodes and 32 edges was obtained. The 11 key proteins were EGFR, CDC42, paxillin (PXN), SLC2A1, glyceraldehyde 3-phosphate dehydrogenase (GAPDH), FGF2, ICAM1, ITGB1, SFN, PGK1, ISG15 (Fig. 4F). Further KEGG enrichment analysis was carried out on these 11 key proteins, and a total of 14 pathways were obtained (Fig. 4G). Table S1 (Please refer to Table S1 for details.) also presents the genes enriched in each pathway. SFN was not enriched, and among the remaining 10 key proteins, CDC42 was enriched 13 times, ITGB1 was enriched 12 times, EGFR was enriched 10 times, and PXN was enriched nine times. These pathways were mainly involved in proteoglycans in cancer, regulation of actin cytoskeleton, HIF-1 signaling pathway, leukocyte transendothelial migration, and human papillomavirus infection.

Since our research was restricted to the genes with increased expression, we did not conduct the same in-depth analysis on the genes with reduced expression as described in the original text. This was due to the research focus of this project, which aimed to understand the specific roles and interactions of the upregulated genes in the context of the biological process.