Transcriptome analysis of HPV16-positive cells expressing intrabodies targeting E6 and E7 oncoproteins to unravel intrabody antitumor activity

Abstract

Background

Human Papillomavirus-associated cancer remains a global health issue despite the availability of prophylactic vaccination. Two single-chain recombinant antibodies specific for the high-risk HPV16-E6 and E7 oncoproteins were previously developed as intrabodies in HPV16-positive cells. The anti-E6 I7NUC and anti-E7 43M2SD intrabodies demonstrated antiproliferative efficacy in vitro and in mouse tumor models by interfering with p53 and pRb tumor suppressors. This study aimed to explore the mechanism behind this antitumor activity.

Methods

HPV16-positive SiHa cells were electroporated with the intrabody DNA plasmids and processed for RNA-seq analysis 6 h after transfection. Differential expression and functional analysis were conducted to identify significantly modulated genes and enriched pathways. RT-qPCR was employed to confirm the modulation of selected genes. The human interactome of E6 and E7, retrieved from BioGRID, was used to check for known interactions.

Results

The anti-E6 I7Nuc modulates 12,820 genes, of which 6,072 are upregulated and 6,748 downregulated, mainly encoding proteins that directly bind to HPV16E6. The anti-E7 43M2SD intrabody modulates 1,174 genes, mostly encoding proteins that interact indirectly with E7 binders. Notably, 3,191 genes in cells expressing I7NUC and 468 in those expressing 43M2SD are novel. The functional analysis revealed over 100 KEGG pathways significantly affected by I7NUC and only one pathway affected by 43M2SD. “Neuroactive ligand-receptor interaction, Proteasome, Nucleocytoplasmic transport, Protein processing in endoplasmic reticulum, Cell cycle, ATP-dependent chromatin remodelling, Cellular senescence, Ubiquitin-mediated proteolysis, Ribosome and Viral carcinogenesis” are the top significant pathways affected by I7NUC, with almost 100% down-modulated DEGs. 43M2SD affected only the “Viral protein interaction with cytokine and cytokine receptor” pathway, with CXCL6, 9, 10 and 11 genes all downregulated. I7NUC determines the transcriptional downregulation of Energy Metabolism pathway enzymes, including SLC2A1, SLC5A1, hexokinase1/2, LDHA, GLS2, SLC7A5, SLC7A11, HIF-1, and a significant percentage of transcripts associated with the phosphorylation-oxidative enzyme complexes.

Conclusions

The anti-E6 and anti-E7 intrabodies interfere with cellular pathways essential for tumor transformation, confirming their great potential as therapeutic agents for HPV16-related diseases. The differences observed in the effects of the two intrabodies at the transcriptional level may be advantageous in clinical applications to different tumor types or grades of premalignant lesions.

Background

Human Papillomavirus (HPV) infection is closely associated with the development of precancerous and cancerous lesions in various parts of the body. Twelve high-risk (HR) genotypes (HPV16, 18, 31, 33, 35, 39, 45, 51, 52, 56, 58 and 59) of the alpha genus are responsible for nearly all cervical cancers, most cancers in the anogenital area, and a growing percentage of Head and Neck Squamous Cervical Carcinoma (HNSCC) [1, 2]. This occurs despite efficacious preventive vaccines being available, due to the incomplete global vaccination coverage [3]. Therefore, it is crucial to develop safe and effective treatments for HPV-related lesions and cancer, especially when caused by genotypes not included in the current vaccine, and for unvaccinated individuals with recurrent HPV-related infections, people with immune deficiencies, people living with HIV, and people not responding to vaccination [4].

We have been working for many years on a potential therapy based on recombinant antibodies targeting the E6 and E7 oncoproteins of the HR HPV16 (16E6 and 16E7, respectively), which alone account for more than 60% of cervical cancers and nearly all (90%) HPV-associated HNSCCs [5]. We developed one anti-16E6 antibody in single-chain format (scFv), I7NUC, and one anti-16E7 scFv, 43M2SD, which we expressed in HPV16-positive cells as intracellular antibodies (intrabodies) [6, 7]. Since the E6 and E7 oncoproteins interact with multiple binders as they shuttle between the cytoplasm and the nucleus [8], we tested in which compartment the two intrabodies were most effective at blocking the oncoproteins. The intrabody coding sequences were cloned into the plasmids of the ScFvExpress series, which carry signals for the expression in the nucleus, endoplasmic reticulum (ER) and secretory compartment [9]. In our studies, we found that the most effective anti-E6 intrabody was I7NUC, carrying a nuclear localization signal (NLS), while the most effective anti-E7 intrabody was 43M2SD, carrying the signal peptide for the secretory compartment and the SEKDEL signal for retention in the ER [10].

The two intrabodies have proven to be effective in hampering cell proliferation and prompting apoptosis in vitro and in preventing or delaying the tumor progression of HPV tumors in preclinical models [6, 11]. The intrabodies can bind to their protein targets within cells, and their activity derives in part from their capacity to hinder the binding of the oncoproteins to their major intracellular targets, p53 for E6 and pRb for E7 [12].

