Upregulation of AKR1C1 and AKR1C3 expression in OPSCC with integrated HPV16 and HPV‐negative tumors is an indicator of poor prognosis

Different studies have shown that HPV16‐positive OPSCC can be subdivided based on integration status (integrated, episomal and mixed forms). Because we showed that integration neither affects the levels of viral genes, nor those of virally disrupted human genes, a genome‐wide screen was performed to identify human genes which expression is influenced by viral integration and have clinical relevance. Thirty‐three fresh‐frozen HPV‐16 positive OPSCC samples with known integration status were analyzed by mRNA expression profiling. Among the genes of interest, Aldo‐keto‐reductases 1C1 and 1C3 (AKR1C1, AKR1C3) were upregulated in tumors with viral integration. Additionally, 141 OPSCC, including 48 HPV‐positive cases, were used to validate protein expression by immunohistochemistry. Results were correlated with clinical and histopathological data. Non‐hierarchical clustering resulted in two main groups differing in mRNA expression patterns, which interestingly corresponded with viral integration status. In OPSCC with integrated viral DNA, often metabolic pathways were deregulated with frequent upregulation of AKR1C1 and AKR1C3 transcripts. Survival analysis of 141 additionally immunostained OPSCC showed unfavorable survival rates for tumors with upregulation of AKR1C1 or AKR1C3 (both p <0.0001), both in HPV‐positive (p ≤0.001) and ‐negative (p ≤0.017) tumors. OPSCC with integrated HPV16 show upregulation of AKR1C1 and AKR1C3 expression, which strongly correlates with poor survival rates. Also in HPV‐negative tumors, upregulation of these proteins correlates with unfavorable outcome. Deregulated AKR1C expression has also been observed in other tumors, making these genes promising candidates to indicate prognosis. In addition, the availability of inhibitors of these gene products may be utilized for drug treatment.


Introduction
Worldwide, head and neck squamous cell carcinoma (HNSCC) is the sixth most common and one of the most lethal types of cancers. Patients frequently develop unfavorable distant metastases and inoperable local and regional recurrences. 1 Even though improvements in detection and clinical treatment (including surgical intervention, radiation and chemotherapy) have been achieved, prognosis of this malignancy with only 40-50% of patients surviving the next 5 years after initial diagnosis remains unfavorable. The most prominent risk factors for the development of HNSCC are excessive tobacco and alcohol use, as well as high-risk human papillomavirus (HPV) infections.
Most frequently, HPV16 has been found to be associated with oropharyngeal squamous cell carcinoma (OPSCC). 2,3 This subgroup of carcinomas shows clinicopathological and molecular characteristics that differ from alcohol-and tobacco-induced carcinomas. 4 A question that still remains to be elucidated is if there is a biological consequence of viral integration. So far, there is little evidence that viral integration may have an impact on prognosis, 5-10 furthermore, in previous work of our group no evidence was found that HPV integration significantly affects viral gene expression or the expression of disrupted human genes as compared to tumors with episomal virus. 3,[11][12][13] Thus, it remains to be identified by which mechanisms HPV integration affects its host cell and if they are of clinical relevance.
Our aim was to perform a genome-wide screen to identify human genes whose expression is influenced by integration of HPV16. For this purpose, we compared HPV16-positive OPSCC harboring extrachromosomal or integrated virus DNA using mRNA microarray expression profiling.
We identified a unique signature of differentially expressed human mRNAs in relation to viral physical state. Selected candidate genes comprised the AKR1C1 and AKR1C3 genes which expression was independently confirmed by RT-qPCR as well as immunohistochemistry using both part of the screening cohort and a new cohort of 141 OPSCC.  Table 1). This cohort was selected out of a previously published group of 75 patients based on the after inclusion criteria: ≥70% tumor cells, sufficient high-quality DNA and RNA (RQI value ≥7; BioRad Experion System, BioRad Laboratories Munich, Germany) and positive HPV-PCR and p16 INK4A Immunostaining. 3 These 33 patients were also analyzed by APOT and / or DIPS PCR for viral integration status, which revealed 9 (27.3%) cases with only integration, 4 (12.1%) cases with both integrated and episomal virus and 20 (60.6%) cases with episomal virus. 3 Cohort 2: Formalin-fixed, paraffin embedded (FFPE) tissue samples from 141 patients with primary OPSCC (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) containing both tumorous and adjacent tumor free tissue were available from the Departments of Otorhinolaryngology, Head and Neck Surgery, and Pathology, University of Cologne. HPV-status (HPV-PCR and p16 INK4A Immunostaining) was available for all 141 patients. 93 patients (66.0%) were HPVnegative, 48 (34.0%) HPV-positive (Table 1).

