Comprehensive Analysis of the PRDXs Family in Head and Neck Squamous Cell Carcinoma

The peroxidase family of peroxiredoxins (PRDXs) plays a vital role in maintaining the intracellular balance of ROS. However, their function in head and neck squamous cell carcinoma (HNSCC) has not been investigated. We therefore explored the value of PRDXs in HNSCC. We found that the expression of PRDX1, PRDX4, and PRDX5 in HNSCC increased while the expression of PRDX2 decreased. Moreover, the high expression of PRDX4/5/6 indicated a poor prognosis. Lower expression of PRDX1/5 was linked to more immune cell infiltration, higher expression of immune-related molecules and a more likely response to anti-PD-1 treatment. Moreover, PRDX5 knockdown inhibited HNSCC cell proliferation, invasion and metastasis and it might promote apoptosis through its antioxidant property. Taken together, our study highlights the potential role of PRDXs in HNSCC. The function of PRDX5 in the development of HNSCC and the formation of the immune microenvironment makes it a promising potential therapeutic target.


INTRODUCTION
Head and neck squamous cell carcinoma (HNSCC) is the sixth most frequent tumor, annually accounting for~830 000 new cases and 430 000 deaths worldwide (1). Despite significant advancements in therapeutic strategies, HNSCC remains a significant health problem as it presents a 5-year overall survival rate of 50% (2). Therefore, more information on the molecular mechanisms responsible for the progress of HNSCC should be obtained to develop more novel therapeutic targets.
The peroxiredoxin (PRDX) family of peroxidases comprises six members (PRDX1 to PRDX6) that constitute a broad range of cellular antioxidant defenders. Based on the number of cysteines, PRDXs are divided into three categories: typical 2-cysteine peroxiredoxins (PRDX1-PRDX4), atypical 2-cysteine peroxiredoxins (PRDX5), and 1-cysteine peroxiredoxins (PRDX6). Although these proteins have different functions in antioxidant protection and cellular redox regulation, they share a common molecular mechanism. PRDXs catalyze the reduction of H2O2, peroxynitrite, and organic hydroperoxides to maintain the intracellular balance of reactive oxygen species (ROS) (3,4). ROS imbalance is closely associated with tumorigenesis and progression. Moderate ROS generation stimulates cell proliferation and differentiation; however, excessive ROS induces oxidative damage to proteins, lipids, and DNA and ultimately leads to cell structure damage and death (5). Thus, the PRDX family has gradually gained attention in the field of cancer research.
According to recent reports, the PRDX family is dysregulated in multiple tumors. PRDX1 is up-regulated in some tumors such as pancreatic ductal adenocarcinoma, nonsmall cell lung cancer, cervical cancer and breast cancer, and promotes the growth and metastasis of cancer cells (6)(7)(8)(9). PRDX2 facilitates the proliferation of colorectal cancer cells and non-small cell lung cancer cells (10,11). PRDX3 and PRDX5, two mitochondrial PRDXs, are associated with unfavorable survival in endometrial cancer and ovarian cancer (12,13). PRDX4 is a metastasis-related marker in oral squamous cell carcinoma, and high PRDX5 expression correlates with a poor prognosis (14). PRDX6 plays a crucial role in mitochondrial dysfunction and ferroptosis (15)(16)(17). However, the role of PRDXs in the occurrence and progression of HNSCC has not yet been elucidated.
Thus, in this study, we focused on the biological roles of the PRDX family in HNSCC. In our results, PRDX5 and PRDX4 were upregulated in HNSCC patients and were associated with a shorter survival time. Patients that showed higher PRDX1/5 expression also had lower immune cell infiltration and lower expression of immune-related molecules. Knocking down PRDX5 in vitro inhibited the proliferation, migration, and invasion of HNSCC and increased apoptosis, which might result from ROS elevation. In addition, patients with low expression of PRDX1/5 was more likely to respond to the PD-1 inhibitor.