The anti-16E6 I7NUC causes partial rescue of p53, resulting in increased cell death by necrosis and/or apoptosis [6]. The anti-E7 43M2SD binds to the E7 N-terminus, where inhibition of HPV-positive SiHa cell proliferation partially depends on the interference with the E7 binding to pRb [13]. The anti-E7 intrabody sequesters E7 in the ER compartment, where the oncoprotein accumulates [11]. The anti-16E6 intrabody accumulates in the nucleus. The histological observation of tumors treated with the two therapeutic intrabodies shows large areas of necrosis, which is instead absent in tumors treated with irrelevant scFvs. Staining of tumor sections for active caspase-3, a major player in the apoptotic process, revealed a high percentage of caspase-3 positive areas in tumors treated with therapeutic scFvs compared to controls [10].

After these investigations on the intrabody mechanism of action, we were interested in clarifying the pathways involved in tumor cells expressing therapeutic intrabodies. This is crucial in view of future applications in clinics and could be useful to identify other targets potentially important for designing novel therapeutic interventions for HPV-associated tumors. Furthermore, knowing the different targets, whether overlapping or not, of intrabodies could allow us to investigate their combined and possibly synergistic use.

Here, we report transcriptomic analysis of HPV16-positive SiHa cells expressing the anti-16E6 I7NUC or the anti-16E7 43M2SD intrabody. Through RNA-Seq and bioinformatics analysis, we identified Differentially Expressed Genes (DEGs) in intrabody-expressing cells in comparison to their control cells expressing the same intracellular localization signal present in the intrabody. We compared the results from transcriptomic analysis with those obtained by RT-qPCR of selected genes. We also performed functional analysis to identify KEGG pathways influenced by the intrabodies and investigated the interactions of proteins encoded by DEGs with the currently known interactome of E6 and E7 proteins. Lastly, we discuss the implications of dysregulation of selected pathways supportive of antitumor intrabody activity.

Materials and methods

Plasmids, cells and electroporation

The anti-16E6 scFv I7NUC, anti-16E7 scFv 43M2SD, scFvExnuc and scFvExSD DNA plasmids [7, 14] were produced in E. coli JM109 cells (PROMEGA), and purified using the EndoFree Plasmid Maxi Kit (QIAGEN) according to the manufacturer’s protocol; DNA was quantified by NanoDrop (Thermo Scientific™) and analysed by agarose gel electrophoresis after Not1 enzymatic restriction for monitoring the integrity.

Human cervical carcinoma SiHa cells (ATCC HTB35) were grown at 37 °C in a humidified atmosphere with 5% CO₂ in Dulbecco’s modified Eagle’s medium (DMEM) containing 10% fetal calf serum (FCS), 100 units/ml penicillin, 100 µg/mL streptomycin, and 2 mM glutamine (GIBCO, GAITHERSBURG, MD, USA).

Cells at 70% confluence were used for electroporation experiments, using the Amaxa^® Nucleofector I instrument and the Cell Line Nucleofector^® Kit V (LONZA) according to the manufacturer’s recommendations. The cells were detached from the plate using TrypLE Express Enzyme (EE, GIBCO) and then rinsed in complete medium. For each condition, 5 × 10^6 cells were centrifuged and resuspended in 100µL of Cell Line Nucleofector Solution V, and 2 µg of the proper DNA plasmid was added. The cell/DNA suspension was transferred into the supplied sterile cuvette provided in the kit and inserted into the Nucleofector device (LONZA), using the T30 program. After electroporation, the sample was resuspended in 400 µL of preheated culture medium and plated on 35-mm culture dishes containing 1 mL of culture medium. The cells were incubated at 37 °C in a humidified atmosphere with 5% CO₂. After 6 h, the transfection was stopped, and the cells were collected and stored at -80 °C until RNA extraction.

RNA extraction and sequencing