Ethics statement
Patient material was used according to the code for proper secondary use of human tissue. The ethics committee of the Medical Faculty of the University of Cologne approved our study (approved protocol no. 11-346). Written, informed consent had been obtained from all patients.

Identification of HPV16 integration status
To identify viral integration status in Cohort 1, both a mRNA based assay (Amplification of Papillomavirus Oncogene Transcripts PCR; APOT-PCR) and a DNA based assay (Detection of Integrated Papillomavirus Sequences-PCR; DIPS-PCR) were performed as described previously. 3

mRNA expression profiling
Total RNA was collected from a subset of 33 samples of appropriate mRNA quality, selected from 75 patients reported in a previous study (Cohort 1 3 ). Samples were analyzed using Agilent Whole Human Genome 4 × 44 K Microarrays, which represent more than 41,000 unique human transcripts. Labelling and hybridizations were performed according to the manufacturer's instructions (Agilent Technologies). Hybridized arrays were scanned using an Axon GenePix 4000B or 4200A scanner. Microarray analysis was performed using GenePix What's new? Oropharyngeal squamous cell carcinoma (OPSCC) is often associated with high-risk human papillomavirus infection, especially HPV16, but the clinical relevance of viral integration in a subset of these tumors remains unclear. Here the authors describe transcriptional signatures of metabolic reprogramming including increased expression of the aldo-keto-reductases AKR1C1 and AKR1C3 proteins in OPSCC cells regardless of HPV status. Tumors with upregulation of these proteins were strongly associated with poor prognosis suggesting prognostic as well as therapeutic implications for OPSCC patients beyond HPV16 infection.    Staging was performed according to AJCC/UICC 8th Edition in Oropharyngeal Squamous Cell Carcinoma.
(1): Total number corresponds to the maximal number of patients analyzed. In the HPV-positive group, data from 129 patients were available.
(2): Relative staining compared to normal epithelium. + means higher expression in tumor cells compared to normal epithelium, -means less or equal staining in tumor cells compared to normal epithelium. Pro 6.0.1.25. For normalization processing, the median array intensity was calculated based on the background-subtracted intensity value for all spots excluding control type spots on the array. The background-subtracted intensity value of each spot was then divided by the median array intensity of each microarray. In the 33 tumor samples, 9 samples harbored exclusively integrated HPV16 DNA and the remaining 24 samples harbored episomal viral DNA irrespective of further integrated viral DNA (20 samples episomal viral DNA, 4 samples episomal and integrated viral DNA). Bioinformatic analysis was performed using Partek genomic suite (Partek, St. Louis, Missouri) using unsupervised clustering. Samples with a differential expression (2.0-fold cut-off filter) and p ≤ 0.05 (student's t-test) were filtered for confidence followed by Benjamini-Hochberg false discovery rate (FDR) correction (FDR ≤ 0.1) and selected for further analysis (Supporting  Information Tables 2 and 3). Interaction networks were calculated using Cytoscape including the plugin ClueGo for analysis of Gene-Ontology (GO) and pathways. 14 Primary mRNA array data have been made publicly available at EMBL-EBI (Accession No. E-MTAB-6350) for use in subsequent analysis.

RT-qPCR expression analysis
Five hundred nanograms of total RNA was reverse transcribed using the iScript cDNA Synthesis Kit (BioRad Laboratories Munich, Germany). qPCR reactions were performed using 1 μL of the 20 μL RT product with iTaq Universal SYBR Green Supermix (BioRad). Amplification was performed for 40 cycles with initial denaturation at 95 C for 5 min, followed by cycles with 15 s at 95 C, 30 s at primer specific temperatures for annealing and 40 s at 72 C followed by melting curve analysis. Specific primers are summarized in Supporting Information Table 4. The detection of the housekeeping gene Hypoxanthine Phosphoribosyltransferase (HPRT) was used for normalization of mRNA levels. Human mRNAs from pancreatic islands, complete pancreas, lung, liver and adrenal gland served as controls (Supporting Information Figs. 1J and 2J).

Immunohistochemistry
Immunohistochemical protein staining on 4 μm thick FFPE tissue sections was carried out as described earlier. 15 Briefly, sections were deparaffinized and subsequently pretreated with 3% H 2 O 2 in methanol for 10 min to quench endogenous peroxidase activity. Antigen retrieval was carried out by heating at 70 C in 0.01 M citrate buffer (pH 6.0) over night. AKR1C1 was detected by mouse monoclonal antibodies (Abcam, clone AT6D10, Cambridge, GB; 1:500 in PBS) and AKR1C3 was detected by mouse monoclonal antibody (Sigma Aldrich, clone NP6.G6.A6, Taufkirchen, Germany, 1:500 in PBS). 3-Nitrotyrosine was detected by mouse monoclonal antibody (Santa Cruz Biotechnol., clone 39B6, Heidelberg, Germany, 1:50 in PBS) and NRF2 by a rabbit polyclonal antibody (Sigma Aldrich, HPA002990, 1:200 in PBS). After incubation with a biotinylated secondary antibody, immunohistochemical detection was performed by an avidinbiotinylated peroxidase complex (ABC) procedure (Vectastain-Elite-ABC kit; Vector, Burlingame, USA). Peroxidase activity was detected using DAB substrate kit (Vector) according to the manufacturers protocol. Sections were counterstained with hematoxylin and mounted in Roti Histokit (Carl Roth, Karlsruhe, Germany). Tonsillar epithelium from routine tonsillectomy, pancreas, liver and placenta served as controls (Supporting Information Figs. 1 A-C and 2 A-C). Antibody specificity was validated by control stainings without primary antibody (Supporting Information Figs. 1 and 2 D-F) and IgG isotype control antibodies (data not shown). Analysis was carried out by four independent observers (CUH, LP, JK and SFP), and in case of interobserver variations, consensus was reached by combined examination of the slides. The staining intensity of the tumor was evaluated in relation to the staining intensity in the adjacent squamous epithelium.

Statistics
Epidemiological data combined with HPV-status and AKR1C1 and AKR1C3 immunohistochemical staining were analyzed using cross-tabulations, χ 2 test and Fisher's exact probability test with the software SPSS 21.0. Overall survival and disease-free survival rates were estimated using the Kaplan-Meier algorithm for incomplete observations. The overall survival time was measured from the date of diagnosis to the last date when the patient was known to be alive (censored) or date of death for any reason (uncensored). The disease-free survival rate was defined as the period of time beginning on the date of diagnosis to the day of the last follow-up examination in which the patient was disease-free (censored), or to the date of local recurrence of the disease or occurrence of regional or distant metastases (uncensored). The log-rank test was used to perform the univariate analysis of the various variables. A significance level of p ≤ 0.05 was considered statistically significant in 2-sided tests. Statistical analyses of microarray and RT-qPCR data was performed using GraphPad Prism 6.0 (GraphPad Software, La Jolla California, USA) using student's t-test and ANOVA.

Results
mRNA expression profiling in relation to integration status mRNA expression profiling was performed on 33 HPV16-positive OPSCC with known viral integration status (cohort 1, Supporting Information Table 1). Unsupervised clustering of mRNA expression data separated the samples into two major clusters (Fig. 1a). When we annotated the viral integration data to these clustered expression results, one cluster matched almost entirely with tumors harboring exclusively host genome integrated HPV16 DNA (tumors with integrated virus), and the other one with tumors containing extrachromosomal HPV16 DNA (tumors with episomal virus). Three of the four tumors with mixed forms of integrated and episomal virus also clustered within the group with solely episomal viral DNA.  Tables 2 and 3). By using the Cytoscape software and ClueGo plugin we identified 24 affected cellular pathways with deregulated mRNAs, including epidermal development and differentiation, hormone regulation and processing, oxidative stress response and metabolic processes of ketones, organic acids and (mono) carboxylic acids (Fig. 1b). Interestingly, four human mRNAs were found to be significantly deregulated in multiple affected pathways (Fig. 1b), i.e. AKR1C1 (upregulated 7.5x, p = 0.009 in OPSCC with integrated viral DNA), AKR1C3 (upregulated 7.0x, p = 0.023), BCL2 (downregulated 3.0x, p = 0.019) and BCL2L10 (upregulated 3.4x, p = 0.005). In contrast, the expression of a number of targets known to be involved in HPV-induced carcinogenesis was not significantly correlated with integration status, such as BCLX, CDKN2A, E2F1, EGFR, PIK3CA, RB1, TP53 and TP63.

