Abstract
Background. Head and neck squamous cell carcinomas (HNSC) are common malignant tumors with a high occurrence and poor prognosis. Tumor protein P73 (TP73) plays an integral role in a wide range of human malignancies, but its gene expression profile, prognostic value, and potential mechanisms in HNSC remain to be comprehensively explored. Objective. This research aimed to elucidate the potential relationship between TP73 and HNSC through bioinformatics analysis. Methods. The Cancer Genome Atlas (TCGA) database was queried to investigate the regulatory role of TP73 in HNSC. The survival probabilities linked to TP73 mRNA were determined via the Kaplan-Meier analysis using R packages. Subsequently, the association of TP73 with several clinical subgroups and immunological subtypes was studied using a cohort from the TCGA-HNSC. Functional analyses were used to identify the potential signaling pathways enriched by the correlated genes of TP73. The relationship between TP73 and immunological aspects, including immune cells, immune inhibitor genes, immune stimulator genes, and tumor immune microenvironment, were investigated. Results. This study showed that the protein and mRNA levels of TP73 in HNSC patients were significantly higher than those in normal tissues. Elevated TP73 expression was related to a better survival outcome in HNSC patients. The TP73 gene was an independent prognostic factor for overall survival in HNSC samples. TP73 was mainly involved in DNA replication, ribosome, apoptosis, mismatch repair, and folate biosynthesis. TP73 was found to be positively correlated with the majority of tumor infiltrating immune cells and immunoinhibitory genes in HNSC. Conclusions. Integrative bioinformatics and statistical analyses displayed that TP73 might serve as a novel marker for the diagnosis and prognosis of HNSC. TP73 modulates immune cells in the tumor microenvironment of HNSC patients, thereby bearing significance for HNSC immunotherapy.
1. Introduction
The p53 gene guard family (including TP53, TP63, and TP73 genes) is involved in tumorigenesis and development by coordinating cell proliferation, death, and differentiation [1]. Unlike frequently mutated TP53, TP73 shows a low mutation rate of 0.6% in primary tumors [2]. Generally, TP73 is considered a tumor suppressor and plays a compensatory role in TP53 mutant tumors. At mild DNA damage stages, TP73 can transcriptionally activate TP53 targets and promote cell cycle arrest, while at irreparable DNA damage stages, it stimulates apoptosis or senescence programs [3].
TP73 exists two functional opposable subtypes, TAp73 and DNp73, imbalance of which is often observed in the process of tumorigenesis [4, 5]. TAp73 behaves as a tumor suppressor and exerts pro-apoptotic effects, while DNp73 lacks N terminus transactivation domain and acts as an oncogene [6]. Therefore, TP73 may suppress or promote tumor growth in different kinds of cancers [7]. According to previous studies, increasing TP73 expression predicted a favorable survival outcome in cervical cancer patients [8], but conversely was associated with an aggressive bladder cancer phenotype [9]. Additionally, aberrant TP73 expression has been linked to hematological malignancies and their poor prognosis [10]. Unregulated TP73 was shown to convert fibroblasts into tumorigenic cells and promote the proliferation of hepatocytes and the progression of early hepatocellular carcinoma [11]. Revealing the molecular mechanisms of TP73 expression in tumors will contribute to understanding TP73 roles, specifically its pro-apoptotic function.
Head and neck squamous cell carcinomas (HNSC) consist of a group of malignancies affecting the mucosal lining in different anatomical regions of the upper aerodigestive tract. HNSC is prone to recurrence and metastasis, so the prognosis is still poor, with only approximately 40%-50% of them surviving within 5 years [12]. The occurrence and development of tumors involve various genetic changes and a variety of physiological changes. Accumulating studies have shown that the TP73 gene is closely implicated in tumor biology. In TP53-mutated HNSC cell lines, the TAP73 tumor suppressor molecule was highly expressed in some cell lines and did not exert its tumor suppressor effect [13]. Immunostaining of oral squamous cell carcinoma (OSCC) revealed that TP73 was a valuable biomarker in diagnosing and monitoring high-risk precancerous lesions in the oral epithelium [14]. However, whether abnormal expression of TP73 occurs in HNSC and its specific molecular mechanisms and translational significance have not yet been elucidated.
The use of bioinformatics in tumorigenesis and development has become an essential mode of oncological research. Our study is the first to use various available databases to bioinformatically analyze the TP73 gene in HNSC patients. The abnormal transcription of TP73 gene in HNSC and its relationship to the clinicopathological parameters among HNSC patients were investigated by leveraging the TCGA database. In addition, survival curve analyses were performed to assess the association between abnormal transcription of TP73 and survival prognosis of HNSC patients. The genes correlated with TP73 were analyzed and functional enrichment analyses were carried out to decipher the activities and pathways that may be affected. Our work shows TP73 could serve as a biomarker of prognosis and displayed possible biological roles in HNSC pathogenesis, thereby making it a highly valuable prospective target for therapy in HNSC.
2. Methods
2.1. Expression of p53 Family in HNSC Patients
With the help of cBioportal database (https://www.cbioportal.org/), genomic characteristics of tumors can be examined from the DNA level by researchers. After logging into the cBioPortal web, we selected “Head and Neck Squamous Cell Carcinoma (TCGA, Firehose Legacy, Nature 2015 and PanCancer Atlas)” in the “Query” section. “TP53, TP63, TP73” was used as the input to query the variation characteristics. The resulting frequency, mutation types, and copy number changes could be viewed within Oncoprint. TP53, TP63, and TP73 3D structure can be displayed through the “Mutation” module. Plots were selected to analyze the relationship between TP73 copy number changes and gene expression levels in HNSC.
2.2. Analysis of Abnormal Levels in Three p53 Family of HNSC
HNSC patients with grade 3 HT-seq data were acquired from TCGA database in FPKM format. Based on TCGA-HNSC data, following analysis of 546 samples, which included 502 HNSC tumors and 44 healthy samples, was performed. Briefly, RNAseq data with FPKM were transformed into log2. The mRNA expression of TP53, TP63, and TP73 in HNSC was analyzed and visualized by using ggplot2 package in R. Samples with no clinical information were excluded. We conducted analysis on both unpaired and paired samples accordingly.
2.3. Correlation Coefficient Analysis of p53 Family Genes in HNSC
Pearson’s correlation coefficient analysis was adopted to identify the correlation of p53 family genes in HNSC. The visualization of a heat map was based on the - and -value analyzed via ggplot2 package in R.
2.4. Abnormal TP73 Expression in Pan-Cancer Containing HNSC
UCSC XENA for RNAseq analysis was processed by Toil method for TCGA as well as GTEx [15]. The TP73 mRNA expression was examined across multiple tumors and visualized by using ggplot2 package in R. The RNAseq data with the format of TPM (transcripts per million reads) were log2 converted and further analyzed by using Mann–Whitney test.
2.5. Evaluation of TP73 mRNA Expression in HNSC
ROC analysis was performed to evaluate the diagnostic value of TP73 in distinguishing the clinical status of HNSC (tumor vs. healthy control) [16]. AUC more significant than 0.9 denotes good accuracy. 0.7-0.9 implies a medium level of precision, 0.5-0.7 indicates a poor level of accuracy, and 0.5 suggests an unintentional outcome.
2.6. Correlation between TP73 Expression and Methylation in HNSC
Pearson’s correlation coefficient was used to investigate the correlation between TP73 and TCGA-HNSC methylation levels. The analyses were based on the beta methylation 450 data and contained 520 RNAseq data of HNSC samples. 49 methylation sites in total were evaluated, three of which contained too many missing values to match.
2.7. TP73 Protein Expression Levels in Cancers and Normal Tissues
We typed TP73 and then selected HEAD AND NECK CANCER in Pathology with the help of Human protein atlas (https://www.proteinatlas.org/). Nasopharynx, thyroid, and tonsils were selected as the representative tissue types of HNSC, and the same antibody (CAB002514) were selected. The TP73 protein levels were also analyzed in other types of cancers and various normal tissues.
2.8. The expression of TP73 in HNSC subgroups divided by clinicopathological features
The mRNA expression level of TP73 in these two groups was examined by using a t-test and visualized by box plots. HNSC tumor samples were divided into 2 groups according to the different clinicopathological variables.
2.9. Investigation of TP73 Prognostic Value by Survival Analysis
Survival heatmaps were visualized via GEPIA2. The TCGA-HNSC database was used to gather TP73 mRNA expression, clinicopathological data, and general information. The chi-square, Fisher’s exact test, and Wilcoxon’s rank-sum were applied. The HNSC patients were divided based on the median level of TP73 expression. Categorical variable level frequencies were then calculated for the low and high TP73 gene expression groups.
Survival data on tumor samples were included in the TCGA-HNSC dataset. Kaplan-Meier (KM) curves were used to compare survival for different TP73 patterns via log-rank test. The statistical analysis was performed based on the survival package in R, and KM plots were visualized by using the survminer package in R. Based on TP73 expression levels, tumor data were classified as low- and high-expressing TP73. Specifically, progress-free interval (PFI), overall survival (OS), and disease-specific survival (DSS) events were evaluated.
2.10. Subgroup Survival Analysis
This part determined whether overexpression of TP73 mRNA had a significant effect on the OS outcome of HNSC cases by survminer package (version 0.4.9). The statistical analyses were performed via Cox regression. Statistical analyses were conducted with R to obtain KM plots.
2.11. Relationship between TP73 Expression and Clinical Characteristics
The HNSC samples were divided into TP73 low- and high-expression groups. The relationship between TP73 expression and 16 clinicopathological features was analyzed by using the Chisq test, Fisher’s test, and Wilcoxon's rank-sum test. Table 1 shows the distribution of TP73 low- and high-expression samples in different subgroups of HNSC patients.
2.12. Univariate and Multivariate Cox Regression for Survival Analysis
An analysis of Cox regression was performed to determine a correlation between clinical variables and prognosis. The results were achieved by applying the coxph function in R.
2.13. Forest, Nomogram, and Calibration Plots
We constructed forest plots in R using the ggplot2 package so that hazard ratio (HR), 95% CI, and value were obtained. When comparing two instances of a binary feature, the HR is used to determine the relative risk of death. HR>1 implies higher, whereas HR<1 suggests lower.
Combining TP73 expression values with clinical variables to construct a predictive nomogram for HNSC has not been reported. Using the TCGA-HNSC dataset, we created a predictive plotting by means of combining characteristics and gene expression. We used the rms and survival package to draw nomogram survival plots in R. The prognostic type chosen was OS. TNM stages, histological stage, smoking, alcohol history, radiation, primary outcome neck dissection, lymphovascular invasion, TP73 expression, and other characteristics were studied.
Each of Calibration diagram was carried out via rms and survival packages. Fitting the actual observed fractional probabilities at three-time points was plotted to evaluate the model’s prediction accuracy for the actual prognostic outcomes. The model has achieved perfect predictions if the 1-, 3-, and 5-year solid lines match the 45° ideal diagonals.
2.14. Construction of Protein-Protein Interaction (PPI) and Gene-Gene Interaction (GGI) Network
A PPI network comprising TP73 coexpressed genes was created based on the STRING database (https://string-db.org/, version 11.5). The following are examples of advanced settings: high confidence (0.700), no more than 20 interactors. Analyses of PPI network’s TSV file were performed by using NetworkAnalyzer tool in Cytoscape software. Central genes were identified as those with a degree of more than 10 via CytoHubba plug-in. Based on TCGA-HNSC data, the expression pattern and OS outcomes of top 10 hub genes in HNSC were evaluated with ggplot2 R language package.
A TP73-based GGI network was built using the GeneMANIA webserver (http://genemania.org). The TP73 gene was used as input and the principal functions in the network with the fewest FDR values were chosen. TP73 and its 20 associated genes make up the network. A GGI network was constructed, and then the network photos and reports were saved.
2.15. Genes Significantly Correlated with TP73 in HNSC
We used the stat package in R in order to perform the correlation analysis regarding TP73. Pearson’s correlation analysis was performed to obtain the significantly correlated genes of TP73 in HNSC. The protein-coding genes were selected from the correlated genes [17].
2.16. The Top 10 Associated Genes in TP73 Plotted in a Heat Map
Following the TP73 correlated genes determination, a list of top 10 genes sorted by cor-Pearson value in descending order was acquired, which were considered the top 10 positive genes. And a list of top 10 negative genes was sorted by cor-Pearson value in ascending order. The expression patterns of these total 20 related genes in HNSC samples were drawn by using a heatmap.
2.17. Functional Enrichment Analysis
A clusterProfiler package (version 3.14.3) was performed to identify the function enriched by significantly correlated genes of TP73. The adjusted values were calculated by Benjamini and Hochberg statistical method. By setting a threshold of adj<0.05 and -value<0.20, GO terms (e.g., BP (biological process), CC (cellular component), and MF (molecular function)) and KEGG pathways strongly enriched by the top 200 positive and top 200 negative TP73-correlated genes were identified. Bubble plots were generated only for the top 30 terms with adjusted values if this threshold setting significantly contained more than 30 terms. Otherwise, all words were utilized for bubble mapping. Bubble plots were created via ggplot2 package in R.
2.18. Gene Set Enrichment Analysis (GSEA)
Log2FC (fold change) values of TP73 genes were obtained for GSEA. The differentiated expressed genes (DEGs) dysregulated between HNSC and healthy control samples were identified by using DESeq2 package (version 1.26.0) in R. [18]. GSEA analysis was carried out by using the clusterProfile package in R [19, 20]. Access was obtained from KEGG, WikiPathways (WP), and Reactome (REAC) databases. Significant terms included , -value<0.25, and . 30 functions, including all of the 10 functions with negative NES value, and top 20 functions with the highest positive NES value, were shown in Table 2 and visualized by mountain plots.
2.19. TP73 Expression and Tumor Immune Infiltrating Cells (TIICs) in HNSC
Pearson’s statistical approach investigated the relationship between TP73 and 24 TIICs in HNSC. The GSVA package (version 1.34.0) in the R was used to conduct the analysis [21]. Prof. Gabriela Bindea’s work provided the genetic markers for the 24 TIICs [22]. The lollipop plot showed the relationship regarding TP73 expression as well as 24 TIICs of HNSC tissues. For the TIICs with statistically significant correlation with TP73 in HNSC, scatter plots regarding these TIICs were displayed.
2.20. The Correlation between TP73 and Surface Markers of TIICs in HNSC
We explored the associations between TP73 and immune marker sets of 16 TIICs based on the TCGA-HNSC database [23]. Accordingly, the cor-Pearson and value were analyzed via ggplot2 in R.
2.21. TP73 and Immunomodulator genes in HNSC
23 immunoinhibitor genes and 42 immunostimulator genes were chosen after thorough research of immunomodulator genes in HNSC [24]. Pearson’s correlation coefficient was utilized to examine the correlation between TP73 expression and immumodulator genes. TP73’s association with each of the statistically significant immunomodulator genes was illustrated using scatter plots. The correlation between TP73 and Estimate-Stromal-Immune score was estimated by the estimate package (version 1.0.13).
3. Results
3.1. Current Research Flowchart
Figure 1 depicts a flowchart for this study’s design visualization. Firstly, an investigation of the TP73 gene in TCGA-pan-cancers and HNSC was conducted. Secondly, to assess the TP73 function in pan-cancer and HNSC, survival analyses were conducted, including KM analysis, univariate, multivariate Cox regression analysis, as well as nomogram plots. Thirdly, the PPI and GGI networks were built to identify the interacted protein and genes. In step 4, Pearson’s analysis was used to determine the top ten positively and negatively associated genes. In step 5, we conducted functional enrichment analysis as well as GSEA for investigating the relevant functions of correlated genes. In step 6, we explored the role TP73 played in tumor immunity by investigating TIICs, the tumor immune microenvironment and immunomodulatory genes.
3.2. The Mutation Status and Expression of p53 Family Genes in TCGA-HNSC
The mutation frequency of the p53 gene family in TCGA-HNSC was depicted in an OncoPrint plot (Figure 2(a)). Alteration of the TP53 gene was observable in 70% of HNSC tumor samples (931/1330), whereas alteration of TP63 in 21% and TP73 in 1.3% (276/1330 and 17/1330). According to the mutation type and site, truncating mutation of TP53 and amplification of TP63 and TP73 genes can also be observed. Figure 2(b) shows the 3D structure of the p53 family genes. Figure 2(c) shows that HNSC samples harboring TP73 deletion (e.g., deep deletion, shadow deletion) exhibited lower mRNA expression, compared with the HNSC samples that exhibit diploid TP73. There was a significant correlation between TP73 copy number value, and mRNA expression in HNSC samples (Figure 2(d): Log2 copy number values: =0.23, =2.64e-7; Figure 2(e): capped relative linear copy number values: =0.23, =1.17e-10). This result indicated that mutations and copy number mutations were not the main causes of changes in TP73 mRNA expression among HNSC patients.
(a)
(b)
(c)
(d)
(e)
3.3. Correlation Coefficient Analysis of p53 Family Genes in HNSC
Results of Figures 3(a) and 3(b) were consistent, showing that the difference in TP53 expression in HNSC and normal tissues was not statistically significant. However, TP63 and TP73 expression levels in HNSC were upregulated than that in normal tissues based on the TCGA data. Figure 3(c) presents that a most highly positive correlation was found between TP73 and TP53 in HNSC with an value of 0.414 ( <0.001). TP63 and TP73 were also significantly correlated ( =0.199, <0.001), while no significant correlation was found between the expression levels of TP53 and TP63 in HNSC ( =0.084). Among these three p53 family genes, elevated TP73 expression was strongly correlated with higher OS regarding the HNSC cohort, while TP53 and TP63 genes were not statistically associated with survival outcomes according to the KM survival curves. Figure 3(d) vividly presents that only the expression difference of TP73 gene has a significant survival probability in HNSC, so this research was mainly focused on it.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
3.4. Different TP73 Expression in Pan-Cancer Containing HNSC
TP73 mRNA was assessed by using TCGA RNAseq data (Figures 3(e) and 3(f)). Apart from HNSC, TP73 was also found to be remarkably higher in numerous pan-cancers, for example, adrenocortical carcinoma (ACC), breast invasive carcinoma (BRCA), esophageal carcinoma (ESCA), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), ovarian serous cystadenocarcinoma (OV), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), stomach adenocarcinoma (STAD), thyroid carcinoma (THCA), thymoma (THYM), and uterine corpus endometrial carcinoma (UCEC). However, TP73 mRNA was statistically decreased in a number of malignancies, including pancreatic adenocarcinoma (PAAD), skin cutaneous melanoma (SKCM), and testicular germ cell tumors (TGCT). Unlike the unpaired sample analyses, the paired sample data showed TP73 has no aberrant differences for BRCA, ESCA, LUAD, PAAD, THCA, and UCEC.
Figures 3(g) and 3(h) prove again the significantly increased TP73 gene expression in HNSC tissues ( in unpaired sample analysis and in paired sample analysis). In addition, the ROC analysis was used to assess the diagnostic value of TP73 in HNSC for its diagnostic potential. TP73 could differentiate HNSC samples from normal tissues with moderate accuracy (AUC =0.737, Figure 3(i)).
3.5. DNA Methylation Analysis of TP73 in HNSCC
The relationship between TP73 and methylation of HNSC patients is summarized in Table 3. TP73 was positively correlated with 10 methylation sites such as cg22614891, cg18021902, cg19135761, cg26128092, and cg07174627. Additionally, TP73 was negatively correlated with 14 methylation sites such as cg11504517, cg05924583, cg16741710, cg24073122, and cg20611911. The other 22 methylation sites had no significant correlation with the expression level of TP73.
3.6. TP73 Protein Expression Levels in Cancers and Normal Tissues
HNSC along with several skin cancers showed moderate to strong TP73 nuclear positivity. Several cervical and a few urothelial cancers exhibited weak to moderate immunoreactivity, while other cancer tissues were negative (Figure 4(a)). Noticeably, 50% HNSC patients showed high/medium TP73 expression. Figure 4(b) displays the protein expression of TP73 in various types of healthy tissues. It showed high expression in tonsil squamous epithelial cells and medium in oral mucosa squamous epithelial cells. In consistent with the elevated mRNA expression (Figures 3(g) and 3(h)), TP73 protein was also highly expressed in HNSC but not detected in normal tissues such as nasopharyngeal thyroid and tonsil (Figure 4(c)).
(a)
(b)
(c)
3.7. Relationship between TP73 Expression and Clinical Features
Five clinical variables were significantly associated with the expression level of TP73 gene (Figure 5). The expression of TP73 was significantly higher in the subgroup tumor patients with response to treatment (PR & CR) compared with the subgroup tumor patients without response to treatment (PD & SD, ). The TP73 mRNA expression levels of HNSC samples with higher histologic grade (G3) were higher than that of samples with lower histological grade (G1: , G2: ). The elevated expression of TP73 was observed in the subgroup tumor patients with younger age (≤60) compared with older age (>60, ). In addition, TP73 mRNA expression levels of HNSC patients without lymphnode neck dissection were higher than that of patients with lymphnode neck dissection (). The TP73 mRNA expression levels of HNSC samples with a T2 stage were elevated than that of samples with a higher T stage (T3: P =0.048, T4: ). However, there was no significant relationship between the other 11 clinical variables (i.e., M stage and lymphovascular invasion) and TP73 gene expression.
3.8. Clinical Characteristics of the TCGA-HNSC Patients
According to Table 1, TP73 expression was significantly correlated with histological grade, lymphnode neck dissection, and DSS and PFI events (P <0.05), while no significant relationships were observed in other 14 variables.
3.9. Survival Analyses of TP73 in Pan-Cancer and Particularly HNSC
Figures 6(a) and 6(b) present the survival heatmaps, showing the prognostic impact of TP73 in various cancers. In ACC, KIRC, and LGG, elevated TP73 was associated with a worse OS, while in ACC, KIRC, LGG, PAAD, and THCA, elevated TP73 was associated with a worse DFS. Noticeably, the upregulated TP73 was associated with a better OS and DFS survival outcomes in HNSC patients.
(a) The prognostic value of TP73 in pan-cancer (OS)
(b) The prognostic value of TP73 in pan-cancer (DFS)
(c) The prognostic value of TP73 in HNSC
(d) Subgroup survival analysis
Figure 6(c) indicates the prognostic values of TP73 gene for HNSC patients. There was a significant difference in survival time distribution of three prognostic types (OS, PFI, and DSS) according to different expression groups of TP73, which suggested the overexpressed TP73 may predict a better OS outcome.
3.10. Subgroup Survival Analyses
Figure 6(d) suggests that high expression of TP73 did not mean better survival for patients with lower T stages (T1 and T2), while for patients with higher stages (T3 and T4), high expression of TP73 represented greater survival probability. Since there were insufficient samples in M1 subgroup, M stage was not used for performing the subgroup survival analyses. The upregulation of TP73 had a significant impact on OS outcomes in HNSC cases with a higher N stage (N2 and N3, ), clinical stage (Stage III and Stage IV, ), and histologic grade (G3 and G4, ). However, upregulation of TP73 did not exert a statistically significant OS outcome in HNSC cases with a lower N stage (N0 and N1, ), clinical stage (Stage I and Stage II, ), and histologic grade (G1 and G2, ), respectively.
3.11. Visualization of Univariate and Multivariate Cox Regressions via Forest Plotting
In Figure 7(a), univariate analyses results showed that various characteristics are risk factors for death outcome in HNSC patients, including M1 stage (), without receiving radiation therapy (), primary therapy outcome of PD & SD (), and with lymphovascular invasion (). Noticeably, an HR of 1.480 for low TP73 expression HNSC patients indicated that the group had a higher risk of death than those high TP73 gene expression ().
(a) Forest plot
(b) Nomogram plot
(c)
(d)
(e)
According to multivariate analyses results, a number of characteristics, for example, high N stage (N1 and N2 and N3) (), not receiving radiation therapy (), and primary therapy outcome of PD and SD (), influenced OS outcomes in HNSC patients. However, the level of TP73 expression did not make a difference in survival outcomes.
3.12. Conduction of a Prognostic Model for Risk Estimation
Figure 7(b) predicts the survival prediction for HNSC samples, based on the TP73 gene as well as clinical factors. The deviation correction line in the 1- and 3-year OS calibration plotting is close to the ideal curve, indicating that the prediction results agree well with the observation results (Figures 7(c) and 7(d)). Noticeably, the predicted result did not agree well with actual survival results according to the calibration curve for predicting 5-year OS (Figure 7(e)).
3.13. PPI Network Plotting
A PPI network of TP73 expressions with 21 nodes and 88 edges was generated, and it has an interaction score >0.7 according to the STRING online database (Figure 8(a)). Apart from TP53, TP73 was also found to interact with other proteins, for instance, MDM2 (MDM2 Proto-Oncogene), YAP1 (Yes1 Associated Transcriptional Regulator), and EP300 (E1A Binding Protein P300). Figure 8(b) displays the expression and prognostic values of top 10 hub genes in HNSC. Among them, 8 hub genes were significantly upregulated in HNSC, including RB1, CDK1, CCNA2, CDK2, MDM2, EP300, RPS27A, and CCNE1. Regarding prognostic values, only the MYC expression level was significantly related to OS outcomes in HNSC (). Elevated MYC expression presented a worse survival probability for HNSC patients.
(a) Protein-protein interaction network
(b) The expression and prognostic values of 10 hub genes in HNSC
(c) Gene-gene interaction network
3.14. GGI Network Plotting
According to Figure 8(c), GGI network contained TP73 gene and its 20 interacted genes. These genes contained CABLES1 (Cdk5 and Abl enzyme substrate 1), TP63 (tumor protein p63), COP1 (COP1 E3 ubiquitin ligase), TP53 (tumor protein p53), YAP1, IGBP1 (immunoglobulin binding protein 1), MDM4 (MDM4 regulator of p53), PPP1R13B (protein phosphatase one regulatory subunit 13B), FBXO45 (F-box protein 45), RCHY1 (ring finger and CHY zinc finger domain containing 1), SMAD2 (SMAD family member 2), PFDN5 (prefoldin subunit 5), ABL1 (ABL proto-oncogene 1, nonreceptor tyrosine kinase), RNF144B (ring finger protein 144B), SIAH1 (siah E3 ubiquitin protein ligase 1), E2F1 (E2F transcription factor 1), EP300 (E1A binding protein p300), DENND4A (DENN domain containing 4A), CHEK2 (checkpoint kinase 2), and HMGB1 (high mobility group box 1). The functional enrichment analysis showed that these 20 interacted genes were significantly enriched in several apoptotic signaling pathways, such as regulation of protein insertion into mitochondrial membrane, protein insertion into mitochondrial membrane, positive regulation of mitochondrial outer membrane permeabilization, intrinsic apoptotic signaling pathway, regulation of mitochondrial outer membrane permeabilization, signal transduction by p53 class mediator, and regulation of cell cycle arrest.
3.15. Heat Map to Visualize the Top Correlated Genes Expression Pattern of TP73 in HNSC
The top correlated genes ranked by the ascending and descending order of cor-Pearson values were obtained and their expression patterns in HNSC were shown in a heatmap. (Figure 9(a)). The top 10 positively correlated genes of TP73 were listed as follows: TXLNA (Taxilin Alpha), SAP30L (SAP30 Like), PRDM15 (PR/SET Domain 15), CDKN2C (Cyclin Dependent Kinase Inhibitor 2C), TSC2 (TSC Complex Subunit 2), SSBP3 (Single Stranded DNA Binding Protein 3), EDARADD (EDAR Associated Death Domain), TCEANC2 (Transcription Elongation Factor A N-Terminal And Central Domain Containing 2), CYB5RL (Cytochrome B5 Reductase Like), and FAM53B (Family With Sequence Similarity 53 Member B). The top 10 negatively correlated genes of TP73 were listed as follows: UPP1 (Uridine Phosphorylase 1), FOSL1 (FOS Like 1, AP-1 Transcription Factor Subunit), TM4SF19 (Transmembrane 4 L Six Family Member 19), PPIF (Peptidylprolyl Isomerase F), RPL22L1 (Ribosomal Protein L22 Like 1), PTS (6-Pyruvoyltetrahydropterin Synthase), KLK8 (Kallikrein Related Peptidase 8), ATP5MPL (ATP Synthase Membrane Subunit J), S100A10 (S100 Calcium Binding Protein A10), and NQO1 (NAD(P)H Quinone Dehydrogenase 1).
(a)
(b)
(c)
(d)
(e)
3.16. Bubble Charts to Identify the Biological Functions of the Significantly TP73-Correlated Genes
Bubble charts (Figure 9(b)) and Table 4 display the results concerning three GO term aspects (GO-BP, -CC, and -MF) and KEGG pathways. The significantly correlated genes of TP73 were primarily enriched in particular apoptotic-related BP, for example, DNA replication, G1/S transition of the mitotic cell cycle, and translational initiation (Figure 9(b)). The significantly correlated genes of TP73 were primarily enriched in the following ribosome-related CC: cytosolic part, ribosome, ribosomal subunit, and cytosolic ribosome (Figure 9(c)). The significantly correlated genes of TP73 were primarily enriched in certain MF, such as structural constituent of ribosome, catalytic activity (acting on DNA), oxidoreductase activity (acting on CH−OH group of donors), oxidoreductase activity (acting on the CH−OH group of donors, NAD or NADP as acceptor), and oxidoreductase activity (acting on NAD(P)H, Figure 9(d)). Five KEGG pathways were enriched, including DNA replication, ribosome, apoptosis, mismatch repair, and folate biosynthesis (Figure 9(e)). TP73 was reported to be related to some gene families, such as the Minichromosome Maintenance Complex family (MCM2, MCM3, MCM5, and MCM6), Ribosomal Protein Small (RPS12, RPS17, RPS25, and 29), and Large (RPL13, RPL24, RPL26, RPL27, RPL36, and RPL39) gene family (Figure 9(e)).
3.17. GSEA Analysis of TP73-Correlated Genes Functional Terms
According to mountain plots (Figure 10) and Table 2, cell cycle signaling was significantly enriched, including Mitotic G1 phase and G1s transition, G2 M checkpoint, and prometaphase and DNA replication pathways (Figure 10(a)). TP73 significantly negatively correlated genes were principally enriched in biological oxidation, insulin signaling pathway, ion channel transport, and RAS signaling pathways (Figure 10(b)).
(a)
(b)
3.18. Exploration of the Correlation between TP73 Expression and Immunity
To further explore the relationship betweenTP73 gene expression and immunity, we analyzed TIICs, tumor microenvironment, and immune-related genes. Figure 11(a) shows that TP73 had a substantial positive correlation with several TIICs, including T helper cells, Treg, Th2 cells, aDC, T cells, Th17 cells, TFH, Tcm, NK cells, B cells, NK CD56dim cells, Cytotoxic cells, and Tem. Additionally, there was a significant negative correlation between TP73 and neutrophils, a TIIC. For the TIICs strongly correlated with TP73 ( value <0.01), scatter plots were depicted (Figure 11(b)).
(a)
(b)
(c)
(d)
(e)
Table 5 showed the relationship between TP73 and marker sets of 16 TIICs in HNSC. Taken Tumor Associated Macrophage (TAM) as an example, the expression of TP73 was significantly positively correlated with a surface marker of TAM—chemokine factor CCL2 ( =0.148, ), and another surface marker of TAM—cytokine IL10 ( =0.134, ). The highest positive correlation was observed between TP73 expression level and a surface marker of Th17 cells—STAT3 ( =0.463, ).
Figure 11(c) additionally explores the relationship between TP73 and the tumor immune microenvironment. The scatter plots found that stromal and estimate scores show no statistical significance (). However, TP73 was positively correlated with the tumor purity of HNSC according to the immune score ( =0.163, ). In addition, 18 immunoinhibitor genes were positively correlated with TP73, including AORA2A, BTLA, CD160, CD244, CD98, CSF1R, CTLA4, HAVCR2, IDO1, IL10, IL10RB, KDR, KR2DL3, LAG3, LGALS9, PDCD1, TGFBR1, and TIGIT (Figure 11(d)). 6 immunostimulatory genes, CD276, HHLA2, IL6, NT5E, PVR, RAET1E, and TNFSF9, were negatively correlated with TP73 (Figure 11(e)).
4. Discussion
As shown in Figure 2(a), the frequency of TP73 genetic mutation was noted only in 1.3% of HNSC tumor samples, which was significantly lower than that for TP53 and TP63. In accordance with our findings, a previous research revealed the infrequent mutation of TP73 in HNSC with a C→A transversion mutation at codon 469 in C-terminal domain [25]. The combination of TP73 exon 2 G4C14-to-A4T14 and p53 intron three variant alleles was found to provide a significantly augmented risk of HNSC (OR =2.22, 95% CI: 1.08-4.56) in a sample drawn from the Italian population. Among participants less than 45 years old, individuals harboring TP73 exon 2 G4A variant allele had a 12.85-fold increased risk of HNSC when compared with people with the homozygous wild-type genotype (95% CI: 2.10-78.74) [26].
Our current research demonstrated that TP73 gene was overexpressed in HNSC tumor tissue as compared to control samples (Figures 3(g) and 3(h)). Previous immunohistochemistry results also identified an elevated expression of TP73 protein in HNSC tumor samples [27]. The current investigation demonstrated that the TP73 expression level was uncorrelated with the M stage representing metastatic status, or with lymphovascular invasion (Figure 5). In a contradictory finding, a HNSC patient cohort revealed a significant correlation between TP73 expression and distant metastasis and perineural/vascular invasion, suggesting an association between TP73 and an aggressive phenotype [27].
As a member of the TP53 family, TP73 has displayed different expression patterns and has diverse associations with prognosis in most human tumors [28, 29]. Previously, TP73 was proved to be an independent survival indicator in some cancers, such as lymphomas, leukemia, and gliomas. In the present study, increased TP73 indicated a better prognostic outcome for HNSC patients (Figure 6(c)). Similarly, downregulation of TP73 was associated with a worse OS in hematological malignancies of myeloid origin [10]. On the contrary, upregulated TP73 was found significantly associated with worse OS and DFS outcomes in laryngeal cancer [30]. Ribeiro showed that OSCC patients who had altered copy number of TP73 genes had a poor prognosis [31]. In gliomas, increased TP73 served as a stand-alone high-risk factor affecting the prognosis, while the downregulation of TP73 was closely related to favorable OS outcome [32].
PPI network analysis identified the top 10 hub genes which interacted with TP73, and eight of these were found significantly upregulated in HNSC. The expression levels of two genes, TP53 and MYC, showed no significant association with HNSC samples (Figure 8(b)). The interplay between TP73 and several associated genes (TP53, MYC, RB1, E2F1, and PPP1R13B) has been evidenced in previous literature. In TP53 mutant tumor cells, TAp73 could initiate the downstream pathway of TP53 and played a compensatory tumor suppressor role [33]. Serving as an oncogene, MYC may enlist TP73 to induce apoptosis in TP53-deficient tumor cells [34]. In addition, recent research indicated that TP73 inhibition was necessary for the anticancer effects of cyclin-dependent kinases RB1, thereby further leading to the cell cycle transition from G1 to S-phase [35]. Aberrations of cell cycle control were observed in HNSC, which resulted in downregulated activity of E2F transcription factors with concomitant enhanced cell cycle progression. TP73 activation via deregulated E2F1 activity may also serve as an anti-tumorigenic safeguard mechanism in HNSC independent of TP53 [36]. Another interacting gene shown in Figure 8(c), PPP1R13B, was found to enhance apoptotic function of TP63 and TP73 by selectively inducing the expression of endogenous TP53 target genes and further inhibiting tumor growth [37].
The results obtained by functional enrichment analysis and GESA analysis (Figures 9 and 10) showed that TP73 gene is mainly enriched in apoptosis-related pathways (e.g., DNA replication and cell cycle), insulin signaling pathway, Ras signaling, BARD1 pathway, regulation of p53 activity, biological oxidation, and ion channel transport. Functionally, TP73 protein can be divided into two subtypes, TAp73 and DNp73. The identical TAp73 isoforms produce markedly different effects on apoptosis and cell cycle arrest and were shown to act independently of TP53, thereby determining its transcriptionally activated apoptosis and growth arrest pathways. For example, TAp73-transfected osteosarcoma cells exhibited high levels of apoptosis induction at the mitochondrial level. In contrast, TAp73-transfected lung cancer cells exhibited high levels of late apoptosis, which was associated with high expression of death receptor pathway genes. In addition, TAp73 affected the cell cycle distribution of induced p21WAF1 mRNA in lung cancer cells but not in osteosarcoma cells [38]. A previous study regarding the involvement of TP73 in regulating reactive oxidative stress (ROS) showed that increased ROS observed in TAp73 knockout mice might result in the accumulation of mutations in proto-oncogenes, consequently contributing to increased genomic instability [39]. Deletion of the DNp73 played a key role in metabolic reprogramming and regression of p53-deficient malignancies, which was achieved by upregulating the IAPP gene. IAPP acted through the calcitonin receptor and RAMP3 to inhibit glycolysis, leading to ROS production as well as apoptosis. Another study proved that DNp73 synergizes with cMyc to promote the proliferation of primary mouse embryonic fibroblasts (MEFs) and inhibit P53-dependent apoptosis of MEFs. Moreover, an interaction between DNp73 and oncogenic ROS-induced fibrosarcomas in nude mice has been shown [40].
The current research also showed that TP73 was significantly positively correlated with TIICs, such as T cells, B cells, and NK cells (Figure 11(a)). Such findings are well supported by previous research. In nontumor lymphocytes, TP73 function has been found to be essential for antigen-induced circulating peripheral T cell death following the activation of T cell receptors in thymocytes [41]. Overexpression of TP73 is a common feature of B-cell chronic lymphocytic leukemia and may be involved in tumorigenesis by altering the ratio between oncogenic and anticancer gene variants of TP73 [42, 43]. NK cell malignancies exhibit specific promoter methylation patterns in which TP73 is consistently involved. These results suggest that TP73 may be an essential target in NK cell tumor transformation, and analysis of its methylation pattern may be a possible molecular tool for the detection of NK cell lymphoma [38]. To our knowledge, there is no report investigating the interplay between TP73 and tumor immune cells in the context of HNSC; thus, the present data provide a novel research direction for future studies.
It is noteworthy to highlight the limitation of the current research. Although we have designed a variety of bioinformatic analyses to illustrate the implication of TP73 in HNSC, we did not perform experimental validation of the findings predicted by this analysis. Future research can include cellular experiments using genetic transfection assays to observe the alteration of TP73 expression in HNSC cell lines. In addition, subsequent research could include the establishment of a subcutaneous tumor model in rats in order to examine the role of TP73 in the proliferation and migration of tumors in vivo. Another possible research direction is a clinical study to validate the correlation between TP73 expression and clinicopathological factors in HNSC. Taken together, our research provides a solid theoretical basis for the implications of TP73 in HNSC and indicates research directions for future studies.
5. Conclusion
The findings identified upregulated TP73 as a valid prognostic predictor for HNSC. Furthermore, the anti-tumor role of TP73 gene in HNSC was linked to mainly five pathways including DNA replication, ribosome, apoptosis, mismatch repair, and folate biosynthesis. This bioinformatics study highlights a possible biomarker role of TP73 and its therapeutic potential for treating HNSC.
Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The author declares that he has no competing interests.
Authors’ Contributions
Yuming Chen (Email: [email protected]; [email protected]) conceptualized the research idea, performed the study design, conducted the bioinformatics analysis, interpreted the results, and wrote the manuscript. The author read and approved the final manuscript.
Acknowledgments
The author acknowledges the support of his colleagues and the valuable insights and suggestions concerning this study offered by the reviewers.