RNA extraction and sequencing were provided as a service by the Laboratory for Advanced Therapy Technologies (LTTA) of the University of Ferrara (Ferrara, Italy). Three distinct transcriptome experiments were conducted using 0.25 µg of total RNA, depleted of ribosomal RNA, extracted from SiHa cells that were electroporated with plasmid DNA expressing anti-E6 16I7NUC or the anti-16E7 43M2SD intrabodies. Parallel transcriptome experiments were carried out with RNA extracts from SiHa cells electroporated with plasmid vectors containing only the intracellular localization signals, which served as a control. Total RNA was extracted from cells using the Maxwell Instrument, an automated nucleic acid purification platform (Promega Corporation, Madison, WI, USA), with the purification kit Maxwell RSC miRNA from Tissue (#AS1460, Promega) according to the manufacturer’s instructions. For RT-qPCR analysis, RNA was purified using the mirVana™ PARIS™ RNA Kit (INVITROGEN) according to the manufacturer’s protocol, which includes the addition of RNase Inhibitor (AMBION) and DNase I (BIOLABS) directly to the Cell Disruption Buffer. RNA was then quantified using a NanoDrop spectrophotometer (Thermo Scientific).

RNA-seq libraries were prepared using the Qiagen QIAseq Fast Select RNA Removal kit for ribosomal RNA depletion and the Qiagen QIAseq Stranded Total RNA Library kit for RNA fragmentation, reverse transcription, second strand synthesis, end-repair, A-addition, and adapters ligation. RNA sequencing was performed according to the Illumina pipeline on a NextSeq 500 Instrument using the NextSeq High-Output kit v2 (150 cycles).

Bioinformatics analysis

Sequence data in FastQ format were checked for quality with FastQC v0.11.9 and then trimmed for adapter and low-quality sequences with Trimmomatic v0.39. Sequences were then mapped versus the human genome (Ensemble GRCh38) using STAR v2.7.9a, counting obtained reads with featureCounts v2.0.3. By RuvSeq v1.24 in R v4.0.4 the batch effect was reduced setting k = 4, next performing differential expression analysis with edgeR v3.32.1, comparing anti-E7 43M2SD- and anti-E6 I7NUC-expressing SiHa versus cells expressing only the respective signals for intracellular localization, and correcting for the False Discovery Rate (FRD < 5%). A threshold of log₂FC > 1 as an absolute value was applied. Overlapping of DEGs was visualized by Venn Diagrams generated with ggVennDiagram v1.2.3 in R, whereas hierarchical clustering was performed with Complex Heatmap v2.20.0 in R.

Functional analysis of enriched KEGG pathways was performed with ClueGo app v.2.5.7 in Cytoscape v3.9.1, setting significance at p < 0.05 for the Benjamini-Hotchberg correction. In Cytoscape, we retrieved from the BioGrid repository the interactome of human proteins binding to HPV proteins, superimposing expression data of DEGs for the two conditions. The STRING app v1.7.1 within Cytoscape was used to retrieve protein-protein interaction networks among proteins encoded by DEGs in the two conditions or by DEGs affected by 43M2SD, and E7 binders retrieved by BioGrid. STRING functional enrichment function was also used to highlight the participation of selected nodes in KEGG pathways by donut visualization. Topological characteristics (Betweeness Centrality, Closeness Centrality and Degree) were calculated by the Network Analyzer tool in Cytoscape. Network size was reduced by keeping nodes having at least two out of three parameters above the 90th percentile for the 43M2SD condition or above the 95th percentile for the I7NUC condition due to the complexity of the network.

The presence of transcripts or proteins for each of the 40 protein-coding DEGs shared by I7NUC and 43M2SD was checked in tissues and cell lines derived from cervical, head, and neck cancers using the Human Protein Atlas, an open-access resource for human proteins (www.proteinatlas.org).

The raw data from the RNA-seq analysis have been uploaded to the NCBI SRA portal and are available using the BioProject ID PRJNA1177673.

Reverse transcriptase quantitative PCR (RT-qPCR)

To confirm transcriptomic results by RT-qPCR, we chose genes that were modulated with at least logFC > 1 as an absolute value. RT-qPCR assays were performed in intrabody-expressing SiHa cells on the ViiA 7 Real-Time PCR System (THERMO FISHER SCIENTIFIC). The expression of the target genes was quantified using the Luna Universal One-Step RT-qPCR kit (BIOLABS), with 20 ng of RNA per reaction. The Luna Universal One-Step RT-qPCR Kit uses real-time SYBR Green I fluorescence dsDNA binding dye to measure DNA amplification after each PCR cycle. The first reverse transcription step was performed at 55 °C for 10 min, followed by denaturation at 95 °C for 1 min and 45 denaturation/extension cycles at 95 °C for 10 s and 60 °C for 1 min. Gene expression levels were assessed using the 2^−ΔΔCT method using housekeeping eukaryotic translation elongation factor 1 Gamma (EEF1G) and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as reference genes for RNAs from I7NUC and 43M2SD expressing cells, respectively [15]. This choice is due to GAPDH being downregulated in the I7NUC transcriptome. The primer pairs used for the analysed genes are provided in Table AF1.

Statistical analysis

To assess statistical significance, the homogeneity of variances between groups was verified using an F-test. Since the calculated F-values consistently fell within the acceptance range at p = 0.05, the assumption of equal variance was considered valid. Consequently, gene expression differences were evaluated using a one-tailed unpaired Student’s t-test assuming equal variances (type 1).

Results

Differentially expressed genes (DEGs) in SiHa cells expressing the anti-E6 I7NUC or the anti-E7 43M2SD intrabody

We investigated differences in the cell transcriptome potentially attributable to the intracellular expression of recombinant antibodies previously developed against 16E6 and 16E7, the main oncoproteins of HPV16. The HPV16-positive SiHa cells were electroporated with plasmids expressing the anti-16E6 I7NUC or the anti-16E7 43M2SD and, in parallel, with plasmids expressing only the respective signals for intracellular localization, Nuclear Localization signal (NLS) and SEKDEL, as controls. Transcriptome analysis revealed a high number of differentially expressed genes (DEGs) in intrabody-expressing cells compared to their controls. In particular, in 43M2SD-expressing cells, we found a total of 1,174 DEGs, of which 574 were upregulated and 600 downregulated; otherwise, in I7NUC-expressing cells, we observed a 10-fold higher number of DEGs, specifically 12,820, of which 6,072 were upregulated and 6,748 were downregulated (Fig. 1, panel A). Notably, 468 genes in 43M2SD- and 3,191 in I7NUC-expressing cells were novel genes that had never been annotated before. Among the DEGs relative to 43M2SD, over a tenth (132) were protein-coding genes, with a similar number of up-and-down-regulated genes. Differently, approximately half of the I7NUC-related DEGs were protein-coding genes, mostly downregulated. As shown in Fig. 1, the number of overlapping DEGs for total annotated transcripts (panel B) and protein-coding genes (panel C) in the two conditions is small. Interestingly, the hierarchical clustering of the 40 shared protein-coding DEGs (Fig. 1, panel D) highlighted some opposite modulations by the I7NUC and 43M2SD intrabodies. Most shared DEGs encode membrane proteins (25/40; 62.25%), among which we find 8 chemoreceptors that are a type of G-protein-coupled receptor belonging to olfactory receptors, the largest gene family in the human genome (400 genes). Of the 25 DEGs showing opposite modulation, 22 are upregulated by I7NUC and downregulated by 43M2SD. The opposite, i.e., the downregulation by I7NUC and the upregulation by 43M2SD, concerns 3 genes, the insulin-like growth factor 2 receptor (IGF2R), the Ankyrin repeat domain 39 (ANKRD39) and the taste 2 receptor member 31 (TAS2R31). The complete list of DEGs detected in the intrabody-expressing cells is in the additional file Table AF2.

Fig. 1

A. Number of differentially expressed genes (DEGs) in HPV16-positive SiHa cells expressing 43M2SD or I7NUC intrabodies. For both total transcripts (TOT) and protein-coding genes, the number of total and up-and-down-regulated genes is indicated. B. Venn Diagrams showing overlapping I7NUC and 43M2SD DEGs considering total annotated transcripts and C) protein-coding genes. D. The hierarchical clustering of the 40 protein-coding DEGs shared by the I7NUC and 43M2SD conditions is shown