Confirmation of AKR1C1, AKR1C3, BCL2 and BCL2L10 mRNA expression by RT-qPCR
To confirm the expression level of the four genes identified by GO-term clustering, RT-qPCR was performed on 27 cases out of cohort 1. Similar to mRNA expression profiling, also the RT-qPCR results demonstrated upregulation of AKR1C1, AKR1C3 and BCL2L10 and downregulation of BCL2 in the OPSCC with integrated virus (Supporting Information Fig. 3A-H).

Upregulation of E6*I transcripts in OPSCC with integrated HPV16
Because studies have shown that metabolic pathways including oxidative stress response may be activated by the virus derived splice variant E6*I, 16 we analyzed the mRNA isolated from Cohort I tumors for E6*I expression by RT-qPCR (Supporting Information Fig. 3 I). E6*I expression showed to be significantly upregulated in OPSCC with exclusively integrated viral genomes (p = 0.001). E6*I expression levels in tumors with mixed forms of integrated and episomal virus (n = 3) were similar to those detected in tumors with exclusively episomal virus.

Immunohistochemical detection of AKR1C1 and AKR1C3
To further evaluate the microarray findings, we focused on examining the expression of AKR1C1 and AKR1C3 proteins, because these enzymes have been reported as potential targets for therapy and are upregulated in epithelial tissue of smokers and chemoresistant tumors. 17 First, we analyzed the expression in a number of normal human control tissues including liver, pancreas, placenta and tonsils (see Supporting Information Figs. 1 and 2 also for some corresponding mRNA expression analysis). For both proteins, high expression levels were detected in liver and exocrine pancreas, pancreatic islands showed both cells with high expression and no or very low expression, and no expression was detected in placenta, which is in agreement with stainings reported in human protein atlas. 18 Both in mRNA and protein expression analysis, AKR1C1 appears to be higher expressed than AKR1C3. In normal tonsils (n = 6), lymphocytes were negative and squamous epithelial keratinocytes showed no or weak predominantly nuclear staining, whereas muscle cells and endothelial Cells exhibited strong staining (Supporting Information Figs. 1 and 2 A-C). Control stainings without primary antibodies and IgG isotype controls were negative in all tissues (Supporting Information Figs. 1 and 2 D-F and data not shown).
Subsequently, we analyzed 24 OPSCC of cohort I with sufficient FFPE tissue available, namely 16 cases with episomal virus and 8 with integration and detected a similar variance in protein expression as was found for mRNA by microarray and RT-qPCR analysis when comparing integration status (Supporting Information Fig. 4). Because we noticed some squamous epithelia with positive staining, we also decided to evaluate staining of both tumors and the adjacent epithelia and to relate the staining intensities to each other in further analysis. In these 24 cases, the resulting protein intensity ratios showed a similar pattern as the RT-qPCR results, i.e. lower expression in OPSCC with episomal HPV16 (Supporting Information Fig. 4 B and C).
To evaluate our findings in a separate large cohort (cohort 2), we immunostained 141 OPSCC and their adjacent epithelia including 48 HPV-positive tumors for AKR1C1 and AKR1C3 (Fig. 2). For AKR1C1, 59 out of 141 OPSCC (41.8%) showed stronger staining in the tumor than in the adjacent epithelium (AKR1C1+). The remaining 82 OPSCC (58.2%) showed equal or less staining than in the adjacent epithelium (AKR1C1-). 38% of the HPV-positive tumors were AKR1C1+ and, remarkably also 43% of the HPV-negative tumors were AKR1C1+. Approximately the same percentages were found for AKR1C3 (Tables 1, 2 and Supporting Information  Table 5).

