# Comparative RNA sequencing reveals that HPV16 E6 abrogates the effect of E6*I on ROS metabolism

Philippe Paget-Bailly, Koceila Meznad, Diane Bruyère, Jérôme Perrard, Michaël Herfs, Alain Jung, Christiane Mougin, Jean-Luc Pretet, Aurélie Baguet

<figure><img src="image_1.png" /></figure>

To cite this version:

Philippe Paget-Bailly, Koceila Meznad, Diane Bruyère, Jérôme Perrard, Michaël Herfs, et al. Comparative RNA sequencing reveals that HPV16 E6 abrogates the effect of E6*I on ROS metabolism. Scientific Reports, 2019, 9 (1), 10.1038/s41598-019-42393-6. hal-03169777

HAL Id: hal-03169777

https://univ-fcomte.hal.science/hal-03169777

Submitted on 3 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_3.png" /></figure>

OPEN

# Comparative RNA sequencing reveals that HPV16 E6 abrogates the effect of E6*I on ROS metabolism

Received: 29 October 2018

Accepted: 27 March 2019

Published online: 11 April 2019

Philippe Paget-Bailly<fnref n="1" />,², Koceila Meznad<fnref n="1" />,², Diane Bruyère⁴, Jérôme Perrard<fnref n="1" />,²,

Michael Herfs⁴, Alain C. Jung⁵, Christiane Mougin<fnref n="1" />,²,³, Jean-Luc Prétet<fnref n="1" />,²,³ & Aurélie Baguet<fnref n="1" />,²

High-risk Human Papillomavirus infections are responsible for anogenital and oropharyngeal cancers. Alternative splicing is an important mechanism controlling HPV16 gene expression. Modulation in the splice pattern leads to polycistronic HPV16 early transcripts encoding a full length E6 oncoprotein or truncated E6 proteins, commonly named E6*. Spliced E6*I transcripts are the most abundant RNAs produced in HPV-related cancers. To date, the biological function of the E6*I isoform remains controversial. In this study, we identified, by RNA sequencing, cellular targets deregulated by E6*I, among which genes related to ROS metabolism. Concomitantly, E6*I-overexpressing cells display high levels of ROS. However, co-overexpression of both E6 and E6*I has no effect on ROS production. In HPV16-infected cells expressing different E6/E6*I levels, we show that the newly identified targets CCL2 and RAC2 are increased by E6*I but decreased by E6 expression, suggesting that E6 abrogates the effect of E6*I. Taken together, these data support the idea that E6*I acts independently of E6 to increase ROS production and that E6 has the ability to counteract the effects of E6*I. This asks the question of how E6*I can be considered separately of E6 in the natural history of HPV16 infection.

Human papillomaviruses (HPVs) are small non-enveloped viruses that present a tropism for squamous epithelium. More than 200 types of HPVs have been described to infect humans<fnref n="1" />. Based on their oncogenic potential, these viruses are classified in high-risk HPV (hrHPV), including HPV16 and HPV18, and low-risk HPV (IrHPV), including HPV6 and HPV11. HPV infections are responsible for cervical intraepithelial lesions that can progress to cancers, but they also cause a large fraction of anal, vulvar, vaginal, penile cancers, and a rising number of oropharyngeal cancers.

HPV16 is the most prevalent type in HPV-associated cancers². Its genome contains a long control region (LCR), 6 open reading frames (ORFs) encoding early (E) proteins under the control of p97 promoter located in the LCR, and 2 ORFs encoding late (L) proteins under the control of p670 promoter located within the E7 ORF³. Viral proteins are produced through the translation of at least 20 polycistronic transcripts obtained by alternative splicing. At least 10 of these transcripts allow the production of the 2 major viral oncoproteins, E6 and E7, but also E6-truncated proteins, E6*I, E6*II and E6^E7⁴,⁵. hrHPV E6 and E7 proteins are consistently expressed in HPV-associated cancers⁶,⁷ and interact with many host cellular proteins. Notably, E6 and E7 proteins target p53 and pRB, respectively, for proteasome-mediated degradation, and thus inactivate these tumor suppressors⁶,⁸,⁹. More than 30 years ago, it was observed that the most abundant HPV16 transcript produced was spliced from the donor site 226 (SD 226) to the acceptor site 409 (SA409), both sites located in E6 ORF¹⁰,¹¹. Interestingly, only hrHPVs harbor these splice sites, indicating that E6 ORF splicing events could be relevant for HPV-driven carcinogenesis¹⁰. Several studies and unpublished data from our laboratory also reported increased levels of spliced E6*I mRNA correlating with the severity of cervical lesions⁵,¹²⁻¹⁴. Concomitantly, it was proposed that the ratio

<footnote n="1">¹EA3181, LabEx LipSTIC ANR-11-LABX-0021, UFR Santé, 19 rue Ambroise Paré, Besançon, France. ²Université Bourgogne Franche Comté, Besançon, France. ³Centre Hospitalier Régional Universitaire, CNR HPV, 3 Bvd Alexandre Fleming, Besançon, France. ⁴Laboratory of Experimental Pathology, GIGA-Cancer, University of Liege, Liege, Belgium. ⁵Université de Strasbourg, Inserm, UMR_S1113, Centre de lutte contre le cancer Paul STRAUSS, Strasbourg, France. Correspondence and requests for materials should be addressed to A.B. (email: aurelie.baguet@univ-fcomte.fr)</footnote>
of E6*I/E6*II spliced variants can be used as a predictive marker of clinical outcome in HPV-related cervical lesions¹³ and oropharyngeal cancers¹⁵.

Even if HPV16 early transcripts detection is used as a tool in screening and investigating HPV-related neoplasia, the biological significance of E6 splicing and ensuing E6*I protein expression remains elusive. It has been proposed that E6 ORF splicing facilitates translation re-initiation of the E7 ORF by increasing the space between E6 and E7 ORF¹⁶,¹⁷. However, other study showed that E7 is preferentially translated from the unspliced E6/E7 transcript rather than from the E6*I/E7 one, suggesting that this splicing event regulates E6 expression but not E7¹⁸,¹⁹. Apart from HPV gene expression regulation, the roles of E6*I ORF product in HPV life cycle and carcinogenesis also remains unclear, although a variety of functions have been reported for the truncated isoform. E6*I inhibits E6-mediated degradation of p53²⁰, causes the degradation of some PDZ proteins, such as Dlg, PATJ and MAGI-1²¹,²², and modulates the expression of a subset of cellular factors involved in stress response, such as aldo-keto reductase genes²³, superoxide dismutase isoform 2 (SOD2), and glutathione peroxidase 1 (GPX1), leading to the accumulation of reactive oxygen species (ROS)²⁴. *In vivo* studies have shown that ectopic E6*I reduces tumor growth of both HPV-positive and HPV-negative cells in a xenograft nude mouse model²⁵. Depending on the cellular context and especially the presence or absence of the E6 protein, E6*I seems implicated in different cellular pathways and its functions remain controversial²⁶.

To get a better understanding on the underlying molecular mechanisms driving HPV-related carcinogenesis, the present study used RNA sequencing technology to analyze the impact of HPV16 E6*I isoform on cellular gene expression.

## Results

**HPV16 E6 and E6*I expression in U-2 OS HPV-negative cell line.** HPV-negative U-2 OS cell line was transfected with expression vectors encoding either all HPV16 E6 isoforms (pXJ40-E6All or pXJ40b-GloΔint-E6All) or E6 (pXJ40-E6 or pXJ40bGloΔint-E6) or E6*I (pXJ40-E6*I) exclusively or with empty vectors (pXJ40 or pXJ40-bGloΔint) (Figs S1 and S2). We compared alternative splicing profiles of E6 transcripts obtained from these different E6-expression vectors. An accumulation of E6 unspliced transcripts produced from the pXJ40bGloΔint-E6All vectors is observed, indicating that the presence of the β-globin intron in the pXJ40 vector favor the splicing between SD226 and SA409 splice sites, and consequently the production of the E6*I RNAs (Fig. S2). Western blot analysis were performed with two E6-specific antibodies. The 6F4 antibody recognizes the N-terminal part of the E6 protein, while the 3F8 antibody targets the C-terminal part of E6 (Fig. 1A). Consequently, the 3F8 antibody does not recognize the E6*I protein and the 6F4 antibody recognizes both E6 isoforms. The full length E6 protein was detected with both E6-directed antibodies in cells transiently expressing the pXJ40-E6All and pXJ40-E6 vectors. No E6*I protein was observed in pXJ40-E6All or pXJ40-E6*I expressing cells. But consistently with its role as a potent E6 antagonist, we observed an increase of p53 protein level in pXJ40-E6*I transfected cells. (Figs 1B and S3A).