Consulting the Human Protein Atlas for the 40 DEGs shared by the two intrabodies, we found that only 24 are present in CC, HNSCC, or cervical cancer-derived cell lines (https://www.proteinatlas.org/, accessed October 1, 2024). The remaining 16 DEGs are uncharacterised in the context of HPV infections and HPV-related diseases. A table with the complete list of the protein-coding DEGs shared by I7NUC and 43M2SD, their description, their up- or down-modulation and their presence or absence in the Human Protein Atlas is shown as an additional file (Table AF3).

We then selected 9 genes identified as differentially expressed for the analysis of their expression levels by RT-qPCR. Specifically, we chose genes that are overexpressed in cervical cancer cells, including PARP14 [16] and GAPDH [17], the WDR83OS gene, potentially involved in cell proliferation [18], the ZNF177 gene, involved in invasiveness [19], and CSF2, CXCL10 and IFNA7, which are involved in inflammation. The transcription of CDKN2A and CDKN1A genes, encoding cell cycle regulatory proteins, was also analyzed. The results of RT-qPCR are shown in comparison with those of DEG analysis in Table 1, and confirm the differential expression for all tested genes, although the results for three of them did not reach statistical significance.

Table 1 DEG analysis and RT-qPCR results for selected genes modulated in SiHa cells expressing either the anti-E6 I7 NUC intrabody (top section) or the anti-E7 43M2SD intrabody, based on at least three independent transfection experiments. The table shows the log2 fold change (Log2 FC) and corresponding p-values from the differential expression analysis (DEG), alongside the relative quantification from RT-qPCR, expressed as –ΔΔCt values with associated p-values. Negative values of Log2 FC or –ΔΔCt indicate downregulation in treated conditions compared to controls; positive values indicate upregulation

Functional enrichment analysis of KEGG pathways affected by the expression of I7NUC and 43M2SD intrabodies

We then performed the functional analysis using the KEGG database to identify pathways significantly enriched in each intrabody condition/treatment.

In I7NUC-expressing cells, we found 107 enriched KEGG pathways, which are shown in an additional file (Table AF4). The top 20 enriched pathways are shown in Fig. 2, panel A, in order of significance. In this list, Neuroactive ligand-receptor interaction is the top significant pathway, but most of the enriched pathways are related to the anti-tumor activity of I7NUC intrabody, being cellular metabolic activity and cellular processes linked to cell proliferation and transformation such as Proteasome, Nucleocytoplasmic transport, Protein processing in endoplasmic reticulum, Cell cycle, ATP-dependent chromatin remodelling, Cellular senescence, Ubiquitin mediated proteolysis, Ribosome, and Viral carcinogenesis. Interestingly, these pathways featured almost 100% down-modulated DEGs.

Fig. 2

A. Top 20 KEGG pathways significantly enriched in I7NUC-expressing cells. Pathways are listed in decreasing order of significance, and bars indicate the total number of featured DEGs as the sum of up- and-down-regulated genes in red and blue, respectively. B. Pathway significantly enriched by 43M2SD expression. The pathway is blue since more than 50% of the featured DEGs are downregulated. Accordingly, downregulated DEGs are in blue

In 43M2SD-expressing cells, only the Viral protein interaction with cytokine and cytokine receptor pathway was enriched, featuring CXCL6, 9, 10 and 11 genes, which were all downregulated as shown in Fig. 2, panel B.

In Fig. 3, we report a network showing the genes shared among the most relevant I7NUC-enriched pathways. With the exception of the Proteasome pathway, the selected pathways shared several genes. The Cellular senescence pathway shared the highest number of genes, specifically 15 DEGs with Viral carcinogenesis and 13 with Cell cycle. In addition, these three pathways shared 12 DEGs, including MDM2, the negative regulator of p53, also shared with Ubiquitin mediated proteolysis. In turn, this pathway shared 9 DEGs with the Protein processing in endoplasmatic reticulum and 2 with the Cell cycle.

Fig. 3

Network of DEGs shared among selected I7NUC-enriched pathways. Pathways in blue mean that more than 50% of the featured DEGs (including those not shown in this figure) are downregulated. Accordingly, downregulated DEGs are in blue

Although Energy metabolism is not among the most prominent pathways in the functional enrichment analysis of DEGs from both intrabodies, it was examined in detail in light of its importance for the activity of the E6 and E7 oncoproteins. In fact, it is well known that E6 and E7 can interfere with several key enzymes involved in glucose transport, glycolysis and oxidative phosphorylation pathways to improve the uptake of nutrients and carbon sources, increasing the production of ATP, NADPH, and FADH2 [20], thus supporting their protumor activity [21]. As reported in Table 2, the expression of I7NUC determines the transcriptional downregulation of both SLC2A1 and SLC5A1 genes coding for the glucose transporter GLUT 1 and the sodium-glucose transporter SGLT 1 proteins, as well as the glycolysis-related enzymes hexokinase 1/2 (HK1/2). The binding of I7NUC to E6 also prevents the increase in LDHA expression induced by E6-mediated p53 degradation [22]. Furthermore, I7NUC-E6 binding prevents the increase in glutaminolysis, lowering GLS2 expression and glutamine uptake by the amino acid transporters Lat-1 (SLC7A5) and XCT (SLC7A11), consequently starving the Krebs cycle. Interestingly, we found in I7NUC-expressing cells a significant downregulation of a subunit of hypoxia-inducible factor 1 (HIF-1), known as an enhancer of transcription, not only of genes involved in the glycolysis pathway but also of genes involved in tumor progression and metastasis [23].

Table 2 Impact of the anti-16E6 (I7NUC) and anti-16E7 (43M2SD) intrabody expression on energy metabolism enzymes. Transcripts that are differentially expressed in the SiHa transcriptome following expression of the I7NUC and 43M2SD intrabodies. The “n” indicates non-differential expression

Lastly, 56 out of 57 protein-coding transcripts encoding the enzymes of the mitochondrial electron transport chain, which are part of the phosphorylation oxidative pathway, were downregulated in the transcriptome of I7NUC-expressing SiHa cells (additional file Table AF5).

Differently, none of the energy metabolism genes was affected by 43M2SD expression except for the upregulation of PKLR, encoding the pyruvate kinase, which catalyses the transphosphorylation of phosphoenolpyruvate into pyruvate and ATP, which is the rate-limiting step of glycolysis.

Interactomes of binders of HPV E6 and E7 and DEG-encoded proteins

To search for evidence that supports the known anti-tumor action of intrabodies, we investigated whether the DEG-encoded proteins represented direct binders of the target proteins of intrabodies. Therefore, we retrieved the human interactome of E6 and E7 from BioGRID [24]. We found that many proteins encoded by I7NUC DEGs are direct binders of E6 (Fig. 4, panel A) and, interestingly, are also direct binders of E7 and other viral oncoproteins such as E5, E1 and E2 (Additional file Figure AF1). The modulated DEG-encoded E6 binders, that are part of the interactome, are all downregulated except Ankyrin Repeat and SOCS Box Containing 4 (ASB4), which are involved in protein degradation. We investigated the pathways in which the E6 binders modulated by I7NUC are involved. We found that they are implicated in Cancer, Apoptosis, TNF signalling and Proteasome pathways and that they share several proteins with some neurological disease pathways (Fig. 4 panel B).

Fig. 4

A. BioGRID interactome demonstrating direct interactions between HPV E6 binders and proteins encoded by I7NUC DEGs. Proteins encoded by I7NUC DEGs are highlighted with shades of red or blue to indicate the up- or down-regulation of the corresponding DEGs, respectively. Grey colour indicates no deregulation. B. E6 binders affected by I7NUC are involved in key cellular pathways. The colour of the pathways reflects the level of significance, while the size indicates the number of featured DEGs

Unlike what was verified for the interactome of I7NUC DEG-encoded proteins, we did not identify any direct E7 binder among the proteins encoded by 43M2SD DEGs using the BioGRID interactome. Therefore, we investigated the possibility of an indirect interaction between the proteins encoded by 43M2SD-modulated DEGs and known E7 binders. The inquiry revealed several interactions with E7 binders showing significant involvement in Viral carcinogenesis, Cell cycle, Cellular senescence, Ubiquitin-mediated proteolysis and HPV infection pathways (Fig. 5).

Fig. 5

Biogrid interactome showing indirect interaction of HPV E7 binders (grey) with proteins encoded by the 43M2SD DEGs; shades of red and blue represent the levels of up-or down-regulation of the corresponding DEGs, respectively. Donuts sections around E7-binding proteins indicate the top five significant KEGG pathways to which they belong, with colours according to the legend

Interaction networks of DEG-encoded proteins

We then tried to identify potential biomarkers, represented by hubs, in the I7NUC and 43M2SD interactomes.

After filtering for topological parameters, the interactome of proteins encoded by the I7NUC-DEGs featured 173 nodes, all of which were downregulated except for two (Fig. 6, panel A). Interestingly, most of the proteins participating in the pathways are transcription factors, such as STAT3, which has the highest number of interactions. Functional enrichment yielded the Viral carcinogenesis pathway as the top significant.

Fig. 6

Protein-protein interaction networks of proteins encoded by DEGs modulated by (A) I7NUC or (B) 43M2SD; shades of red and blue indicate the levels of up-or down-regulation of the respective DEGs. Donuts sections around proteins indicate the top five significant KEGG pathways to which the DEGs belong, with colours according to the legend. C. Intersection of the two interactomes

The interactome for 43M2SD featured 25 nodes, including induced and repressed ones (Fig. 6, panel B). The top three interacting hubs with a degree of 10 were CXCL9, CXCL10 and CXCL11.

To verify possible overlaps, we extracted the intersection of the two interactomes, identifying a core of seven potentially interesting proteins common to the two conditions, including CXCL10 and CXCL11 as the most relevant (Fig. 6, panel C).

Discussion

While the causal relationship between HR HPV E6 and E7 oncoproteins and cancer development is well established, and the pathways involved are broadly studied, we do not yet have a clear picture of the targets to hit to achieve the most efficacious therapy for cervical cancer. In this context, transcriptomic analysis can help unravel the details of pathways leading to cancer evolution, thus favouring the development of therapeutic strategies.

In this perspective, and in light of our previous studies on the efficacy of anti-16E6 I7NUC and anti-16E7 43M2SD intrabodies in combating HPV16-associated tumors, we performed transcriptomic analysis of SiHa cells expressing or not the intrabodies. The aim was to focus on the mechanism underlying the action of intrabodies and to possibly identify new targets for effective antitumor action.

Notable differences were found between the two therapeutic intrabodies in terms of total DEGs, number of protein-encoding genes affected (6045 for I7NUC and 87 for 43M2SD) and direction of this modulation. In fact, the I7NUC protein-encoding DEGs represent approximately half of the total, the majority of which are downregulated, while those related to 43M2SD are only a tenth of the total, with a comparable number of up- and-down-regulations.

The clear numerical difference of DEGs related to the two intrabodies is likely due to their intracellular localization. I7NUC primarily resides in the nucleus, where it co-localizes with the E6 protein [6] and where mRNA synthesis occurs. Conversely, 43M2SD localizes in the ER by virtue of its SEKDEL tag and is able to relocate the E7 oncoprotein in this cytoplasmic compartment [11].

In agreement with I7NUC intrabody targeting, we found that the third most enriched KEGG pathway was Nucleocytoplasmic transport, in which 75.9% of the genes were downregulated. In particular, genes for importin and exportin are downregulated, as well as genes forming the Nup358 complex, the central channel and the symmetrical nups, which might result in a severely reduced translocation through the nuclear pore complexes (NPC) [25]. In addition, the expression of genes involved in the exon-junction complex is downregulated, with a possible impact on the export of mRNA from the nucleus to the cytoplasm and consequent impairment of mRNA surveillance [26]. Nucleocytoplasmic transport is linked to the Ubiquitin mediated proteolysis and the Protein processing in endoplasmic reticulum pathways by two shared genes, UBE2I and SEC13, respectively. These two last pathways, which have in common several DEGs, are strictly related, regulating protein trafficking and degradation. In the Protein processing in endoplasmic reticulum pathway (second top enriched pathway), I7NUC affected 67.65% of the genes, downregulating over 96% of them. Such extensive downregulation severely affects protein translocation, folding, post-translational modifications, shuttling to and from the Golgi compartment, and targeting as well as the unfolded protein response [27]. Interestingly, the caspase 12 (CASP12) is the only gene to be upregulated in this pathway. CASP12 is specifically activated in cells experiencing ER stress (ERS) too severe and prolonged and was found to play a role in ER stress-mediated apoptosis [28]. Furthermore, the impairment of UPR, coupled with the concomitant weakening of the Ubiquitin mediated proteolysis pathway, can block the degradation of misfolded proteins, leading to the harmful accumulation of proteins within the cell [29]. In addition, the Proteasome pathway, which represents the last step of the protein degradation process, is impaired, with I7NUC down-regulating 67.7% of the genes including those involved in the immunoproteasome. Overall, impairment of Ubiquitin mediated proteolysis, Protein processing in endoplasmic reticulum and Proteasome pathways may result in high toxicity and represent a major intrabody weapon to kill the cells that express it. In support of this hypothesis, some recent therapeutic approaches against cancer are based on the inhibition of the proteasome [30]. Noteworthy, the Ubiquitin mediated proteolysis and the Proteasome pathways were significantly enriched in the I7NUC DEGs-encoded direct binders of E6, further indicating these pathways as relevant in the mechanism of I7NUC intrabody action. These results, in addition to confirming the previous observation on the ability of I7NUC to bind E6 [6], lead to the hypothesis that the presence of I7NUC is able to block the interaction of the oncoprotein with its numerous cellular targets, thus interrupting cellular pathways that depend on these bonds.

Specifically, one of the main protumor actions that E6 exerts in HPV-infected cells is to induce degradation of the p53 tumor suppressor by forming a complex with the cellular E3-ubiquitin-protein ligase (UBE3A, also known as E6AP), which then binds to and ubiquinates p53 [8]. The I7NUC expression was found to prevent p53 degradation [6]. By repressing the ubiquitination and the proteasome machinery, the I7 NUC intrabody can support p53 rescue and ultimately restore the tumor suppressor functions. The Viral carcinogenesis pathway appeared to have a central role in the dysregulation induced by the I7NUC intrabody being top enriched also in the interaction network among proteins encoded by DEGs. In addition, 39% of the downregulated genes are implicated in the HPV infection pathway, and several featured DEGs are shared with Cell cycle and Cellular senescence pathways. This confirms a strong alteration of the cell cycle, and one could postulate that cells expressing I7NUC enter senescence less than untreated controls. Indeed, recent studies indicate that the senescence-associated secretory phenotype (SASP), linked in particular to senescent stromal cells, contributes to tumorigenesis by creating an inflammatory environment promoting proliferation and metastasis [31]. Consequently, by reducing senescence, I7NUC may counteract the occurrence of an inflammatory environment and play a role against tumor progression. The broad antitumor effect of I7NUC is also documented by the downregulation of a high number of genes belonging to energy metabolism, which impairs the production of ATP, NADHPH, and FADH2 [20]. Functional analysis on SiHa cells expressing 43M2SD identified only the enrichment of Viral protein interaction with cytokine and cytokine receptor pathway, featuring CXCL6, CXCL9, CXCL10 and CXCL11 chemokine-coding genes, all downregulated. Chemokines are a group of low molecular weight peptides whose main function, together with other cytokines, is to regulate recruitment, activation, phenotype, and function of immune cells by controlling their location and interactions in the tumor microenvironment (TME), a crucial component in cervical cancer biology [32]. Chemokines can also regulate various tumor cell properties, including proliferation and invasiveness, as well as neoangiogenesis in stromal cells, leading to tumor onset and development [33]. Thus, the downregulated expression of these chemokines may be associated with the known antiproliferative and proapoptotic effect of 43M2SD [10]. There are numerous studies to support this hypothesis. In a recent study using public gene expression databases, the comparison of chemokine levels between cervical cancers and adjacent tissues biopsies highlighted that in the tumor tissues, the expression of CXCL1/3/5/6/8/9/10/11/13/16/17 is upregulated while the expression of CXCL12/14 is downregulated [34]. Overexpression of different chemokines, including CXCL10 and CXCL11, was also reported in HPV16-positive SiHa cells and in human keratinocytes transduced with the HPV16 E6/E7 oncoproteins [35]. Accordingly, although CXCL10 was generally considered a chemokine with angiostatic and antitumor properties [36], some evidence suggests that the increased expression of CXCL10 and its corresponding receptor CXCR3 may be associated with advanced human cancers. The tumor-promoting activity of CXCL10 could depend on the splice variant of its receptor CXCR3. It has, in fact, been demonstrated that the isoform known as CXCR3-A, found in most cell types including keratinocytes, coupled to Gαi protein, enhances cell proliferation and invasiveness [37] by activation of the ERK1/2, p38/MAPK, JNK and PI3-kinase/Akt signalling pathways [38, 39]. Furthermore, via microarray analysis, CXCL10 secretion was correlated with upregulated exosomal expression of PD-L1 by fibroblasts in HPV-positive patients, contributing to viral immune escape and carcinogenesis [40]. Similarly, CXCL9 binds to the CXCR3 receptor, and an increasing body of evidence demonstrates that it can act as a tumor promoter in multiple types of cancer, including cervical cancer [41]. In addition, CXCL9 was found among the serum proteins that progressively increase with disease stage in patients with cervical squamous cell carcinoma [42], suggesting that this chemokine could be involved in the development of this pathology. Also, higher expression levels of CXCL11 were found in tumor tissue compared to adjacent tissues in CC samples, and the expression of this chemokine is usually associated with cancer promotion in terms of proliferation, migration, angiogenesis, and T-cell infiltration [43, 44]. A relationship between high levels of CXCL11 and PD-L1 expression was also reported [43].

The reduction of CXCL6 levels inhibited cervical cancer cell growth and metastatic potential, and reduced tumor growth in vivo [45].

The reduced transcription of chemokines highlighted by functional analysis of 43M2SD-DEGs may be ascribed to the abnormal E7 localization with consequent abrogation of its activity at the nuclear level. BioGRID analysis of indirect interactions of HPV E7 binders with 43M2SD-encoded proteins revealed essential proteins with signalling and/or transcriptional function for chemokines [46]. For example, E7 was found to bind to the TNF receptor-associated factors TRAF2, 3 and 5. These are known as signal transducers of the c-Jun N-terminal kinase (JNK) and the IκB kinase (IKK), leading to activation of the AP1 (Activator protein 1) and nuclear factor (NF)-κB transcription factors [47]. Furthermore, E7 was found to also bind to Jun and Rel, belonging to AP1 and NF-κB, respectively, suggesting that the reduced chemokines transcription in 43M2SD-expressing cells may depend not only on the interference with signal transduction mechanisms but also on the inhibition of the translocation of transcription factors into the nucleus.

Interestingly, CXCL9, CXCL10 and CXCL11 represent relevant hubs not only in the interactome of proteins encoded by DEGs affected by 43M2SD but also in the network obtained by the intersection of I7NUC and 43M2SD interactome, suggesting a role for these chemokines also in the mechanism of I7NUC. Most of the nodes identified in the I7NUC interaction network are transcription factors, all of which are downregulated except IFNG and ALB.

In contrast to many DEG-encoded proteins affected by I7NUC and identified as direct binders of E6, none of the DEG-encoded proteins modulated by the 43M2SD intrabody were binders of E7. However, the analysis of the interactome of E7-binders evidenced that also for this intrabody the same pathways involved in the activity of I7NUC could be involved, such as the Viral carcinogenesis, Cell cycle, Cellular senescence, Ubiquitin mediated proteolysis as well as HPV infection pathways, whose impairment is certainly compatible with the 43M2SD mechanism of action.

Quite expectedly, due to the complex relationships among HPV oncoproteins, we found that many proteins encoded by the anti-E6 I7NUC DEGs are also interactors of the HPV E1, E2, E5, and even E7 oncoproteins. This aspect might even offer a benefit for the therapeutic use of I7NUC and will be explored in further detail in future studies.

Although they have different cellular targeting and mechanisms, the two intrabodies exhibit some similarities, such as shared interaction hubs and possible pathways underlying their mechanisms, identified by both direct and indirect evidence. Interestingly, forty protein-coding genes were modulated by both intrabodies. Of these, twenty-four are present in CC, HNSCC and cervical cancer-derived cell lines according to the Human Protein Atlas. The remaining sixteen are genes that have not been correlated to HPV infections or HPV-related diseases. No specific articles linking these DEGs to HPV16 have been published so far, to the best of our knowledge.

Most shared DEGs encode membrane proteins (25/40, 62.25%), including eight chemoreceptors that are G-proteins, part of the large olfactory receptor gene family. Notably, twenty-five of the shared genes are regulated in opposite directions by the two intrabodies.

The opposite modulations should be considered whenever therapeutic intrabodies are used simultaneously, a possibility we previously suggested to explore the potential synergism of intrabodies, as discordant regulation could hinder any synergy. However, tests in animal models have shown that both intrabodies used individually show protection of mice by delaying the tumor onset in all animals, whereas 60% of them remained tumor-free for the entire duration of the experiment [10].

Confirmation of the expression of several genes identified as modulated by transcriptomic analysis was confirmed using RT-qPCR. However, this technique showed different sensitivity compared to transcriptomic analysis, and the significant results could only be obtained for some tested genes.

The transcriptomic analysis presented here confirmed at the gene expression level the involvement of the anti-E6 and-E7 intrabodies in subverting cellular pathways essential for tumor transformation, thus strengthening the possibility of their use as therapeutic drugs for HPV16-related diseases. However, the differences detected in the effects of the two intrabodies at the transcription level deserve further investigation, which could guide their use in clinical applications depending on the type of tumor or the grade of the lesion.

Conclusions

The transcriptome analysis of HPV16-positive SiHa cells expressing either the anti-16E6-I7NUC or the anti-16E7-43M2SD intrabodies shows that the anti-16E6 I7NUC intrabody, which contains a nuclear localisation signal, affects ten times more genes than the anti-16E7 43M2SD intrabody targeted in the ER compartment. The genes affected are involved in various cellular pathways essential for tumor transformation, and most of them are downregulated. By clarifying the mechanism of action of the two intrabodies, our data strengthen their antitumor potential as therapeutic agents for HPV16-related diseases. The differences observed in the effects of the two intrabodies at the transcriptional level could even prove advantageous for their clinical use, especially if employed for different types of tumors or grades of premalignant lesions.