Subgroup analysis of AKR1C stained samples for 3-NT and NRF2
To validate the association between mRNA array data and oxidative stress, we immunostained n = 24 randomly chosen HPV-positive and -negative samples from cohort II against 3-Nitrotyrosine (3-NT), which is generated by Peroxynitrite and the oxidative stress related transcription factor Nuclear Factor (Erythroid-derived 2)-Like 2 (NF2L2/NRF2) 19,20 (Supporting Information Fig. 5). In HPV-positive samples, both 3-NT and NRF2 showed similar staining intensities compared to AKR1C3 staining (p < 0.0001) (Supporting Information Fig. 5 and Supporting Information Table 6). In HPV-negative samples some exceptions were found, i.e. one discordant staining pattern for 3-NT and two for NRF2 in comparison with AKR1C3 (3-NT p = 0.011; NRF2 p = 0.197).

Subgroup analysis of AKR1C stained samples for viral integration
In order to examine if the AKR1C1+/C3+ tumors in the HPV-positive group also show a higher number of cases with viral integration, we studied 20 cases of the HPV-positive OPSCC from cohort II of which sufficient fresh frozen tissue was available for viral integration analysis. APOT-PCR showed 7 cases with integration, 4 with a mixed pattern, and 9 with episomal viral DNA. Five of seven OPSCC with integrated virus were AKR1C1+/C3+, whereas all 9 cases with episomal virus were negative. This is in line with the results from cohort I. The cases with mixed viral physical status showed both positive (n = 3) and negative (n = 1) staining.
Correlation of AKR1C1 and AKR1C3 with clinicopathological data. By correlating the immunohistochemical results with clinicopathological data we demonstrated a highly significant association of AKR1C1+ and AKR1C3+ tumors with tumor relapse and patient death (p < 0.0001) ( Table 2 and Supporting  Information Table 5). Although, there is a strong association for both proteins with HPV-status, this is less evident for history of alcohol (AKR1C1 p = 0.093, AKR1C3 p = 0.044) and tobacco (AKR1C1 p = 0.280, AKR1C3 0.093) consumption (Table 1). Furthermore, HPV-status was associated with known parameters such as a lower T-stage, less tumor relapse and death (Table 2 and Supporting Information Table 5).