Data Source
Level 3 RNASeq, miRNASeq data and the corresponding clinical information of the HNSCC datasets were obtained from The Cancer Genome Atlas database (TCGA) (https://portal.gdc. cancer.gov/). RNASeq data was transformed into transcripts per kilobase million (TPM) values.

Gene Set Variation Analysis and Functional State Signatures
We divided TCGA HNSCC cohort into different groups (PRDXs-high and PRDXs-low) based on the "survminer" R package. To understand the biological pathways in which PRDXs could be involved, we performed Gene Set Variation Analysis (GSVA) using "GSVA" R packages. The gene sets of "hallmark" was downloaded from the "msigdbr" R package. An adjusted P-value of less than 0.01 was considered significant. In addition, CancerSEA (http://biocc.hrbmu.edu.cn/CancerSEA/ home.jsp) was used to evaluate the relationship between PRDXs and 14 functional states based on the GSE103322 data set, including angiogenesis, apoptosis, cell cycle, differentiation, DNA damage, DNA repair, EMT, hypoxia, inflammation, invasion, metastasis, proliferation, quiescence, and stemness.

Tumor Immunology Analysis
We obtained tumor purity via ESTIMATE algorithm (https:// bioinformatics.mdanderson.org/estimate/). We quantified the immune cell infiltration using single-sample gene set enrichment analysis (ssGSEA) based on 23 types of immune cells signatures (18). The differences of immune-related molecules were also assessed in the different groups of PRDXs.

Methylation, Copy Number Variation and lncRNA Analysis of PRDX in HNSCC
Methylation and copy number variation (CNV) of the PRDXs in HNSCC datasets was assessed and visualized based on MEXPRESS (https://mexpress.be/). The relationship between lncRNA and PRDXs was assessed using spearman correlation analysis. |Correlation coefficient| ≥ 0.3 and adjusted P-value < 0.05 were defined as significant.

Immunotherapy Response Prediction
To indirectly predict the immunotherapy response of different expression of PRDX1/5/6, we assessed the similarity of gene expression profiles between PRDX1/5/6 and immunotherapytreated melanoma patients based on subclass mapping (19,20).

Western Blotting
Total proteins of HNSCC cell lines were extracted with RIPA lysis buffer and subjected to 10% sodium dodecyl sulfatepolyacrylamide gel electrophoresis (SDS-PAGE). The proteins were transferred to a 0.22-µm PVDF membrane (Millipore, United States) and blocked with 5% non-fat skimmed milk at RT for 1 h. Next, the membranes were incubated with primary antibodies against PRDX5 (Abclonal, China) and GAPDH (Proteintech, China) at 4°C overnight. After that, the membranes were incubated with species-specific secondary antibodies at RT for 1 hour. Finally, the bands were detected using the enhanced chemiluminescence system (Bio-Rad, United States).

Cell Counting Kit-8 (CCK-8) Assay
HN6 and HSC4 cells with transfected siRNAs were seeded at 1500 cells/200 µl per well in 96-well plates. CCK-8 kit (Beyotime, China) was applied to assess cell viability at different time points according to the protocol of the manufacturer.

Colony Formation Assay
HN6 and HSC4 cells were transfected with siRNAs, seeded in 6well plates (500 cells per well), and cultured for 14 days. Then, colonies were fixed by methanol and stained with 0.5% crystal violet solution, and visible colonies were manually counted.

Wound Healing Assay
HN6 and HSC4 cells transfected with siRNAs were grown in 6well plates and cultured at 37°C. When the cells grew to 90% confluence, a 200 µl tip was used to create a thin scratch. The supernatant cells were washed with PBS. Subsequently, the cells were subjected to a serum-free medium for starvation. The area of the scratch was recorded at 0 h and 24 h with an inverted microscope (ZEISS, Germany).

Transwell Assay
Cell invasion assays were performed in transwell chambers (8 µm pores, Corning, United States) precoated with Matrigel (BD biosciences, United States). In detail, a total of 200 µl serum-free medium containing 8 × 105 cells were added in the upper chambers, and 600 µl medium containing 20% FBS was placed in the lower chambers. After 24 h of incubation, the invading cells were fixed with 4% paraformaldehyde for 30 minutes, stained with crystal violet for 1 h, and counted under an inverted microscope (ZEISS, Germany).

Apoptosis Assay
Apoptosis of HN6 and HSC4 cells was detected by Muse Annexin V & Dead Cell Assay Kit (KeyGEN, China) according to the manufacture protocol. HN6 and HSC4 cells transfected with siRNAs were starved in serum-free medium and resuspended with Annexin V buffer. Then, the cells were incubated with Annexin V-FITC and PI in the dark. A flow cytometer (CytoFlex) was used to assess apoptotic cells.

Mitochondrial ROS (mtROS) Measurement
We employed MitoSOX Red reagent (Thermo Fisher, USA) to determinate mtROS according to the manufacturer's instructions. HN6 and HSC4 cells transfected with siRNA were incubated with 5 µM MitoSox Red dye at 37°C for 10 minutes and then washed with PBS three times. The fluorescence signal of the labeled samples was analyzed via a flow cytometer (CytoFlex).

High Expression of PRDX4 and PRDX5 Predicts HNSCC Prognosis
First, we analyzed the expression levels of PRDXs (PRDX1-6) between HNSCC and normal controls in the TCGA database.
Compared to the normal tissues, PRDX1, PRDX4 and PRDX5 were significantly upregulated, while PRDX2 was significantly downregulated in HNSCC tumors ( Figure 1A). Similar results were found in paired HNSCC and normal tissues ( Figure 1B). Kaplan-Meier survival analysis indicated that higher PRDX3, PRDX4, PRDX5 and PRDX6 expression levels were correlated with poor prognosis ( Figure 1C). Additionally, we found that PRDX3/4/5/6 could be an independent risk factor for the survival of HNSCC (Table 1). High expression of PRDX3/4/5/6 had a significantly higher mortality.

Gene Set Variation Analysis and Functional State Signatures
GSVA analysis was performed to explore the possible roles of PRDXs in HNSCC patients. PRDX1, PRDX5 and PRDX6 were positively associated with DNA repair, MYC targets V1/V2, mTORC1 signaling, oxidative phosphorylation and fatty acid metabolism. Similarity, PRDX1, PRDX5 and PRDX6 were negatively associated with heme metabolism, KRAS signaling up, IL6_JAK_STAT3 signaling, complement, interferon gamma response and inflammatory response. PRDX2 and PRDX3 were positively associated with E2F targets, G2M checkpoint and spermatogenesis. PRDX3 and PRDX4 were negatively associated with heme metabolism, estrogen response, p53 pathway and cholesterol homeostasis ( Figure 2A).
Next, a single cell transcriptome cohort (GSE103322) was used to assess the relationship between PRDXs and 14 functional state signatures in HNSCC tumor cells based on the CancerSEA database. PRDX1 was positively correlated with cell cycle, DNA repair and stemness while PRDX1 was negatively correlated with angiogenesis, apoptosis, differentiation, DNA damage, EMT, hypoxia, inflammation, invasion, metastasis, proliferation and quiescence. PRDX2 was positively correlated with angiogenesis, cell cycle, differentiation, DNA damage, DNA repair, EMT, hypoxia, invasion, metastasis, proliferation and quiescence, while it was negatively correlated with stemness. PRDX3 was positively correlated with apoptosis, cell cycle, differentiation, DNA damage, DNA repair, invasion, metastasis, proliferation and stemness, while it was negatively correlated with angiogenesis, inflammation and quiescence. PRDX4 was positively correlated with angiogenesis, cell cycle, DNA damage, DNA repair, EMT, hypoxia, invasion, metastasis and proliferation, while it was negatively correlated with inflammation and stemness. PRDX5 was positively correlated with angiogenesis, apoptosis, differentiation, EMT, hypoxia, inflammation, invasion and metastasis, while it was negatively correlated with cell cycle, DNA damage, DNA repair and stemness. PRDX6 was positively correlated with cell cycle, DNA damage and DNA repair, while it was negatively correlated with angiogenesis, apoptosis, differentiation, EMT, hypoxia, inflammation, invasion and metastasis, quiescence and stemness ( Figure 2B).

The Relationship Between Immune Microenvironment and PRDXs
Considering that both PRDXs were associated with the immunerelated pathways, the relationship between the immune microenvironment and PRDXs was analyzed. Patients with low PRDX1/5/6 expression showed higher stromal and immune scores than PRDX1/5/6-high patients (Figures 3A, B). PRDX2-low patients gained more stromal score, while a higher immune score was not observed in the PRDX2-low patients.
Thus, we further explored the difference in tumor microenvironment between PRDX1/5/6-low and PRDX1/56high. Consistently with the higher immune score in the PRDX1/5/6-low group, we found increased immune cell infiltration in the PRDX1/5/6-low group (Figures 3C-E).

Methylation, CNV and lncRNA Analysis of PRDX in HNSCC
To explore the mechanisms regulating the expression of PRDXs in HNSCC, we assessed the relationship between PRDXs and methylation, CNV and lncRNA. PRDX1 was negatively correlated with 9 methylation probes (from r = -0.087 to r = -0.317). PRDX2 showed a relatively strong correlation with methylation (from r = -0.186 to r =-0.610). PRDX3 was positively or negatively associated with methylation probes. PRDX5 was negatively correlated with 4 methylation probes (from r = -0.127 to r = -0.177). PRDX6 was also negatively associated with methylation. Only one methylation probe was positively associated with PRDX4. The expression of PRDXs was significantly associated with CNV, especially PRDX6 (r = 0.371) and PRDX5 (r = 0.366) ( Figure 5A). In addition, the correlation between lncRNA and PRDXs was shown in Figure 5B and Table S1.

PRDX1/PRDX5-Low Patients Are More Sensitive to Immune Therapy
Because of the negative correlation between PRDX1/5/6 and immune score, we further explored the response of different PRDX1/5/6 expression to immunotherapy. The results indicated that PRDX1/5-low was more likely to respond to anti-PD-1 ( Figure 6).

Silencing PRDX5 Inhibits HNSCC Cell Proliferation, Migration, and Invasion and Induces Cell Apoptosis
PRDX4 has been reported to promote progression of oral squamous cell carcinoma (14); however, the biological role of PRDX5 in HNSCC is unknown. Firstly, we compared the mRNA levels of PRDX5 in a human normal epithelial cell line (NOK) and in HNSCC cell lines. qPCR confirmed that PRDX5 was upregulated in multiple HNSCC cell lines ( Figure 7A). To elaborate on the biological function of PRDX5, PRDX5 was knocked down using siRNAs in HSC4 and HN6 cell lines ( Figures 7B, C). CCK8 and colony formation assays were used to test the effect of PRDX5 on the proliferation of HNSCC. Silencing PRDX5 in HSC4 and HN6 cells reduced cell viability and colony numbers ( Figures 7D, E). The wound-healing assay showed that silencing PRDX5 significantly decreased cell migration ( Figure 7F). Consistently, the transwell invasion assay showed that the number of cells in the siRNA group that passed through the Matrigel membrane was lower than that of the NC group ( Figure 7G), suggesting decreased invasion. Meanwhile, PRDX5 knockdown provoked apoptosis of HNSCC cells ( Figure 7H). Furthermore, we found that silencing PRDX5 increased mtROS levels in HNSCC cells ( Figure 7I). Thus, these results demonstrated that PRDX5 could promote HNSCC progression through its antioxidant property.

DISCUSSION
Rapidly growing cancerous cells need to maintain high levels of anabolism by ingesting large amounts of nutrients, leading to excessive ROS production and subsequent oxidative stress. To avoid oxidative stress, the expression of antioxidant enzyme systems, such as PRDXs, is tightly regulated (21,22). Although PRDXs has been analyzed in pan-cancer level, the detailed role of PRDXs in prognosis, expression-regulatory relationships, functional state signatures, immune microenvironment and the response to immune treatment in HNSCC is lacking. Thus, we further explored the effect of PRDXs family in HNSCC. In this study, we found that the expression of PRDX1/4/5 increased in HNSCC patients while the expression of PRDX2 decreased. Higher expression of PRDX4/5/6 was correlated with a poor prognosis. In addition, low expression of PRDX1/5 tended to be more likely to response to anti-PD-1 treatment. PRDX2 was associated with prognosis in multiple cancers, such as colorectal cancer (11), hepatocellular carcinoma (23) and esophageal squamous cell carcinoma (24). Though PRDX2 did not predict HNSCC prognosis, we still found that PRDX2 was positively associated with aggressive phenotype, such as angiogenesis, EMT, invasion, metastasis and proliferation, which was in line with oral squamous cell carcinoma (OSCC). Overexpression of PRDX2 could promote OSCC cell proliferation, cell-cycle progression, migration and inhibit apoptosis (25). PRDX3 acts as an oncogene in hepatocellular carcinoma (26), prostate cancer (27) and renal cell carcinoma (28). Knockdown of PRDX3 is associated with decreased ATP synthase and cellular ATP levels, resulting in slower cell growth (26). In our study, we also found that PRDX3 might invoke HNSCC proliferation.
PRDX4 is the only member of the PRDXs family located in the endoplasmic reticulum, and its expression is upregulated in multiple tumors, including OSCC, lung cancer, prostate cancer, colorectal cancer, ovarian cancer and breast cancer (29). We also found that the upregulated PRDX4 might play an oncogene function in provoking proliferation, metastasis, and invasion, leading to a poor prognosis of HNSCC patients, which is confirmed in OSCC (14). The silencing of PRDX4 is related to the generation of ROS in the endoplasmic reticulum, which diffuses to the nucleus and induces DNA damage (30). In addition, PRDX4 can reduce the expression of adhesion molecules (31), thereby increasing the possibility of metastasis and invasion. These may partially explain the involvement of PRDX4 in the proliferation, metastasis and invasion of HNSCC.
PRDX6 is the only 1-cysteine PRDX located in the cytoplasm and translocates to acidic organelles upon its phosphorylation (32). Loss of PRDX6 leads to mitochondrial dysfunction, low growth rate, decreased migration and invasion abilities and ferroptosis (1,15,16,33,34). In our study, patients with higher PRDX6 expression had shorter survival times, suggesting that PRDX6 might act as an oncogene in HNSCC. However, we also found that PRDX6 was negatively correlated with angiogenesis, EMT, invasion and metastasis. Thus, the function of PRDX6 in HNSCC needs to be further explored.
PRDX5 is the only atypical 2-Cys PRDXs, mainly located in mitochondrion. PRDX5 is elevated and drives tumorigenic phenotype in colon cancer (35), non-small cell lung cancer (36), and gastric cancer (37). In our study, we also found that PRDX5 could act as an oncogene in HNSCC. Knocking down of PRDX5 significantly increased mitochondrial ROS levels, promoted apoptosis and suppressed proliferation, migration and invasion of HNSCC cells. The relationship between excess ROS and aggressive phenotypes has been demonstrated in multiple tumors (22). Owing to the relationship between ROS and immune, we also explore the connection between immune microenvironment and PRDXs. Compared with PRDX1/5-high patients, we found an increase in increased immune cell infiltration, HLA molecules, chemokines, cytolytic activity-related genes and immune checkpoint molecules in PRDX1/5-low patients, indicating that PRDX1/5-low exhibited a hot tumor state. Several possible mechanisms may contribute to this phenotype. For one thing, in the PRDX1/5-low group, a high level of ROS can induce the production of chemokines (for example, CXCL9), thereby recruiting T cells into the tumor (38). For another thing, high ROS levels can lead to cellular oxidative damage and immunogenic cell death in the tumor, providing a potential antigenic stimulation for the immune system, and further attract immune cells into the tumor (39). Additionally, tumor-specific MHC-II expression is correlated with the activation of IFNG pathway, higher levels of Th1 cytokines, and CD4 + /CD8 + lymphocytes infiltration (40). Therefore, it is reasonable to find that patients with low PRDX1/5 have a more promising response to anti-PD-1 therapy. Furthermore, PRDX5 actives the Nrf2 signaling pathway, which could protect tumor cells from antitumor immunity (36,41). Therefore, our results highlight the value of targeting PRDX1/5 to inhibit the progression of HNSCC and turn "cold" tumors into "hot" tumors. The box plots of immune checkpoints between PRDX1/5/6 high and PRDX1/5/6 low patients. ns, not significant; *p-value < 0.05; **p-value < 0.01; ***p-value < 0.001; ****p-value < 0.0001.
Although this study systematically analyzed the value of PRDXs and verified the function of PRDX5 in HNSCC, some limitations are acknowledged. The PRDXs expression and survival analysis is mainly based on the data retrieved from the TCGA database, and therefore multi-centered cohorts are required to confirm the results presented here. In addition, it is necessary to elucidate the in-depth mechanism by which PRDX5 promotes the aggressive phenotype of HNSCC.

CONCLUSION
In summary, we found that PRDX1/2/4/5 were dysregulated in HNSCC and that higher expression of PRDX4/5/6 was associated with a shorter survival time. The possible function of PRDX5 in the development of HNSCC and the regulation of the immune microenvironment makes this protein a promising therapeutic target.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
JX and BC designed the study. RC performed all the bioinformatics analysis described here. WZ performed in vitro experiments. RC, HZ, and WZ wrote and edited the manuscript. LW, XC, and XR collected and examined the data. JX and BC supervised the project. All authors read and approved the final manuscript.