MG-132 treatment, a proteasome inhibitor, was used on transfected cells to reveal a specific band with the 6F4 antibody that corresponds to the E6*I isoform (Fig. 1C, lane 8 and Fig. S3B). This indicate that an E6*I protein can be produced from the pXJ40-E6*I vector but also suggest that low levels of E6*I protein can be explain, at least in part, by its rapid turnover in these cells.

**Establishment of stable cell lines expressing HPV16 E6 isoforms.** Stable cell lines expressing either all E6 isoforms or E6*I were generated from the U-2 OS cell line. Cellular clones were isolated upon selection of cells resistant to G418: U-CTL (clones #1 and #3, transfected with pXJ40-empty), U-E6All (clones #1, #6 and #7, transfected with pXJ40-E6All), and U-E6*I (clones #1, #6 and #14, transfected with pXJ40-E6*I). Then, these cellular clones were tested for E6 and E6*I expression. As expected, E6 protein was observed in U-E6All stable cell lines compared to U-CTL and U-E6*I. A decrease of the wild-type p53 protein level is observed in these cells, suggesting that U-E6All cells express of a functional E6 protein (Fig. 2A and S4A). Next, we compared E6 expression from U-E6All clones to naturally-infected HPV16 cell lines such as SiHa, Ca Ski and W12 clones. We observe a significantly higher expression of E6 in U-E6All clones. Interestingly, this high expression does not result in a greater degradation of p53 compared to SiHa and Ca Ski cells. We observe low amount of E6 in the three W12 clones tested but still managing variable levels of p53 between the three clones tested (Figs 2B and S4B). Because E6*I protein was hardly detectable, E6 and E6*I RNA levels were simultaneously monitored by RT-PCR. The RNA profile of U-E6All stable cell lines shows a low amount of full-length E6 RNA (unspliced) and E6*II RNA (SD226/SA556 spliced RNA) in comparison with the high amount of E6*I RNA (SD226/SA409 spliced RNA). As expected, U-E6*I stable cell lines express only E6*I RNA and U-CTL stable cell lines, like parental U-2 OS (U-CTL#par), did not express E6 or E6*I (Figs 2C and S5A). Similar splicing profiles of E6 transcripts are detected between U-E6All stable cell lines, HPV16-positive cancer cell lines (SCC90, SiHa, Ca Ski) and in most head and neck squamous cell carcinoma (HNSCC) tumor samples tested (Figs 2D and S5B). Also, in the naturally HPV16-infected cervical keratinocyte cells (W12 cells), E6*I transcript is the most abundant splice product compared to E6 and E6*II transcripts (Figs 2E and S5C).

**E6 and E6*I expression differentially alters phenotypic characteristics relevant to cancer cells.** To evaluate phenotypic characteristics relevant to cancer cells in our cellular models, functional assays were performed with U-CTL, U-E6All and U-E6*I cell lines. First, cell proliferation was assessed by measuring the number of cells once a day over 4 days. Results show that all cell lines had the same growth rate (Fig. 3A). In contrast, colony formation assays, performed to monitor the ability of a single cell to grow into a colony, show a significant decrease in the number of colonies for U-E6All compared with U-CTL. We observe that expression of E6*I alone has a stronger inhibitory effect on colony formation (Fig. 3B,C). E6 and E6*I from HPV16 have both been shown to induce oxidative stress in different cellular models²⁴,²⁷. Relative production and release of H₂O₂ in medium
<figure><img src="image_4.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_6.png" /></figure>

**Figure 1.** Transient transfection of E6 expression vectors in U-2 OS HPV-negative cell line. (A) Schematic representation of the E6 and E6*I protein isoforms with indicated regions recognized by the two E6-specific antibodies. The 6F4 antibody can recognize the N-terminal part of the E6 and E6*I proteins and the 3F8 antibody can target the C-terminal part of E6 protein. (B) Western blotting analysis showing p53, E6 and E6*I protein levels in U-2 OS cells transiently transfected with the pXJ40 empty vector (CTL), pXJ40-E6All vector (E6All), pXJ40-E6 vector (E6) or pXJ40-E6*I vector (E6*I). After 48-hour incubation, cells were treated either with DMSO or 10µM MG132 proteasome inhibitor for 16 hours. The membrane was then incubated with the two E6 antibodies. For the detection of E6*I (corresponding to the first 43 amino acids of E6 protein), only the 6F4 antibody can be used. "nd" referred to a band at 11 kDa corresponding to a cleavage product observed in U-E6All cells treated with MG132. Full-length blots are presented in Figure S3.

from cells, reported by the oxidation of amplex red reagent to fluorescent resorufin, was used as indicator of ROS production. We observe a 20% significant increase in H₂O₂ released from U-E6*I clones relative to control. However, no change in H₂O₂ release from U-E6All clones is observed (Fig. 3D). It was previously observed that HPV16 E6 increases the radio sensitivity of HPV-negative cancer cells²⁸. To confirm this in our cellular models, the radiation sensitivity was measured *via* clonogenicity assays following irradiation. In all cell lines, survival decreases with increasing irradiation dose (Fig. 3E). We observed that only U-E6All cell lines show a significant decrease of survival after irradiation compared to U-CTL and U-E6*I cell lines (Fig. 3F). These data suggest that E6 and E6*I can have isoform-specific functions in the cells.

**Transcriptome profiles reveal that E6*I deregulates cellular targets independently of E6.** To assess cellular targets potentially deregulated by E6*I, the whole transcriptome of stable cell lines – U-E6All (clones #1, #6 and #7), U-E6*I (clones #1, #6 and #14) and U-CTL (clones #1, #3 and #par) – were analyzed by RNA sequencing. First, unsupervised hierarchical clustering analysis based on RNA sequencing datasets reveals 2 clusters. One of these clusters contains the 3 clones of U-E6All cells (Fig. 4A and S6). Venn diagram representation of deregulated genes in U-E6All and U-E6*I compared with U-CTL show that 419 transcripts were significantly deregulated (FDR < 0.05) in U-E6All (Fig. 4B), of which 264 were upregulated and 155 were down-regulated (Fig. 4C and Table S1). The expression of 53 transcripts were deregulated in U-E6*I cell lines (Fig. 4B), of which 45 were up-regulated and 8 were down-regulated (Fig. 4C and Table S2). An overlap of 11 deregulated transcripts was observed between U-E6All and U-E6*I cell lines. These data suggest that E6*I modulates cellular gene expression independently of E6 but also indicate that the effect of E6*I on transcriptome is disrupted when co-expressed with E6.

**GO terms and KEGG pathway analysis of deregulated transcripts in U-E6All and U-E6*I cells.** Classification using GO (Gene Ontology) enrichment and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis showed an over-representation of U-E6All deregulated transcripts in enriched terms,
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

**Figure 2.** Expression profiles of HPV16 E6 and E6*I in stable cell lines. (A) Western blotting analysis showing p53, E6 and β-actin protein levels in U-2 OS parental cell line (U-CTL#par) and in stably transfected cells with pXJ40 empty vector (U-CTL #1 and #3) or pXJ40-E6All vector (U-E6All #1, #6 and #7) or pXJ40-E6*I vector (U-E6*I #1, #6 and #14). β-actin is used as loading control. Full-length blots are presented in Figure S4. (B) Western blotting analysis comparing protein levels of E6 and p53 between U-2 OS cellular clones (U-CTL#par, U-E6All #1, #6 and #7) and naturally HPV16-infected cells (W12_20861, W12_20862, W12_20863, SiHa and Ca Ski cell lines). β-actin is used as loading control. (C) RT-PCR analysis of alternatively spliced HPV16 E6 transcripts in U-E6All and U-E6*I stably transfected cells. As expected no E6 transcripts were observed in U-2 OS parental cell line (U-CTL#par) and in U-CTL #1 and #3 stably transfected cells. The detection of the vector in stably transfected cells was done by pXJ40 vector specific primers and RPLP0 was used as loading control. Full-length gels are presented in Figure S5. (D) RT-PCR analysis of alternatively spliced E6 transcripts in HPV16 + cell lines SCC090, SiHa, Ca Ski, in stable cell lines U-E6All #1, #6 and #7 and in HPV16 + HNSCC tumor samples. RPLP0 was used as loading control. Note that a band can be observed at 400 bp and could correspond to heteroduplexes between the different amplicons generated during PCR reaction. Full-length gels are presented in Figure S6. (E) RT-PCR showing alternatively spliced E6 transcripts in naturally HPV16-infected cells (W12_20861, W12_20862, W12_20863).

such as single-organism cellular process (GO: 0044763), protein-metabolic process (GO: 0019538), metabolic pathways (KEGG: hsa01100), ribosome (KEGG: hsa03010) and pathways in cancer (KEGG: hsa05200) (Fig. 5A). Based on protein-coding transcripts deregulated in U-E6All cells, a protein interaction network was constructed using the STRING database (Fig. 5B). Many of the deregulated transcripts identified in U-E6All cells have been previously identified in HPV-infected cells and tumor samples such as transcripts encoding ribosomal subunits, SYCP2, NUP210 and DICER1<sup>29-33</sup>.

Gene ontology and KEGG pathway analysis of U-E6*I cell lines revealed 6 gene classes deregulated by ectopic expression of E6*I: Extracellular Matrix Organization (GO:0030198), Response to decreased oxygen levels (GO: 0036293), MAPK signaling pathway (KEGG: hsa04010), Focal adhesion (KEGG: hsa04510), Cytokine-cytokine receptor interaction (KEGG: hsa04060) and Pathway in cancer (KEGG: hsa05200) (Fig. 6A). In addition, the STRING database reported interactions between 10 protein-coding transcripts based on text mining, protein-protein interaction and co-expression data. Furthermore, these data show that the core interactors belong
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 3.** HPV16 E6 and E6*I expression differentially alters phenotypic characteristics relevant to cancer cells. (A) Cell proliferation assay with stably transfected cells U-CTL (black), U-E6All (red) and U-E6*I (blue) was performed to monitor daily live cells for 4 days. Cell counts are represented as mean +/− S.D. of 3 independent experiments. (B) Representative images of the colony formation assay realized with U-CTL#par, U-CTL#3, U-E6All#7 and U-E6*I#1. (C) Quantification of the colony formation assay realized with U-CTL (#par, #1, #3), U-E6All (#1, #6 and #7) and U-E6*I (#1, #6 and #14) cell lines. Foci counts are represented as mean +/− S.D. of 3 independent experiments. (D) Detection of H₂O₂ released from U-CTL (#par, #1, #3), U-E6All (#1, #6 and #7) and U-E6*I (#1, #6 and #14) cell lines in 1 h using the Amplex Red assay. H₂O₂ relative levels are represented as mean +/− S.D. of 3 independent experiments. (E/F) Clonogenic radio sensitivity assay performed with U-CTL, U-E6All and U-E6*I cell lines. Cell colonies (>50 cells) were counted 10 days following irradiation (0 to 8 Gy). Mann-Whitney test were used to assess statistical significance (*p* < 0.05; **p* < 0.01; ***p* < 0.001).</figcaption></figure>

to the extracellular matrix organization and/or oxidative stress response pathways (Fig. 6B). These data reveal specific pathways deregulated by E6*I in an E6-independent manner.

**RNA sequencing transcripts validation by RT-qPCR.** To validate RNA-seq data, expression of several transcripts was confirmed by RT-qPCR. We selected SYCP2, NUP210 and DICER1 transcripts deregulated in U-E6All cells and already known to be modulated by E6 expression (Fig. 7A), the 11 transcripts deregulated in both U-E6All and U-E6*I cells (Fig. 7B) and 13 transcripts only deregulated in U-E6*I cells including the 6 transcripts assigned to “response to decreased oxygen levels” (GO: 0036293) (Fig. 7C). Except for the lincRNA (long intergenic non-coding RNA) named CH507-210p18.1, fold changes estimated by RT-qPCR and RNA-sequencing are consistent for all transcripts tested. Pearson correlation coefficients (R) were calculated. We found high
<figure><img src="image_18.png" /><figcaption>**Figure 4.** Transcriptome profiles reveal that E6*I deregulates cellular targets independently of E6. (A) Unsupervised hierarchical clustering of RNA-seq data from the U-CTL, U-E6All and U-E6*I cell lines. (B) Venn diagram showing the number of deregulated genes in the 2 pairwise comparisons U-E6All vs. U-CTL and U-E6*I vs. U-CTL. (C) Number of deregulated transcripts in U-E6All and U-E6*I cell lines compared to U-CTL.</figcaption></figure>

correlations (R = 0.996 for U-E6All-deregulated transcripts, R = 0.853 for both U-E6All- and U-E6*I-deregulated transcripts and R = 0.993 for U-E6*I-deregulated transcripts) indicating a strong concordance between RNA-seq and RT-qPCR data.

**Expression pattern of E6 or E6*I deregulated transcripts in W12 cell line.** To validate the biological significance of the cellular transcripts found deregulated by E6*I, we performed experiments on W12 cells, a naturally HPV16-infected cervical keratinocyte cell line. First, absolute quantification was done to measure the copy number of E6 and E6*I mRNAs in three undifferentiated W12 clones (W12_20861, W12_20862, W12_20863). The W12_20861 clone expresses the highest level of E6*I transcript (approximately 4.10⁶ copies/µL) compared to W12_20862 and W12_20863 clones (approximately 500 000 copies/µL) (Fig. 8A). Unspliced E6 transcripts are highly expressed in the W12_20863 clone (approximately 1.10⁶ copies/µL) compared to W12_20861 and W12_20862 clones (approximately 230 000 and 14 000 copies/µL, respectively). In each W12 clone tested, we observed a strong correlation between RNA and protein levels for E6 with high level of E6 protein and low level of p53 protein in W12_20863 clone compared to W12_20861 and W12_20862. The high levels of E6*I mRNA and E7 protein in W12_20861 clone is also consistent with the lowest level of pRB protein (Figs 8B and S7). Overall, each clone possesses a unique expression profile of HPV transcripts and proteins allowing us to correlate those profiles to the expression of cellular targets identify by RNA-seq. So, the expression levels of transcripts modulated by E6 such as SYCP2, p21 and TNFAIP6 was measured by RT-qPCR in each W12 clone. As expected, SYCP2 and TNFAIP6 transcripts are upregulated while p21, a p53 target gene, is downregulated in the W12_20863 clone expressing the highest level of E6 (Fig. 8C).

Next, we monitored the RNA levels of CCL2, RAC2 and PDGFB transcripts modulated by E6*I. These transcripts encode proteins implicated in the activation of the NADPH Oxidase complex-2 (Nox2), required for Nox-dependent ROS production. There were higher levels of CCL2 and RAC2 in W12 clone expressing high level of E6*I (W12_20861). Nevertheless, a decreased expression of these transcripts is observed in cells expressing high level of E6 (W12_20863) suggesting that E6 could have the ability to counteract the effects of E6*I. Altogether, these results show that E6*I expression could drive the expression of genes implicated in oxidative stress metabolism in the W12 cell lines but also indicate that E6 can abrogate the effect of E6*I on ROS metabolism.

Discussion

To our knowledge, this study is the first report assessing the impact of HPV16 E6*I on cellular gene expression using RNA sequencing. The aim of this study was to uncover target genes specifically modulated by E6*I independently of the full length E6 protein. To this end, stable cell lines expressing either all E6 protein isoforms (U-E6All) or E6*I (U-E6*I) exclusively, were generated from U-2 OS. The U-2 OS cell line represents an adequate model for studying HPV gene expression and HPV-associated carcinogenesis because it expresses a wild-type p53, which can be targeted by E6. Furthermore, it was shown that viral gene expression of both alpha and beta papillomaviruses in U-2 OS is very similar to that described in keratinocytes³⁴. Immunoblotting experiments allowed to confirm the expression of a functional E6 by loss of p53 in U-E6All cell lines³⁵. E6 early transcripts were also efficiently spliced in U-E6All cell lines, yielding high amount of E6*I transcripts as observed in HPV16-positive cell lines and HNSCC clinical samples³⁶. However, although the E6*I protein was detected in transient transfection, no protein was observed in stable cell lines, probably due to the extremely short half-life of the E6*I isoform. This was already documented in a previous study, which reported the detection of an HPV16 E6*I fusion protein following transient transfection but failed to detect E6*I in stable cell lines³⁷. In the present study, we observe an accumulation of endogenous p53 protein in the HPV16 E6*I-expressing U-2 OS cells, suggesting that E6*I could have opposite effect of E6 on p53 protein level. But, we also show that E6*I expression causes an accumulation of ROS, indicating that the increasing level of p53 observed can be mediated by the ROS production. Also, emerging roles of non-coding RNAs in gene expression regulation and the noticeable difficulties to detect the E6*I protein by western blotting encouraged us to consider that the transcriptome deregulation
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 6.** Enrichment and pathway analysis of deregulated genes in U-E6*I cell lines compared with U-CTL. (A) Gene Ontology and KEGG pathway analysis of differentially expressed transcripts in U-E6*I cell lines. (B) Protein-protein interaction network, corresponding to the protein-encoding genes deregulated in U-E6*I cell lines, generated using STRING database⁷⁵. Line thickness indicates the strength of data support. In blue and red are represented proteins involved in extracellular matrix organization and response to decreased oxygen levels, respectively.</figcaption></figure>

expressing E6 isoforms or only E6*I²⁴, we report that stable expression of E6*I increases ROS production but no change is observed in cells co-expressing E6. Finally, radio sensitivity assays showed that U-E6All cells are more sensitive to radiation than both E6*I and control cells. These latter results could be partially explained by the TP53 wild-type status of our cellular model given the fact that E6 was shown to disrupt the p53-mediated cellular responses to DNA damage⁴¹. Altogether, these data suggest that E6 isoforms have a specific impact on cellular phenotype and that the impact of E6*I is modulated when co-expressed with E6.

Recently, a whole transcriptome study using differentiated keratinocytes showed that HPV16 infection modulates genes involved in the regulation of cellular matrix adhesion, ROS biosynthesis and inflammatory response⁴². In our study, RNA sequencing identified 419 transcripts significantly deregulated in HPV-negative cellular model expressing E6 isoforms. Among them, a substantial number was already described in the literature as deregulated in HPV16 cellular models or in HPV16-positive clinical samples such as the down-regulation of a large set of transcripts encoding ribosomal subunits²⁹. Hence, KEGG pathway and GO enrichment analysis showed viral gene expression and protein synthesis related terms enriched. Our HPV-negative U-2 OS cellular model highlights the implication of E6 in the downregulation of ribosomal subunits transcripts. In addition, we confirm the upregulation of SYCP2, Nup210 and DICER1 transcripts, which were previously shown to be increased in HPV16-positive clinical samples from cervical²⁹ and head and neck cancers, and HPV-negative cell lines expressing E6 and E7 proteins³⁰,³³.

Comparative RNA sequencing analysis also revealed that E6*I expression was sufficient to modulate the cellular transcriptome. KEGG pathway and GO enrichment annotations showed that E6*I-deregulated transcripts were assigned to extracellular matrix organization, response to decreased oxygen levels, mitogen activated protein
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 7.** Correlation between RNA-sequencing and RT-qPCR results. RT-qPCR analyses of (A) 3 transcripts deregulated in U-E6All cell lines (U-CTL vs U-E6All), (B) the 11 transcripts deregulated in both U-E6All and U-E6*I cell lines (U-CTL vs U-E6*I) and (C) 13 transcripts deregulated in U-E6*I cell lines (U-CTL vs U-E6*I). Fold changes are represented in x- and y-axes for RNA-seq and RT-qPCR, respectively. For each condition, the data represent mean of the 3 clones. For RT-qPCR, RPLP0 was used as housekeeping gene.</figcaption></figure>

kinase (MAPK), focal adhesion, cytokine-cytokine receptor interaction and cancer pathways. STRING analysis indicates that most transcripts implicated in extracellular matrix and response to decreased oxygen levels are predicted to form a significantly enriched core of interactors. A previous study, as well as ROS measurement in our model, indicate that E6*I increases levels of ROS²⁴. Thus, we focused downstream analyses on E6*I-deregulated
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>**Figure 8.** Expression levels of E6- and E6*I-deregulated transcripts in W12 clones. (A) RT-qPCR analysis showing E6 and E6*I absolute mRNA levels (copies per µL of pure cDNA) in W12_20861, W12_20862 and W12_20863 clones. pXJ40-E6 and pXJ40-E6*I diluted in salmon sperm DNA (50 ng/µL) were used to generate standard curves. (B) Western blot analysis showing HPV16 E6 and E7, p53 and pRb protein levels in W12 clones. β-actin is used as loading control. RT-qPCR analyses showing. Full-length gels are presented in Figure S7. (C) SYCP2, TNFAIP6 and p21, and (D) CCL2, RAC2 and PDGFB relative mRNA expression in W12 clones. RPLP0 was used as housekeeping gene. RT-qPCR are represented as means ±/− S.D. of 3 independent experiments.</figcaption></figure>

transcripts implicated in the latter ontology, namely CRYAB<sup>43</sup>, CCL2<sup>44</sup>, CTGF<sup>45</sup>, PDGFB<sup>46,47</sup>, STC1<sup>48</sup>, ARNT2<sup>49</sup>, MMP13<sup>50,51</sup>, MMP3<sup>52,53</sup>, RAC2<sup>54,55</sup> and TNFAIP6<sup>56,57</sup> were also found associated with ROS metabolism. These E6*I-deregulated genes have been validated using RT-qPCR. The high correlations found demonstrate the robustness of fold changes estimated from RNA-seq.

To further investigate the effects of E6*I expression observed in the U-2 OS model, we validated the cellular target genes involved in ROS metabolism in the W12 HPV16-infected keratinocyte cell line. These cells, derived from a low-grade cervical lesion, has been extensively used for the study of HPV infection, life cycle and carcinogenesis. Many clones have been established harboring episomal or integrated type I and II HPV16 genomes to study all those aspects<sup>58</sup>. Interestingly, each clone expresses different levels of E6 and E6*I mRNAs allowing us to get a better understanding on how E6*I potentially modulates the cellular targets identified by RNA-seq and on the existing interplay between E6 and E6*I. Quantification of HPV16 transcripts and proteins in W12 clones revealed a close correlation between mRNA levels and subsequent protein levels as it has been previously reported<sup>59</sup>. TNFAIP6 was found upregulated in both U-E6All and U-E6*I cells by RNA-seq. Nevertheless,
TNFAIP6 RNA levels correlate with E6 expression in the W12. Finally, we confirm that RAC2, CCL2 and PDGFB mRNA levels correlate with E6*I expression in the W12 cellular model. RAC2 is a Rho-like small GTPase, implicated in the activation of the NADPH Oxidase complex-2 (Nox2), catalyzing the production of ROS. It has been shown that ROS production induced by Rac2 leads to DNA damage in leukemia stem cells causing genome instability suspected to drive the rapid accumulation of mutations in progenitor cells⁵⁴. Low levels of RAC2 have also been linked to lower ROS production and DNA damage in quiescent cells exposed to radiation⁵⁵. Interestingly, PDGF growth factor is also closely related to the Nox complex. PDGF interaction with its receptor was shown to activate ROS production by the Nox complex in hepatic⁴⁷ and epithelial model⁴⁶. Finally, CCL2 increases production of ROS by regulating gp91phox, a subunit of the Nox complex⁶⁰. CCL2 expression has also been shown to be reduced when the Nox complex is impaired and produce low levels of ROS⁶¹. The expression of RAC2, CCL2 and PDGFB induced by E6*I could drive Nox2 complex activation, thus representing a new pathway by which E6*I is able to increase cellular ROS production.

Oxidative stress is crucial in carcinogenesis because of its role in inflammation, mutagenesis, initiation and neoplastic progression⁶². Cellular antioxidant defenses exist to maintain cellular homeostasis and viruses are well known to induce oxidative stress by interfering with those defenses⁶³,⁶⁴. Studies have linked HPV16 E6*I expression to changes in the activity of antioxidant enzymes, increased level of ROS and accumulation of DNA damage²⁴,⁶⁵,⁶⁶. In our study, RNA sequencing did not reveal any change in SOD and GPX mRNA levels when E6*I is expressed, probably due to a post-translational regulation of these enzymes by E6*I. In the context of hrHPV-driven carcinogenesis, it has been proposed that oxidative stress induced by E6*I could cause genome instability thus facilitating the integration of HPV genomes into the host cell genome⁶⁷,⁶⁸. This hypothesis is consistent with the correlation between the severity of cervical lesions and the increasing levels of spliced E6*I mRNA⁵,¹²⁻¹⁴,²⁶.

In conclusion, our study reports that HPV16 E6*I induces modulation of the U-2 OS cellular transcriptome using RNA sequencing. We identified cellular genes involved in ROS metabolism that were deregulated by the ectopic expression of E6*I in our model and that are also found deregulated by endogenous E6*I expression in the W12 keratinocyte model. In agreement with previous studies, E6*I can promote viral genome integration into host genome⁶⁹ by inducing ROS and subsequently increasing genome instability²⁴. Here, we highlight the NADPH oxidase pathway as a new potential target of HPV16 E6*I that could explain the increased ROS production. Further studies will be needed to fully understand how oxidative stress related genes modulated by E6*I could affect HPV-driven carcinogenesis and life cycle. This study also demonstrate that E6*I-cellular transcriptome modulation is altered in the presence of E6 suggesting, among others things, a rivalry between the two protein isoforms. The overlap of only 11 transcripts between U-E6All and U-E6*I conditions suggests that E6 is able to partially negate the effects of E6*I. The modulation of ROS production in U-2 OS, and ROS-related genes expression in U-2 OS and W12 clones depending on E6/E6*I levels is in agreement with this interplay. Similarly, it has been observed that the effects of ectopic E6*I on antioxidant enzymes²⁴ seems to be lost when E6*I is co-expressed with E6⁷⁰. Other paradoxical observations were made concerning the expression of E6 and E6*I. The first aspect of this interplay was reported on the opposite effects of both isoforms on p53 degradation and subsequent apoptosis of E6-expressing cells²⁰. More recently, the protection of E6-expressing cells from TNF was only seen when E6 or E6*I are expressed alone but not when both isoforms are expressed³⁷. These observations ask the question of how E6*I can be considered separately of E6 in the natural history of HPV16 infection.

## Materials and Methods

**Cell culture.** Human osteosarcoma derived cell line U-2 OS, HPV16-positive carcinoma derived cell lines SiHa (ATCC® HTB-35™), SCC090 (ATCC® CRL-3239™) and U-2 OS stable cell lines were cultured in DMEM medium (Lonza, Verviers, Belgium). U-2 OS stably transfected clones were routinely cultured in the presence of 625µg/mL G418 to maintain selection. HPV16-positive carcinoma derived cell line Ca Ski (ATCC® CRL-1550™) was cultured in RPMI-1640 medium (Lonza, Verviers, Belgium). All media were supplemented with 10% Fetal Calf Serum (Eurobio, Courtaboeuf, France). Undifferentiated W12 20861, 20862 and 20863 clones were cultured in monolayer on Mitomycin C treated mouse 3T3 cells as previously described⁵⁸. Cell lines were incubated at 37°C in a humidified atmosphere with 5% of CO₂.

**Plasmids, transient and stable transfections.** pXJ40 (empty) and pXJ40-E6All vectors, containing the HPV16 E6 coding sequence, have been generously provided by Murielle Masson (ESBS, Illkirch, France). The pXJ40βGloΔint vector was generated from the pXJ40 vector, using Q5Site-Directed Mutagenesis Kit (New England Biolabs) and the following primers 5'-CTCCTGGGCAACGTG-3' and 5'-CCTGAAGTTCTCAGGATCG-3'. The pXJ40-E6*I vector, lacking the intron 1 located in the E6 ORF corresponding to the nucleotides from 227 to 408 was generated from the pXJ40-E6All vector (Primers: 5'-GTGTATTAACTGTCAAAAGCC-3' and 5'-CTCACGTCGCAGTAACTG-3'). Transient transfections of the U-2 OS cell line were performed using JetPEI reagent (Polyplus-transfection, Illkirch, France) according to manufacturer's instructions. Cells were incubated in the presence or absence of 10 µM of proteasome inhibitor MG132 for 16 h and harvested after 48 h incubation post-transfection. U-2 OS stable cell lines were obtained using JetPEI reagent by co-transfecting with a ratio of 1:9 respectively, the pcDNA3.1-2Flag vector (neomycin resistance) and corresponding pXJ40 vectors (pXJ-CTL, pXJ40-E6All and pXJ40-E6*I). G418 resistant stable clones were then selected after 3 weeks under selection with 800 µg/mL G418 (Euromedex).

**Cell proliferation and clonogenicity assays.** Cell proliferation assays were done with 10,000 cells in 12-well plates by counting cells each day for 4 days using the EVE™ automatic cell counter (Nano EnTek). Clonogenicity assays were performed by seeding 500 cells of each U-CTL, U-E6All and U-E6*I cell lines in 10 cm dishes. After 14 days of culture, colony were counted after crystal violet staining. Clonogenicity radio sensitivity
assays were performed by seeding 300 cells in 6-well plates 24 h before gamma irradiation (0 to 8 Gy) using Gammacell® 40 Exactor (Best Theratronics, Ltd). After treatment, cells were cultured for 10 days then stained with crystal violet. Cell colonies presenting >50 cells were finally counted. For each biological replicate (n = 3), 2 wells were analyzed for colonies counting.

**Measurement of Hydrogen peroxide (H₂O₂).** Measurement of H₂O₂ released from stably transfected U-2 OS cells was carried using the Amplex® Red Hydrogen Peroxide/Peroxidase Assay Kit (Invitrogen) following manufacturer instructions. Four hours after seeding 800,000 cells of each clone in 35 mm dishes, cells were incubated for 1 h in 500 µL PBS1X. Then, 50 µL was transferred in an opaque 96 well plate in triplicate and mixed with an equal amount of Amplex Red reagent (50 µM Amplex Red and 0.1 U/mL HRP final concentrations). After 1 h incubation at room temperature in the dark, fluorescence was measured (560ex/590em nm) on a Synergy™ H1 microplate reader (Biotek).

**Immunoblot assays.** U-2 OS, SiHa, Ca Ski or W12 cells were lysed using RIPA buffer (50 mM Tris-HCl pH7.4, 150 mM NaCl, 1% NP40, 0.5% Na deoxycholate, 1 mM EDTA) supplemented with 30 µg/mL of anti-protease mixture (Roche Diagnostics) followed by sonication. After centrifugation at 10,000 × g for 10 min at 4°C, protein concentrations were determined in supernatant using Bio-Rad Protein Assay (Bio-Rad) according to manufacturer's instructions. Forty µg of protein in Laemml ysis buffer were then loaded onto a sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) gel. After the transfer of proteins onto 0.2 µm Amersham™ Hybond™ PVDF membrane (GE Healthcare) and blocking 1 hour with 5% nonfat milk in Tris-buffered Saline-Tween20, membranes were probed with primary antibodies for 2 hours (except cell signaling antibodies on night) then with HRP-conjugated secondary antibodies (see below) for 1 h. Chemiluminescent signal was detected using PierceECL 2 Western Blotting Substrate (ThermoScientific).

**Antibodies.** Anti-HPV16E6 antibodies 1E-6F4 (E6 N-ter) and 2E-3F8 (E6 C-ter) were obtained from Euromedex. Anti-p53 (DO-1) and anti-pRb (4H1) antibodies were obtained from Cell Signaling. Anti-β-actin (AC-15) antibody was obtained from Sigma-Aldrich. Anti-HPV16E7 (NM2) antibody was obtained from Santa Cruz Biotechnology. HRP-conjugated goat anti-mouse secondary antibody was purchased from BD Pharmingen. Primary and secondary antibodies have been diluted respectively to 1:1000 and 1:5000 in blocking buffer (except for β-actin at 1:20,000 and 1:30,000).

**RT-PCR and RT-qPCR.** Total RNAs were extracted using Ribozol reagent (Amresco) and then 500 µg were retro-transcribed into cDNA using Maxima first strand cDNA synthesis kit (ThermoScientific) according to manufacturer instructions. For PCR, cDNAs were amplified using DreamTaq polymerase (Thermo Scientific) according to manufacturer's instructions. PCR products were analyzed on a 2% (w/v) agarose gel prepared with Tris-Borate-EDTA buffer. Quantitative PCR was performed using Power SYBR Green PCR mix (Life Technologies) according to manufacturer instructions and analyzed on a StepOnePlus Real-Time PCR system (Applied Biosystems). Serial dilutions of pXJ40-E6 and pXJ40-E6*I from 10⁷ to 10¹ copies/µL diluted in salmon sperm DNA (50 ng/µL) (Invitrogen) were used as template to generate standard curve for absolute quantification of E6 and E6*I mRNA copies respectively in W12 clones. Primer sequences used are listed in Supplementary Table 3. RPLP0 expression levels was used as internal normalization standards. The 2⁻ΔΔCt method was used for relative mRNA quantification.

**RNA sequencing.** RNA-seq was carried out to identify differentially expressed transcripts in U-E6All or U-E6*I cell lines compared to U-CTL cell lines. For each clone, total RNAs were extracted using Ribozol reagent (Amresco) according to manufacturer's instructions. cDNA libraries were produced with 1 µg of total RNA using TrueSeq Stranded mRNA Sample preparation kit (Illumina, Inc.). Sequencing was carried on 75 bp in single read on the NextSeq500 platform using the NextSeq500 High Output v2 kit (Illumina, Inc.).

**Reads alignment.** Raw reads quality was assessed using the fastQC software v0.11.3 (Babraham Institute, UK). Reads from each sample were then aligned to the reference human genome (GRCh38, release 87) using TopHat v2.1.1 (John Hopkins University, USA)⁷¹ with a minimum mapping quality score of 30. For each sample, an average of 45 million reads were uniquely mapped.

**Differential expression analysis of cellular transcripts deregulated by E6 and/or E6*I.** Cufflinks (Trapnell Laboratory, University of Washington, USA) was used for transcripts assembly and to assess transcripts abundance, which was normalized between all 9 samples using the fragments per kilobase per million mapped reads (FPKM) method to correct transcript length and library size bias⁷². Differentially expressed transcripts between control condition U-CTL and U-E6All or U-E6*I conditions were finally identified using Cuffdiff (Trapnell Laboratory). Transcripts with a false discovery rate (FDR) < 0.05 were considered significantly deregulated. No fold change cutoff was applied to the data. Nevertheless, more than 95% of transcripts with FDR < 0.05 had a fold change >2.

**Accession numbers.** RNA sequencing data in the form of BAM files have been submitted to SRA@ncbi.nih.gov under the BioProject PRJNA421890 and BioSamples SAMN08159315 (U-CTL #par); SAMN08159316 (U-CTL #1); SAMN08159317 (U-CTL #3); SAMN08159318 (U-E6All #1); SAMN08159319 (U-E6All #6); SAMN08159320 (U-E6All #7); SAMN08159321(U-E6*I #1); SAMN08159323 (U-E6*I #6); SAMN08159324 (U-E6*I #14).
**Functional analysis of deregulated genes.** Enrichment analyses were performed using Gene Ontology (GO) database v1.2 on PANTHER v13.0<sup>73</sup>, Kyoto Encyclopedia of Genes and Genomes (KEGG) release 81.0<sup>74</sup> and STRING v10.5<sup>75</sup>.

**Statistics.** Statistical analyses were carried out using XLSTAT v2017.4 (Addinsoft, Paris, France). Pearson's correlation coefficients (r) were used to investigate correlation between RNA sequencing and RT-qPCR analyses. Mann-Whitney tests were used to analyze statistical significance of proliferation, clonogenicity and H₂O₂ assays, and RT-qPCR experiments (*p* < 0.05; **p** < 0.01; ***p** < 0.001). Proliferation, clonogenicity, H₂O₂ and RT-qPCR graphs represent the mean of 3 experiments and the error bars represent the standard deviation.

**References**

1. Schiffman, M. et al. The role of human papillomavirus (HPV) genotyping in cervical cancer screening: A large-scale evaluation of the cobas HPV test. *Cancer Epidemiol. Biomark. Prev. Publ. Am. Assoc. Cancer Res. Cosponsored Am. Soc. Prev. Oncol.* **24**, 1304–1310 (2015).

2. de Martel, C., Plummer, M., Vignat, J. & Franceschi, S. Worldwide burden of cancer attributable to HPV by site, country and HPV type. *Int. J. Cancer* **141**, 664–670 (2017).

3. Graham, S. V. Keratinocyte Differentiation-Dependent Human Papillomavirus Gene Regulation. *Viruses* **9**, 245 (2017).

4. Johansson, C. & Schwartz, S. Regulation of human papillomavirus gene expression by splicing and polyadenylation. *Nat. Rev. Microbiol.* **11**, 239–251 (2013).

5. Chen, J. et al. Mapping of HPV transcripts in four human cervical lesions using RNAseq suggests quantitative rearrangements during carcinogenic progression. *Virology* **462**, 14–24 (2014).

6. Moody, C. A. & Laimins, L. A. Human papillomavirus oncoproteins: pathways to transformation. *Nat. Rev. Cancer* **10**, 550 (2010).

7. McLaughlin-Drubin, M. E. & Münger, K. The human papillomavirus E7 oncoprotein. *Virology* **384**, 335–344 (2009).

8. Dyson, N., Howley, P. M., Münger, K. & Harlow, E. The Human Papilloma Virus-16 E7 Oncoprotein is Able to Bind to the Retinoblastoma Gene Product. *Science* **243**, 934–937 (1989).

9. Werness, B. A., Levine, A. J. & Howley, P. M. Association of human papillomavirus types 16 and 18 E6 proteins with p53. *Science* **248**, 76–80 (1990).

10. Smotkin, D., Prokoph, H. & Wettstein, F. O. Oncogenic and nononcogenic human genital papillomaviruses generate the E7 mRNA by different mechanisms. *J. Virol.* **63**, 1441–1447 (1989).

11. Cornelissen, M. T. E. et al. Uniformity of the Splicing Pattern of the E6/E7 Transcripts in Human Papillomavirus Type 16-transformed Human Fibroblasts, Human Cervical Premalignant Lesions and Carcinomas. *J. Gen. Virol.* **71**, 1243–1246 (1990).

12. Kösel, S., Burggraf, S., Engelhardt, W. & Olgemöller, B. Increased levels of HPV16 E6*I transcripts in high-grade cervical cytology and histology (CIN II+) detected by rapid real-time RT-PCR amplification. *Cytopathology* **18**, 290–299 (2007).

13. Cricca, M. et al. Molecular analysis of HPV 16 E6I/E6II spliced mRNAs and correlation with the viral physical state and the grade of the cervical lesion. *J. Med. Virol.* **81**, 1276–1282 (2009).

14. Lin, K. et al. E6-associated transcription patterns in human papilloma virus 16-positive cervical tissues. *Oncol. Lett.* **9**, 478–482 (2015).

15. Hong, A. et al. E6 viral protein ratio correlates with outcomes in human papillomavirus related oropharyngeal cancer. *Cancer Biol. Ther.* **17**, 181–187 (2016).

16. Tang, S., Tao, M., McCoy, J. P. & Zheng, Z.-M. The E7 Oncoprotein Is Translated from Spliced E6*I Transcripts in High-Risk Human Papillomavirus Type 16- or Type 18-Positive Cervical Cancer Cell Lines via Translation Reinitiation. *J. Virol.* **80**, 4249–4263 (2006).

17. Zheng, Z.-M. Regulation of alternative RNA splicing by exon definition and exon sequences in viral and mammalian gene expression. *J. Biomed. Sci.* **11**, 278–294 (2004).

18. Stacey, S. N. et al. Leaky scanning is the predominant mechanism for translation of human papillomavirus type 16 E7 oncoprotein from E6/E7 bicistronic mRNA. *J. Virol.* **74**, 7284–7297 (2000).

19. del Moral-Hernández, O. et al. The HPV-16 E7 oncoprotein is expressed mainly from the unspliced E6/E7 transcript in cervical carcinoma C33-A cells. *Arch. Virol.* **155**, 1959–1970 (2010).

20. Pim, D., Massimi, P. & Banks, L. Alternatively spliced HPV-18 E6* protein inhibits E6 mediated degradation of p53 and suppresses transformed cell growth. *Oncogene* **15**, 257–264 (1997).

21. Storrs, C. H. & Silverstein, S. J. PAT, a Tight Junction-Associated PDZ Protein, Is a Novel Degradation Target of High-Risk Human Papillomavirus E6 and the Alternatively Spliced Isoform 18 E6*. *J. Virol.* **81**, 4080–4090 (2007).

22. Pim, D., Tomaić, V. & Banks, L. The Human Papillomavirus (HPV) E6* Proteins from High-Risk, Mucosal HPVs Can Direct Degradation of Cellular Proteins in the Absence of Full-Length E6 Protein. *J. Virol.* **83**, 9863–9874 (2009).

23. Wanichwatanadecha, P., Sirisrimangkorn, S., Kaewprag, J. & Ponglikitmongkol, M. Transactivation activity of human papillomavirus type 16 E6*I on aldo-keto reductase genes enhances chemoresistance in cervical cancer cells. *J. Gen. Virol.* **93**, 1081–1092 (2012).

24. Williams, V. M., Filippova, M., Filippov, V., Payne, K. J. & Duerksen-Hughes, P. Human Papillomavirus Type 16 E6* Induces Oxidative Stress and DNA Damage. *J. Virol.* **88**, 6751–6761 (2014).

25. Filippova, M. et al. The small splice variant of HPV16 E6, E6N, reduces tumor formation in cervical carcinoma xenografts. *Virology* **450–451**, 153–164 (2014).

26. Olmedo-Nieva, L., Muñoz-Bello, J. O., Contreras-Paredes, A. & Lizano, M. The Role of E6 Spliced Isoforms (E6*) in Human Papillomavirus-Induced Carcinogenesis. *Viruses* **10** (2018).

27. Liu, Y., Tergaonkar, V., Krishna, S. & Androphy, E. J. Human Papillomavirus Type 16 E6-enhanced Susceptibility of L929 Cells to Tumor Necrosis Factor α Correlates with Increased Accumulation of Reactive Oxygen Species. *J. Biol. Chem.* **274**, 24819–24827 (1999).

28. Shin, H.-J. et al. Human papillomavirus 16 E6 increases the radiosensitivity of p53-mutated cervical cancer cells, associated with up-regulation of aurora A. *Int. J. Radiat. Biol.* **86**, 769–779 (2010).

29. den Boon, J. A. et al. Molecular transitions from papillomavirus infection to cervical precancer and cancer: Role of stromal estrogen receptor signaling. *Proc. Natl. Acad. Sci. USA* **112**, E3255–E3264 (2015).

30. Pyeon, D. et al. Fundamental Differences in Cell Cycle Deregulation in Human Papillomavirus-Positive and Human Papillomavirus-Negative Head/Neck and Cervical Cancers. *Cancer Res.* **67**, 4605–4619 (2007).

31. Masterson, L. et al. Deregulation of SYCP2 predicts early stage human papillomavirus-positive oropharyngeal carcinoma: A prospective whole transcriptome analysis. *Cancer Sci.* **106**, 1568–1575 (2015).

32. Rajkumar, T. et al. Identification and validation of genes involved in cervical tumourigenesis. *BMC Cancer* **11**, 80 (2011).

33. Harden, M. & Munger, K. Perturbation of DROSHA and DICER expression by Human Papillomavirus 16 Oncoproteins. *Virology* **507**, 192–198 (2017).

34. Isok-Paas, H., Männik, A., Ustav, E. & Ustav, M. The transcription map of HPV11 in U2OS cells adequately reflects the initial and stable replication phases of the viral genome. *Virol. J.* **12**, 59 (2015).

35. Scheffner, M., Werness, B. A., Huibregtse, J. M., Levine, A. J. & Howley, P. M. The E6 oncoprotein encoded by human papillomavirus types 16 and 18 promotes the degradation of p53. *Cell* **63**, 1129–1136 (1990).
36. Jung, A. C. et al. Biological and clinical relevance of transcriptionally active human papillomavirus (HPV) infection in oropharynx squamous cell carcinoma. *Int. J. Cancer* **126**, 1882–1894 (2010).

37. Filippova, M. et al. Complexes of Human Papillomavirus Type 16 E6 Proteins Form Pseudo-Death-Inducing Signaling Complex Structures during Tumor Necrosis Factor-Mediated Apoptosis. *J. Virol.* **83**, 210–227 (2009).

38. Czech, B. & Hannon, G. J. Small RNA sorting: matchmaking for Argonautes. *Nat. Rev. Genet.* **12**, 19 (2011).

39. Sedman, S. A. et al. The full-length E6 protein of human papillomavirus type 16 has transforming and trans-activating activities and cooperates with E7 to immortalize keratinocytes in culture. *J. Virol.* **65**, 4860–4866 (1991).

40. Ristriani, T., Fournane, S., Orfanoudakis, G., Travé, G. & Masson, M. A single-codon mutation converts HPV16 E6 oncoprotein into a potential tumor suppressor, which induces p53-dependent senescence of HPV-positive HeLa cervical cancer cells. *Oncogene* **28**, 762–772 (2009).

41. Kessis, T. D. et al. Human papillomavirus 16 E6 expression disrupts the p53-mediated cellular response to DNA damage. *Proc. Natl. Acad. Sci.* **90**, 3988–3992 (1993).

42. Klymenko, T. et al. RNA-Seq Analysis of Differentiated Keratinocytes Reveals a Massive Response to Late Events during Human Papillomavirus 16 Infection, Including Loss of Epithelial Barrier Function. *J. Virol.* **91**, e01001–17 (2017).

43. van de Schootbrugge, C. et al. Effect of hypoxia on the expression of αB-crystallin in head and neck squamous cell carcinoma. *BMC Cancer* **14**, 252 (2014).

44. Li, X. et al. A CCL2/ROS autoregulation loop is critical for cancer-associated fibroblasts-enhanced tumor growth of oral squamous cell carcinoma. *Carcinogenesis* **35**, 1362–1370 (2014).

45. Murphy-Marshman, H. et al. Antioxidants and NOX1/NOX4 inhibition blocks TGFβ1-induced CCN2 and α-SMA expression in dermal and gingival fibroblasts. *PLoS ONE* **12** (2017).

46. Chen, K. C.-W., Zhou, Y., Zhang, W. & Lou, M. F. Control of PDGF-induced reactive oxygen species (ROS) generation and signal transduction in human lens epithelial cells. *Mol. Vis.* **13**, 374–387 (2007).

47. Adachi, T. et al. NAD(P)H oxidase plays a crucial role in PDGF-induced proliferation of hepatic stellate cells. *Hepatology* **41**, 1272–1281 (2005).

48. Ohkouchi, S., Kikuchi, T., Ebina, M. & Nukiwa, T. Novel functions of stanniocalcin-1 (STC1) through uncoupling protein 2 (UCP2) up-regulation; promoting survival of cancer cells under oxidative stress and inducing the uncoupling respiration (Warburg effect). *Eur. Respir. J.* **38**, p1490 (2011).

49. Keith, B., Adelman, D. M. & Simon, M. C. Targeted mutation of the murine arylhydrocarbon receptor nuclear translocator 2 (Arnt2) gene reveals partial redundancy with Arnt. *Proc. Natl. Acad. Sci. USA* **98**, 6692–6697 (2001).

50. Waldron, A. L. et al. Oxidative stress-dependent MMP-13 activity underlies glucose neurotoxicity. *J. Diabetes Complications* **32**, 249–257 (2018).

51. Del Carlo, M., Schwartz, D., Erickson, E. A. & Loeser, R. F. The Endogenous Production of Reactive Oxygen Species Is Required For Stimulation of Human Articular Chondrocyte MMP Production by Fibronectin Fragments. *Free Radic. Biol. Med.* **42**, 1350–1358 (2007).

52. Radisky, D. C. et al. Rac1b and reactive oxygen species mediate MMP-3-induced EMT and genomic instability. *Nature* **436**, 123–127 (2005).

53. Choi, D.-H. et al. Matrix Metalloproteinase-3 Causes Dopaminergic Neuronal Death through Nox1-Regenerated Oxidative Stress. *PLoS ONE* **9** (2014).

54. Nieborowska-Skorska, M. et al. Rac2-MRC-cIII-generated ROS cause genomic instability in chronic myeloid leukemia stem cells and primitive progenitors. *Blood* **119**, 4253–4263 (2012).

55. Pei, H. et al. RAC2-P38 MAPK-dependent NADPH oxidase activity is associated with the resistance of quiescent cells to ionizing radiation. *Cell Cycle* **16**, 113–122 (2017).

56. Watanabe, R. et al. Emerging Roles of Tumor Necrosis Factor-Stimulated Gene-6 in the Pathophysiology and Treatment of Atherosclerosis. *Int. J. Mol. Sci.* **19** (2018).

57. Shkolnik, K. et al. Reactive oxygen species are indispensable in ovulation. *Proc. Natl. Acad. Sci. USA* **108**, 1462–1467 (2011).

58. Jeon, S., Allen-Hoffmann, B. L. & Lambert, P. F. Integration of human papillomavirus type 16 into the human genome correlates with a selective growth advantage of cells. *J. Virol.* **69**, 2989–2997 (1995).

59. Scarpini, C. G., Groves, I. J., Pett, M. R., Ward, D. & Coleman, N. Virus transcript levels and cell growth rates after naturally occurring HPV16 integration events in basal cervical keratinocytes. *J. Pathol.* **233**, 281–293 (2014).

60. Chun, E. et al. CCL2 Promotes Colorectal Carcinogenesis by Enhancing Polymorphonuclear Myeloid-Derived Suppressor Cell Population and Function. *Cell Rep.* **12**, 244–257 (2015).

61. Chen, X.-L., Zhang, Q., Zhao, R. & Medford, R. M. Superoxide, H2O2, and iron are required for TNF-α-induced MCP-1 gene expression in endothelial cells: role of Rac1 and NADPH oxidase. *Am. J. Physiol.-Heart Circ. Physiol.* **286**, H1001–H1007 (2004).

62. Reuter, S., Gupta, S. C., Chaturvedi, M. M. & Aggarwal, B. B. Oxidative stress, inflammation, and cancer: how are they linked? *Free Radic. Biol. Med.* **49**, 1603–1616 (2010).

63. Paracha, U. Z. et al. Oxidative stress and hepatitis C virus. *ViroJ. Virol.* **10**, 251 (2013).

64. Reshi, M. L., Su, Y.-C. & Hong, J.-R. RNA Viruses: ROS-Mediated Cell Death. *International Journal of Cell Biology*, https://doi.org/10.1155/2014/467452 (2014).

65. Mouret, S. et al. E6* oncoprotein expression of human papillomavirus type-16 determines different ultraviolet sensitivity related to glutathione and glutathione peroxidase antioxidant defence. *Exp. Dermatol.* **14**, 401–410 (2005).

66. Evans, W. et al. Overexpression of HPV16 E6* Alters β-Integrin and Mitochondrial Dysfunction Pathways in Cervical Cancer Cells. *Cancer Genomics. Proteomics* **13**, 259–273 (2016).

67. Williams, V. M., Filippova, M., Soto, U. & Duerksen-Hughes, P. J. HPV-DNA integration and carcinogenesis: putative roles for inflammation and oxidative stress. *Future Virol.* **6**, 45–57 (2010).

68. Camini, F. C., da Silva Caetano, C. C., Almeida, L. T. & de Brito Magalhães, C. L. Implications of oxidative stress on viral pathogenesis. *Arch. Virol.* **162**, 907–917 (2017).

69. Chen Wongworawat, Y., Filippova, M., Williams, V. M., Filippov, V. & Duerksen-Hughes, P. J. Chronic oxidative stress increases the integration frequency of foreign DNA and human papillomavirus 16 in human keratinocytes. *Am. J. Cancer Res.* **6**, 764–780 (2016).

70. Cruz-Gregorio, A. et al. Human Papillomavirus Types 16 and 18 Early-expressed Proteins Differentially Modulate the Cellular Redox State and DNA Damage. *Int. J. Biol. Sci.* **14**, 21–35 (2018).

71. Trapnell, C., Pachter, L. & Salzberg, S. L. TopHat: discovering splice junctions with RNA-Seq. *Bioinformatics* **25**, 1105–1111 (2009).

72. Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L. & Wold, B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. *Nat. Methods* **5**, 621–628 (2008).

73. Mi, H. et al. PANTHER version 11: expanded annotation data from Gene Ontology and Reactome pathways, and data analysis tool enhancements. *Nucleic Acids Res.* **45**, D183–D189 (2017).

74. Ogata, H. et al. KEGG: Kyoto Encyclopedia of Genes and Genomes. *Nucleic Acids Res.* **27**, 29–34 (1999).

75. Szklarczyk, D. et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. *Nucleic Acids Res.* **45**, D362–D368 (2017).
**Acknowledgements**

We thank Muriel Masson for kindly providing the pXJ40-E6All vector, Julie Frejaville and Anne Peigney for technical assistance, Gilles Travé for helpful advice and discussions and Nathalie Dourdin for carefully reading the manuscript. We also wish to thanks Paul F. Lambert for kindly providing W12 cell clones. K.M. and P.P. had fellowships from the French ministry of research and technology. J.P. had fellowships from the Conseil Régional de Franche-Comté. This work was supported by Cancéropôle Est (CE), Ligue Contre le Cancer (Comité du Doubs and CCIR-GE), Université de Franche-Comté (UFC) and Région Franche-Comté.

**Author Contributions**

A.B. conceived the project. P.P.B. performed the bioinformatics analysis, the biochemical assays and the functional assays, K.M. realized the cloning and established the stable cell lines. D.B. performed radio sensitivity tests, A.J. provided the clinical samples, P.P.B. K.M., D.B., J.P., M.H., C.M., J.L.P. and A.B. analyzed the data. P.P.B. and A.B. wrote the paper.

**Additional Information**

**Supplementary information** accompanies this paper at https://doi.org/10.1038/s41598-019-42393-6.

**Competing Interests:** The authors declare no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_32.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2019