AKR1C1/C3 protein expression, HPV and survival in OPSCC
AKR1C1/C3 protein expression did correlate with both worse overall (OS) (Fig. 3) and disease-free survival (DFS) rates (p < 0.0001 in all cases) as did HPV-status (p < 0.0001 in both cases) and low T staging (p = 0.001) and lacking smoking history (p = 0.039) (Supporting Information Table 5 and Supporting Information Fig. 7).
We also correlated separately AKR1C1 and AKR1C3 staining in the tumor and in the adjacent epithelium in relation to survival. This revealed that particularly low expression in the tumor adjacent epithelium is already associated with unfavorable prognosis (AKR1C1 p = 0.0104; AKR1C3 p = 0.0245), and in that combination with higher expression in the adjacent tumor, this association becomes more prominent (p < 0.0001 for both AKR1C1 and AKR1C3) (Fig. 3, Supporting Information Fig. 6, Supporting Information Table 5 and Supporting Information  Table 6).

AKR1C1 / C3 protein expression and treatment modalities
Since the AKR1C enzymes are reported to be Phase I detoxifying enzymes especially involved in the turnover of platin drugs, we correlated AKR1C staining intensities with treatment modalities (Supporting Information Table 5 and Supporting Information Fig. 7). Consequently, AKR1C1+ / C3+ turned out to be a negative predictor in correlation to Chemo-and Radiotherapy both in OS and DSF. Most prominent, low AKR1C1 and AKR1C3 expression turned out to be a highly significant predictors for favorable outcome in surgical  (Table 2 and Supporting Information Fig. 7) for OS and DFS. T-stage, HPV and AKR1C3 turned out to be independent significant prognostic predictors.

Discussion
It remains unclear if HPV integration in OPSCC affects the levels of viral and/or HPV-disrupted human gene transcripts. In previous work, we did not find clear evidence for this hypothesis, 3,21 since constitutive, rather than a high level of viral and virally disrupted gene transcripts was observed in the tumors. Therefore, in our study it was our aim to perform a genome-wide screen to identify human genes which expression is influenced by integration of HPV16. By comparing HPV16-positive OPSCC harboring episomal or integrated virus using mRNA microarray expression profiling, we identified a unique signature of differentially expressed human mRNAs in relation to viral physical state. Upregulated AKR1C1 and AKR1C3 expression was observed in HPV16-positive OPSCC with viral integration and, interestingly, associated with poor prognosis both in HPV-positive and -negative tumors.
Despite some literature reporting that in OPSCC with integrated HPV there may be a direct effect of viral protein expression on the host genome and human gene expression, 9,22 in most tumors this effect is unclear. Here we provide new data showing expression signatures associated with HPV16 physical status in OPSCC. Interestingly, the tumors with viral integration showed deregulated expression of genes involved in metabolic pathways (Fig. 1b), frequently including upregulated AKR1C1 and/or AKR1C3 expression. Furthermore, this upregulated expression was also observed in a subgroup of HPVnegative OPSCC (see further below). Increased expression of AKR1C genes has previously been reported by others. Martinez et al. compared gene expression in a small group of HPVpositive vs. -negative OPSCC by mRNA expression profiling, finding higher expression of AKR1C3 in HPV-negative tumors. 23 In a protein expression study on uterine cervical cancer (n = 145), Ueda et al. found a significant correlation between AKR1C1 (termed by alternative nomenclature as DDH1, Dihydrodiol Dehydrogenase Type I) and tumors infected with HPV types 16 and 18, FIGO stage, lymph node involvement and patients' survival. 24 AKR1C1 and AKR1C3 upregulation has also been observed after transfection of the HPV-negative uterine cervical carcinoma cell line C33A with a HPV16-E6*I construct. 25 The translated truncated HPV16-E6 product E6*I appeared to have a direct promoting effect on AKR1C1 and AKR1C3 expression via the SP1 binding sites in their promotor regions.
Our observation that the AKR1C1 and AKR1C3 genes are upregulated in OPSCC with integrated HPV16 do fit with these experiments, since we and others also detected E6*I mRNA upregulation in tumors with integrated HPV16. 26,27 Like tobacco smoke, E6*I can promote reactive oxygen species (ROS) activated pathways, including activation of NRF2. 16 Known consequences of NRF2 activation are e.g. activated PI3K-AKT signaling, deficiency in autophagy, impaired DNA damage response, enhanced cell proliferation, and interestingly, metabolic reprogramming and chemotherapeutic drug modification resulting in resistance 17,28,29 (Fig. 4).
AKR1C1 and AKR1C3 belong to the aldo-keto reductase (AKR) superfamily of NAD(P)H dependent oxidoreductases. Both enzymes have multiple cellular functions in e.g. regulating prostaglandin, steroid hormone and retinoid metabolism and moreover are Phase I detoxifying enzymes involved in modifying toxic substances, such as chemotherapeutic drugs and tobacco smoke components. 17,30,31 Besides direct upregulation by E6*I, the expression of both enzymes can be initiated under the presence of ROS by the Keap1/Cul3/Nrf2 system via direct binding of NRF2 to the antioxidant response element (ARE) located in the promotor regions of the AKR1Cs and in several additional oxidative stress related genes which proved to be deregulated in our study (Fig. 4). Interestingly, AKR1C1 and AKR1C3 generate several feedback loops controlling their own expression. On the one hand, through the production of retinol from retinaldehyde, so that the concentration of retinoic acid, a known inhibitor of NRF2 function, is lowered, which leads to additional NRF2 activation. 32 On the other hand, in response to tobacco smoke, they exacerbate the carcinogenicity of polycyclic aromatic hydrocarbons (PAH) by oxidizing trans-hydrodiol functional groups (leading to e.g. benzo[a]pyrene-7,8-dione [BPQ]), which in turn binds AREs and further stimulates AKR1C1 and AKR1C3 expression. In this way PAH transhydrodiols would enhance their own genotoxicity by inducing expression of the AKR1Cs. 17 Because tobacco smoke thus can result in upregulation of AKR1C1 and AKR1C3 and 86% of patients in our cohort of 141 OPSCC were smokers (cohort II), we assessed if there was an association between these parameters by analyzing expression in the tumors with their adjacent epithelium. Interestingly, only in the adjacent epithelium there was prognostic difference between high and low stained samples (Supporting Information Fig. 6). This prognostic effect even became higher when the normal epithelium staining intensity was at least one grade higher than the tumor staining intensity. From these data and expression of NRF2 and 3-NT (Supporting Information Fig. 5 and Supporting Information Table 6, Fig. 4) we conclude that high staining in the epithelium is an indicator for a better capacity to handle oxidative stress. As a consequence, higher staining in the tumor than in the epithelium is overwhelming this individual capacity to deal with oxidative stress resulting in a more aggressive tumor. Our approach of semi-quantitative normalization of tumor staining intensity by comparing it with adjacent epithelium might be furthermore supported by a study from Namani et al., where RNAseq data from the TCGA HNSCC cohort were normalized against normal tissue expression and subsequently lead to a signature of overexpressed oxidative stress response genes including AKR1C1 and AKR1C3. 33 For confirmation, we analyzed Host genome integration of HPV16 leads to upregulation of the viral E6 spliced isoform E6*I. Like tobacco smoke, E6*I can promote reactive oxygen species (ROS) signaling by influencing transcriptional factor expression like NRF2 (Nuclear Factor Like 2) to modulate the expression of antioxidant enzymes including AKR1C1 and AKR1C3. 25 The AKR1Cs 1 and 3 are known to catalyze the production of retinol isoforms from retinaldehyde and subsequently lowering retinoic acid concentrations. 37 Retinoic acid is a known inhibitor of NRF2 function, therefore decreased concentrations of retinoic acid lead to additional NRF2 activation. 32 Furthermore, NRF2 controlled NQO1 (NAD(P)H dehydrogenase [quinone] 1) stabilizes p53, preventing it from degradation. 47 Several additional NRF2 activated oxidative stress response genes harboring ARE elements were found to be upregulated by our mRNA array analysis including CD36 (cluster of differentiation 36, also known as fatty acid translocase), DUOX1 (Dual Oxidase 1), GPX3 (Glutathione Peroxidase 3), MPO (Myeloperoxidase), SOD2 (Superoxide Dismutase 2), and SRXN1 (Sulfiredoxin 1). 17,48,49 In HPV16 integration positive OPSCC, HPV16-E6 might furthermore be blocked by its isoform E6*I. 5 As a consequence, we found BCL2 (B-cell lymphoma 2) to be downregulated. Through increased expression of its antagonist BCL2L10, putative positive effects on apoptosis and autophagy might be directly antagonized. Known consequences of NRF2 activation are e.g. activated PIK3-AKT signaling, metabolic reprogramming, enhanced cell proliferation, deficiency in autophagy, and interestingly, resistance to chemotherapy as well as impaired DNA damage response. The latter is supposed to promote integration of HPV16 episomes into the human genome.Green = upregulation, red = downregulation in mRNA array analysis.
[Color figure can be viewed at wileyonlinelibrary.com] RNAseq data from the TCGA HNSCC cohort with both tumor and corresponding normal tissue expression data available (n = 43). Again, overexpression of AKR1C1 and AKR1C3 in the tumor compared to corresponding normal tissue correlated with unfavorable survival (AKR1C1 p = 0.0406; AKR1C3 p = 0.0232) (data not shown).
The presence of a differential expression of AKR1C1 and AKR1C3 between tumor and epithelium might be an optimal readout for this metabolic reprogramming pattern. 17,34 We found no correlation between smoking and upregulation of AKR1C1 and AKR1C3 expression. An additional mechanism promoting NRF2 and subsequently AKR1C1 and AKR1C3 activation is through methylation and mutation of genes coding for proteins of the oxidative stress system (e.g. CUL3, KEAP1, NRF2, RBX1), especially found in HPV-negative HNSCC. 34,35 Consequently, in HPV-negative samples we found 18.2% cases with discordant NRF2 staining intensities compared to AKR1C staining (Supporting Information Table 6). Interestingly, in an independent study comparing methylation profiles of HPVpositive and -negative OPSCC, methylation of the ALDH1A2 gene and subsequent downregulation of gene expression also was related to retinol production, activation of oxidative stress response and an unfavorable prognosis 36,37 (Fig. 4). This is in agreement with our findings presented here.
By univariate analysis we found several parameters to be associated with unfavorable prognosis in OPSCC, including upregulation of AKR1C1 and AKR1C3 along with known factors such as HPV, T-stage, N-stage and smoking. HPV, T-stage, AKR1C1 and AKR1C3 turned out to be significant also in multivariate analysis. Moreover, AKR1C1 and AKR1C3 expression turned out to be strong indicators of prognosis in both HPVpositive and -negative OPSCC. Interestingly, smoking proved not to be an independent variable in multivariate analysis, which might be explained by the fact that AKR1C1 and AKR1C3 are indicators of oxidative stress including, but not restricted to tobacco smoke components. 31 If these results could be confirmed in larger studies, than AKR1C1 and/or AKR1C3 may be considered to be included in prediction models in OPSCC and HNSCC. [38][39][40] Low risk groups can then profit from de-intensification of treatment protocols, 41 whereas intermediate and high risk groups could be selected for other therapeutic options, such as inhibitors of the PI3K and NRF2 pathways including the AKR1C proteins. 17,42,43 Potential treatment options addressing oxidative stress related pathways involving NRF2, AKR1C and PI3K upregulation were recently discussed. [44][45][46] However, due to the high expression of AKR1C enzymes in normal tissues like liver and pancreas, it will be important to investigate potential adverse events of these drugs in these tissues.
In conclusion, we provide further evidence that viral host genome integration of HPV16 has a biological effect on the host cell, e. g. by deregulating metabolic pathways including upregulating AKR1C1 and AKR1C3 expression. Increased expression of these proteins was also found to have impact on prognosis in HPV-negative tumors. Therefore, our findings may impact prognosis as well as increase options for treatment of OPSCC in general.