HPV-associated penile cancer: Impact of copy number alterations in miRNA/mRNA interactions and potential druggable targets

Abstract

BACKGROUND:

Penile cancer (PeCa) is a rare disease, but its incidence has increased worldwide, mostly in HPV ${}^{+}$ patients. Nevertheless, there is still no targeted treatment for this carcinoma.

OBJECTIVE:

To predict the main signaling pathways involved in penile tumorigenesis and its potential drug targets.

METHODS:

Genome-wide copy number profiling was performed in 28 PeCa. Integration analysis of CNAs and miRNAs and mRNA targets was performed by DIANA-TarBase v.8. The potential impact of the miRNAs/target genes on biological pathways was assessed by DIANA-miRPath v.3.0. For each miRNA, KEGG pathways were generated based on the tarbase and microT-CDS algorithms. Pharmaco-miR was used to identify associations between miRNAs and their target genes to predict druggable targets.

RESULTS:

269 miRNAs and 2,395 genes were mapped in cytobands with CNAs. The comparison of the miRNAs mapped at these cytobands and the miRNAs that were predicted to regulate the genes also mapped in these regions, resulted in a set of common 35 miRNAs and 292 genes. Enrichment pathway revealed their involvement in five top signaling pathways. EGFR and COX2 were identified as potential druggable targets.

CONCLUSION:

Our data indicate the potential use of EGFR and COX2 inhibitors as a target treatment for PeCa patients.

Keywords

Genomic profiling molecular target biomarkers human papillomavirus penile tumor

1. Introduction

Penile cancer (PeCa) is a rare malignancy; however, it is a public health concern in developing countries. In Brazil, it represents 2% of all male cancers, occurring more frequently in the North and Northeast regions. Maranhão state (Northeast region), where the patients of this study were obtained, presents with the highest PeCa rates in Latin America and worldwide [1]. The patients in this region are usually diagnosed with advanced tumour stages, and lymph node involvement, which invariably confers poor prognosis. Human papillomavirus (HPV) has been linked to nearly 50% of PeCa cases worldwide [2, 3]. In our region however, our previous study, revealed a higher frequency of HPV-related PeCa, in nearly 97% of the cases [4].

The available surgical treatment for advanced PeCa involves partial or total amputations of the penis and lymphadenectomy, which are associated with post-surgical complications and high morbidities [5, 6]. The combined use of chemotherapy was observed to marginally increase the survival rate of the patients [7, 8, 9, 10], with however the inherent toxicity of the chemo-agents. In addition, the rarity of this tumour makes it difficult to conduct large-scale clinical trials, so there are still no effective treatment strategies for advanced PeCa patients. Therefore, effective and less toxic treatments, such as the ones based on targeted therapy, are critically needed to treat PeCa.

MicroRNAs (miRNAs), non-coding RNA molecules, have emerged as one of the promising molecular markers for therapeutic use in several types of cancers [11, 12, 13]. These molecules are preferentially mapped in genomic regions with copy number gains and/or losses, as well as in regions of fragile sites [14, 15]. Further studies have confirmed that miRNAs play important roles in oncogenesis by regulating the expression of genes located within the regions affected by CNA, and/or miRNAs expression changes can be driven by CNAs within their own genomic regions [16, 17, 18].

A large number of integration strategies for distinct “omics” platforms have been described and have significantly contributed to the characterization of tumour molecular landscapes, such as the ones performed by the International Cancer Genome Consortium (ICGC) and The Cancer Genome Atlas (TCGA). Most of these strategies are based on computational pipelines, where the data is gathered from independent and public avaliable resources. In PeCa there is relatively limited information of the somatic alterations that occur at the distinct molecular levels and consequently of the integrated molecular signatures.

Therefore, and considering the prevalence of PeCa in developing countries, as well as the increasing number of cases worldwide, our main objective was to determine the predictive impact of copy number alterations in the miRNA/mRNA interactions, and to identify the potentially affected functional pathways in penile tumorigenesis. In addition, we aimed to determine whether the observed alterations and/or the signaling pathways involved, could be potentially targeted by the current clinically approved drugs.

2. Materials and methods

2.1 Experimental subjects

2.1.1 Patient accrual, social and clinical-histopathological variables

Tumour tissues of 28 PeCa patients, collected at the Hospital Aldenora Bello (São Luís, Maranhão, Brazil) were included in this study. The fresh tumour tissue samples were obtained at the diagnosis and before any cancer treatment. The clinical and histopathological variables were obtained from patients’ medical records. All patients signed consent form, which was approved by the Ethics Committee on Human Beings at the Federal University of Maranhão and by the National Research Ethics Commission (CONEP), Brazil (CAAE protocol number: 46371515.5.0000.5087).

The mean age at diagnosis of the patients was 64.6 $\pm$ 15.2, ranging from 32 to 85 years old. Most of the patients were married (75%) and presented low level of education (only 17.9% presented high school degree); 50% self-declared current smoker or declared to have been a smoker in the past 5 years, and 63.3% declared to regularly drink alcohol.

The tumours were classified as Squamous Cell Carcinoma (SCC), presenting an average size of 3.6 $\pm$ 2.2 cm, ranging from 0.5 to 10 cm. In most of the patients (39.3%), the tumours were restricted to the glans. The most frequent histological subtypes were condylomatous and usual, present in 42.9% and 39.3% of the cases, respectively. Tumour grades II and III were the most frequent, present in 59.6% and 25.9% of the patients, respectively. Lymph node and perineural invasion were observed at diagnosis in 16% and 20% of the patients, respectively. Phimosis was observed in 54.5% of the patients, most of which (91.7%) with over 40 years old. Most of the patients (75%) were submitted to partial penectomy. This clinical-histopathological information from the patients is presented in Table 1.

Table 1
Clinical-histopathological characterization of the PeCa patients evaluated in this study ( $n=$ 28)

Case	Age	Phymosis ocurrence	HPV subtype	Histological subtype	Predominant lesion	Tu. size ${}^{\rm b}$	Tu. grade	Tu. stage	Lymphatic invasion	Perineural invasion
01	62	No	16	Usual	Ulcerated	2.5	II	pT3	Yes	Yes
02	80	N/A	16	Condilomatous	Ulcerated	10	II	pT2	No	No
03	81	Yes	18,53	Condilomatous	Ulcerated	2.1	II	pT1	No	No
04	76	Yes	30	Condilomatous	Ulcerated	3	I	pT1b	Yes	No
05	74	Yes	73,74	Usual and basaloid	Verrucous	0.8	I	pT2	No	No
06	78	N/A	16,06	Usual and basaloid	Vegetative	4	III	pT1b	No	No
07	57	No	58	Usual and condilomatous	Nodular ulcerated	3	III	pT2	No	No
08	85	No	16	Condilomatous	Nodular ulcerated	4	III	pT3	No	Yes
09	83	No	16	Usual	Ulcerated	4.5	III	pT3	No	Yes
10	40	No	59,66	Condilomatous	Ulcer-vegetative	3	II	pT2	No	No
11	68	N/A	16,66	Condilomatous	Verrucous	5.5	I	pT3	No	No
12	41	Yes	56	Condilomatous	Vegetative	3.5	II	pT2	No	No
13	74	No	11,35,59	Condilomatous	Vegetative	9	II	pT3	No	No
14	65	N/A	16,74	Condilomatous	Vegetative	3.5	II	pT2	No	No
15	51	N/A	16,35,59	Usual	Ulcer-vegetative	3	II	pT2	Yes	No
16	32	No	16,44,74	Usual and condilomatous	Ulcer-vegetative	3.7	II	pT2	No	No
17	69	No	16,30	Condilomatous	Vegetative	3.5	I	pT2	No	No
18	37	Yes	Positive ${}^{\rm a}$	Usual	Nodular ulcerated	5.0	III	pT3	No	Yes
19	54	Yes	Positive ${}^{\rm a}$	Usual	Vegetative	0.8	II	pT2	Yes	No
20	44	Yes	16	Usual	Ulcerated	1.7	II	pT1	No	No
21	51	Yes	16	Basaloid	Vegetative	N/A	II	pT2	N/A	N/A
22	71	Yes	Positive ${}^{\rm a}$	Usual	Ulcerated	2	III	Tis	No	No
23	74	No	Positive ${}^{\rm a}$	Condilomatous	Ulcerated	1.8	II	pT1a	No	No
24	83	Yes	35,59	Usual	Ulcer-vegetative	N/A	II	pT2	N/A	N/A
25	77	N/A	Positive ${}^{\rm a}$	Usual	Ulcerated	0.5	ni	N/A	No	No
26	73	Yes	16	Usual	Ulcerated	N/A	II	pT1	N/A	N/A
27	71	No	16	Usual	Ulcerated	3.2	III	pT3	No	No
28	59	Yes	Positive ${}^{\rm a}$	Condilomatous	Verrucous	5.5	II	pT3	No	Yes

${}^{\rm a}$ No information on HPV genotype; ${}^{\rm b}$ Tu $=$ tumor (mm); N/A $=$ not available.

2.1.2 HPV genotyping

HPV genotyping was performed for all the patients by nested-PCR (according to BioGenetics Molecular Technologies (Uberlândia, Minas Gerais, Brazil) (patent number BR102017004615.0), and DNA sequencing, according to our previous protocol [4]. All the patients were HPV positive, of which 95.5% presented high-risk subtypes. HPV16 was the most common subtype, present in 63.6% of cases, followed by subtypes 35, 59, 74 (13.6%, each); 30, 44, 66 (9%, each); 06, 18, 53, 58 and 73 (4.5%, each). Multiple infections were detected in 50% of the cases (Table 1).

2.1.3 Copy number data of the PeCa tumours

All the patients were profiled for genome-wide copy number alterations (CNAs) by array-Comparative Genomic Hybridization (aCGH) analysis, as reported by Macedo et al. [4]. The aCGH analysis showed a total number of 2,314 CNAs, with a mean of 82.6 $\pm$ 53.5 CNAs per case (ranged from 23 to 224). Most of the CNAs observed were gains (56.9%/1,316), followed by losses (37.3%/864), amplifications (4.0%/92), and deletions (1.8%/42).

2.2 Integration analysis of CNA and miRNAs and mRNA targets

This analysis was performed by two different approaches: i. By the mapping of the genes and miRNAs at the main regions with CNA, called here CNA-genes and CNA-miRNAs, respectively, and ii. By the identification of common gene targets of the selected miRNAs above, that may be affected by both CNAs and miRNA expression alterations. The workflow of the study design and corresponding results are presented in the results section (Fig. 1).

Figure 1.

Study design and corresponding results of the aCGH and miRNA data integration approaches, followed by the functional enrichment pathway and druggable targets predicted analysis. In one approach (269 CNA-miRNAs), 1,355 interactions were identified involving 65 CNA-miRNAs and 898 target genes (TarBase v.8.0). MiRPath v.3.0 was used for the enrichment of pathways considering two prediction algorithms, tarbase and microT-CDS. Five common pathways (14 CNA-miRNAs and 25 target genes) were identified comparing the two pathway methods. In addition, tarbase enrichment revealed an HPV-related carcinogenesis pathway (7 CNA-miRNAs and 13 target genes). In a second approach (2,395 CNA-genes), 1,358 interactions were identified involving 348 miRNAs and 292 CNA-genes. Comparing the interactions with the first approach, 115 common interactions were revealed, that is, between 35 CNA-miRNAs and 78 CNA-genes. Pathways enrichment was conducted with these interactions where it revealed 13 top pathways (10 CNA-miRNA and 13 CNA-genes). The Pharmarco-miR program was used to identify drugs targeting miRNA targets in both approaches.

2.2.1 Mapping of miRNAs in genomic regions with CNAs

Initially, the complete aCGH data for all the 28 patients were analyzed to identify the most affected cytobands with CNAs. Only the CNAs present in at least 25% of the cases were considered. Next, the annotation of genes and miRNAs located in the cytobands with CNAs were obtained from the Agilent Cytogenomics v.7.0 interval base reports (Agilent Technologies, Inc.). The location of each miRNA was further confirmed using miRBase v.22 (http://www.mirbase.org) [19]. The genomic regions that were sites of HPV integration were accessed by the HPVbase platform (http://crdd. osdd.net/servers/hpvbase/index.html), according to Kumar Guptar and Kumar [20].

2.2.2 Association of genomic regions with CNAs and miRNAs with clinical and histopathological parameters from the patients

The association analysis of each clinical-histopathological variable and copy number alterations was performed using bivariate tests and $X^{2}$ Pearson or Fisher’s Exact test, when applicable. The significance level was $P<$ 0.05. The statistical tests were performed using the Graphpad Prism v.8.4.3 software.

2.2.3 Identification of mRNA targets

The targets of the miRNAs mapped in the regions with CNA (CNA-miRNAs) were identified by the DIANA-TarBase v.8 (carolina.imis.athena-innovation. gr/diana_tools/web/index.php?r=tarbasev8/index) [21]. Only target genes presenting miRNA Target Gene (miTG) score $\geqslant$ 0.9 were included in this analysis. In addition, we identified the miRNAs that potentially regulate the CNA-genes and selected the common miRNAs/target genes of both integration analysis.

2.2.4 Enrichment pathways analysis

To assess the potential impact of the miRNAs/target genes (that were mapped in the main cytobands selected above) on biological pathways associated with cancer, we used DIANA-miRPath v.3.0 tools (http://snf-515788.vm.okeanos.grnet.gr/) based on $p$ -values (FDR correction), according to Vlachos et al. [22]. The KEGG (Kyoto Encyclopedia of Genes and Genomes) was used for pathways enrichment analysis ( $p<$ 0.05). For each miRNA, KEGG pathways were generated based on the tarbase and microT-CDS algorithms of the DIANA tools. Top pathways were identified using the “pathway union” function.

Table 2
Cytobands with CNAs in PeCa ( $n=$ 28) analyzed by aCGH in association with clinical-histopathological parameters

Cytoband	CNA event	HPV integration site (genotype) ${}^{\rm a}$	Variables	$P$ -value
2p16.3	Gain	Yes (16, 18)	Tumor pT1	0.0130
			Ulcerated	0.0338
			Tumor size $>$ 5.1 cm	0.0238
2p12-p11.2	Gain	Yes (16)	Subtype usual	0.0472
			Vegetative	0.0461
			Tumor size $>$ 5.1 cm	0.0204
3q26.1	Gain	Yes (45)	Ulcer-vegetative	0.0126
			Verrucous	0.0343
4p14-p13	Loss	Yes (16)	Tumor size $>$ 5.1 cm	0.0198
4q13.2	Loss/deletion	Yes (16)	Ulcer-vegetative	0.0017
7p22.3-p11	Gain	Yes (16, 45)	Perineural invasion	0.0075
7q21.11	Gain	Yes (16, 35)	Tumor pT3	0.0381
			Tumor size $>$ 5.1 cm	0.0409
8p23.1	loss/deletion	Yes (16)	Vegetative	0.0233
9p21	Gain	Yes (16,18)	Tumor pT1	0.0078
11q24-q25	Gain	Yes (16)	Tumor pT1	0.0243
14q11.2	Loss	Yes (16, 18)	Perineural invasion	0.0038
14q12	Gain	Yes (16)	HPV16 and multiple infections	0.0044
15q11.2-q13.3	Gain	Yes (18)	Vegetative	0.0490
			Tumor pT1	0.0297
15q26.2-q26.3	Gain	Yes (16)	Age $>$ 40 years	0.0330
			Grade I	0.0022
22q11.21	Gain	Yes (33)	Perineural invasion	0.0412
			Tumor pT2	0.0404

${}^{\rm a}$ genotypes identified in previous studies according to HPVbase. aCGH performed by Cytogenomics v.5.0 software; Agilent Technologies Inc, CA, algorithm ADM-2; threshold $>$ 6.

2.3 Analysis of druggable targets

Pharmaco-miR (www.pharmaco-mir.org) was used to identify associations between miRNAs and their target genes to predict druggable targets, as described in the literature [23]. This analysis was performed for: i. mRNA targets identified by both algorithms (tarbase and microT-CDS); ii. mRNA targets related to HPV infection, and iii. miRNAs/mRNA target genes involved in the the top KEGG pathways identified. PathVisio v.3.3 [24] and Biorender (https://biorender.com/) were used to rebuild KEGG pathways including the miRNAs, their respective gene targets and corresponding predicted drugs.

3. Results

The overall results of this study with the corresponding study design are presented in Fig. 1.

3.1 Copy number alterations and their association with clinical and histopathological parameters

Genome-wide copy number analysis in the 28 PeCa patients analyzed by aCGH revealed 38 altered cytobands with copy number alterations (CNAs), from which 89.5% were HPV integration sites in the host genome. The association analysis of the main affected cytobands with CNAs and the clinical-histopathological parameters from the patients revealed 15 cytobands (Table 2). All these altered cytobands were already described as HPV integration sites according to HPVbase. Gains on 14q12 presented association with HPV16 and multiple infections, while gains on 15q26.2-q26.3 were associated with age $>$ 40 years. Gains on 2p16.3, 9p21, 11q24-q25 and 15q11.2-q13.3 showed association with early stages of the primary tumour (pT1). On the other hand, gains on 7q21.11, 7p22.3-p11 and 22q11.21, and loss on 14q11.2 and 4p14-13 were associated with intermediate stages of the tumour (perineural invasion, pT2, pT3 and tumour size $>$ 5.1 cm).

3.2 Mapping of genes and miRNAs in genomic regions with CNA

According to the annotations of the CytoGenomics base interval report, 30 of the 38 affected cytobands harbored miRNAs (Fig. 1, Supplementary Table 1). Interestingly, 28 of these 30 cytobands (93.3%) are HPV integrated sites according to HPVbase. A total of 269 miRNAs were mapped in these cytobands (CNA-miRNAs), 211 in regions affected by copy number gains (78.4%), 27 by losses (10%), 13 by amplification (4.8%), and 18 by losses and deletions (6.7%).

A second integration approach was based on the identification of miRNAs that potentially regulated the genes mapped at the main affected cytobands: 2,395 genes were identified to be mapped at the 38 cytobands with CNAs (CNA-genes). Based on the TarBase/Diana Tools, considering miTG $\geqslant$ 0.9, 292 of these genes were predicted to be regulated by 348 miRNAs. The comparison of the miRNAs mapped at the main cytobands and the miRNAs that were predicted to regulate the genes also mapped in these regions, resulted in 35 miRNAs (Fig. 1, Supplementary Table 2A).

3.3 Prediction of mRNA targets in regions with CNA

Target prediction analysis (considering miTG $\geqslant$ 0.9) revealed that 65 out of 269 miRNAs mapped in the affected cytobands, potentially regulate 898 gene targets (Fig. 1, Supplementary Table 2B). MiR-30d, mapped at 8q24.22 presented the highest number of targets (229 genes), followed by miR-30b (mapped at 8q24.22), and miR-548d (mapped at 8q24.13) with 224 and 134 target genes, respectively. MiR-151a-3p (8q24.22), miR-649 (22q11.21), miR-1269a (22q11.21), miR-3150b-3p (8q22.1), miR-4286 (8p23.1), miR-4435 (2p11.2), miR-4473 (9p21), and miR-4802-5p (4p14) presented the lowest number of targets (1 gene, each).

In addition, the identification of the genes mapped at the affected cytobands and that were predicted to be regulated by the 35 miRNAs mapped at the same cytobands, resulted in 78 mRNA targets (Fig. 1). From these miRNAs, miR-30b-3p, miR-30d-5p and miR-548d-3p were the ones to present the highest number of target genes, 19, 19 and 10, respectively.

3.4 Pathway enrichment analysis

The main biological functions related to the 65 miRNAs mapped at the 30 cytobands (CNA-miRNAs) were evaluated by KEGG pathway enrichment analysis using two distinct algorithms. The tarbase algorithm revealed 66 pathways, among the top 13: adherens junction (hsa04520), cell cycle (hsa04110), ECM-receptor interaction (hsa04512), fatty acid biosynthesis (hsa00061), fatty acid metabolism (hsa01212), glioma (hsa05214), hippo signaling pathway (hsa04390), lysine degradation (hsa00310), mucin type O-Glycan biosynthesis (hsa00512), prion diseases (hsa05020), proteoglycans in cancer (hsa05205), steroid biosynthesis (hsa00100), and viral carcinogenesis (hsa05203). The microT-CDS prediction algorithm revealed 51 main KEGG pathways, among the top 13: axon guidance (hsa04360), foxO signaling pathway (hsa04068), GABAergic synapse (hsa04727), hippo signaling pathway (hsa04390), long-term depression (hsa04730), lysine degradation (hsa00310), morphine addiction (hsa05032), mucin type O-Glycan biosynthesis (hsa00512), prion diseases (hsa05020), proteoglycans in cancer (hsa05205), signaling pathways regulating pluripotency of stem cells (hsa04550), steroid hormone biosynthesis (hsa00140), and TGF-beta signaling pathway (hsa04350). The significant pathways revealed by both of the miRNA prediction algorithms used are shown in Fig. 2 and Supplementary Table 3A and B.

Figure 2.

Targeted KEGG pathway clusters/heatmap regulated by the miRNAs mapped at the cytobands mostly affected by CNAs in the cases studied. A: Top pathways identified by the tarbase algorithm; B: Top pathways identified by the microT-CDS algorithm ( $P<$ 0.05; DIANA/miRPath v.3). Red colors indicate a stronger role of the miRNAs on a given pathways as compared to the lighter colors. Blue arrows indicate the miRNAs and corresponding signaling pathways that were identified by both algorithms (tarbase and microT-CDS).

Table 3

Target genes identified as potentially regulated by miRNAs mapped in regions with gains or amplifications and also described as negatively controlled by HPV oncoproteins (DIANA-miRPath v.3.0; KEGG: hsa05203)

Gene	HPV oncoprotein	CNA-miRNA	$P$ -value	Biological consequence
ATP6V0D1	E5	miR-148a-3p	0.0001	Inhibition of V-ATPase-dependent endosomal acidification and increased recycling of EGFR-receptor to the cell surface
CCNA2	E7	miR-148a-3p
miR-4429
miR-548a-3p
miR-548d-5p	0.0001
0.0048
0.0152
0.0262	Deregulation of the cell cycle by inactivating cyclin A2 (CCNA2)
CDK2	E7	miR-548a-3p	0.0152	Alteration of cell cycle
CDKN1B	E7	miR-30d-3p
miR-148a-3p
miR-548a-3p
miR-550a-3-5p	0.0007
0.0001
0.0152
0.0499	Disruption of the cell cycle checkpoint by inactivating P27 (CDKN1B)
CHD4	E7	miR-148a-3p
miR-185-5p	0.0001
0.023	Increased cell growth
EP300	E6	miR-30d-3p
miR-148a-3p	0.0007
0.0001	Inhibition of P53 activation through the E6-P300 complex (EP300)
IRF3	E6	miR-148a-3p	0.0001	Inactivation of apoptotic activity promoted by P53 by the association E6-IRF3 (IRF3)
JUN	E7	miR-4429	0.0048	Disturbance of the transcription factor c-JUN (JUN)
PKM	E7	miR-148a-3p
miR-4429	0.0001
0.0048	Increased cell growth
RB1	E7	miR-30d-3p
miR-148a-3p	0.0007
0.0001	Inactivation of tumor suppressor protein pRB
RBL2	E7	miR-550a-3-5p	0.0499	Increased cell proliferation by inactivating P130 (RBL2)
TBP	E7	miR-148a-3p
miR-548a-3p	0.0001
0.0152	Inhibition of transcription factors
UBR4	E7	miR-148a-3p
miR-4429	0.0001
0.0048	Alteration of cell membrane morphogenesis and cell migration by the interaction of E7 with P600 (UBR4)

DIANA-miRPath v.3.0; KEGG: HPV-mediated viral carcinogenesis pathway – hsa05203.

Notably, five common pathways were identified by both algorithms: hippo signaling pathway (miR-548d-5p, miR-548k, miR-548n and miR-4517), lysine degradation (miR-339-5p, miR-548d-3p and miR-624-3p), mucin type O-Glycan biosynthesis (miR-30b-5p, miR-30d-5p and miR-5680), prion diseases (miR-30d-3p and miR-148a-3p) and proteoglycans in cancer (miR-100-5p, miR-148a-3p and miR-217). Furthermore, considering that all the tumours were HPV positive, we also investigated the HPV-associated carcinogenesis pathway, in which 13 genes were identified as potentially regulated by miRNAs mapped in regions with gains or amplifications (tarbase algorithm). Interestingly, these genes have also been described as negatively controlled by HPV oncoproteins (KEGG- hsa05203) (Supplementary Fig. 1, Table 3).

Pathway enrichment analysis using the two algorithms (tarbase and microT-CDS) was also conducted considering the common genes ( $n=$ 78) mapped at the main cytobands and that were also target of the miRNAs ( $n=$ 35) mapped at the same affected cytobands. In this analysis however, no common pathways were identified between these two prediction algorithms. Therefore, the pathways enrichment analysis results of the tarbase algorithm were chosen, considering that they were based on experimental validations of miRNA/mRNA target interactions. This analysis revealed 63 pathways, which top 13 involved 13 genes (CCNE2, E2F5, EPAS1, INHBA, LYN, PDE7A, PPP1CB, PPP3R1, RAC1, SEMA3D, SOS1, YWHAQ, and YWHAZ) (Table 4). Of note, these genes are also involved in molecular pathways directly related to cancer, such as the P53 (hsa04115) and TGF- $\beta$ (hsa04350) signaling pathways, proteoglycans in cancer (hsa05205), and viral carcinogenesis (hsa05203).

Table 4

Top 13 signaling pathways with the involvement of 13 genes observed with gains of copy number and also targeted by miRNAs mapped at the main affected cytobands (Tarbase algorithm)

CNA-gene	Cytoband	CNA-miRNA	Top pathways
CCNE2	8q11.1-24.3	miR-30b-5p
miR-30d-5p	Oocyte meiosis, pathways in cancer, P53 signaling pathway
E2F5	8q11.1-24.3	miR-1-3p	Cell cycle
EPAS1	2p25.3-p11.1	miR-148a-3p	Pathways in cancer
INHBA	7p22.3-p11	miR-548d-3p	Signaling pathways regulating pluripotency of stem cells and TGF-beta signaling pathway
LYN	8q11.1-24.3	miR-185-5p	Viral carcinogenesis
PDE7A	8q11.1-24.3	miR-30b-5p
miR-30d-5p	Fatty acid biosynthesis and metabolism and pathways in cancer
PPP1CB	2p25.3-p11.1	miR-148a-3p	Hippo signaling pathway, oocyte meiosis and proteoglycans in cancer
PPP3R1	2p25.3-p11.1	miR-298
miR-30b-5p
miR-30d-5p	Axon guidance and oocyte meiosis
RAC1	7p22.3-p11	miR-589-3p	Adherens junction, axon guidance, pathways in cancer, proteoglycans in cancer and viral carcinogenesis
SEMA3D	7q21.11	miR-548a-3p	Viral carcinogenesis
SOS1	2p25.3-p11.1	miR-148a-3p
miR-548d-5p	FoxO signaling pathway, gap junction, glioma, pathways in cancer and proteoglycans
in cancer
YWHAQ	2p25.3-p11.1	miR-1-3p	Cell cycle, hippo signaling pathway, oocyte meiosis and viral carcinogenesis
YWHAZ	8q11.1-24.3	miR-1-3p	Cell cycle, hippo signaling pathway, oocyte meiosis and viral carcinogenesis

Figure 3.

Integrative network analysis showing the main pathways and the druggable targets predicted by Pharmaco-miR analysis. The network was generated by PathVisio v3.3 and Biorender softwares using the data generated through the Diana-Tools, KEGG and Pharmaco-miR.

3.5 Drug targets prediction

Finally, drug prediction analysis was performed using the five common pathways identified from CNA-miRNA (hippo signaling, lysine degradation, O-glycan mucin biosynthesis, prion, and proteoglycan diseases in cancer), as well as using the 13 main common pathways identified from the CNA-genes, and CNA-miRNAs analysis. Moreover, considering the high prevalence of HPV infection in the PeCa patients of this study, the HPV-mediated carcinogenesis pathway was also used for drug prediction. This analysis revealed drugs currently used to treat PeCa, such as Cisplatin, which acts on the PARD6B (regulated by miR-548k/548d-5p), and RB1 (regulated by miR-30d-3p and miR-148a-3p) genes. In addition, the MTOR (regulated by miR-100-5p) and SOS1 (regulated by miR-148a-3 and miR548k/548d-5p) pathways were found to be susceptible candidates to pharmacological inhibition by Doxorubicin and Imatinib, respectively. Interestingly, Cetuximab (an EGFR inhibitor) and Celecoxib (selective cyclooxygenase-2 inhibitor (COX2) were detected as alternative therapies for PeCa (Fig. 3 and Supplementary Table 4).

4. Discussion

PeCa is a rare disease in developed countries and, despite the recent increase in the number of molecular studies performed in these tumours [25, 26, 27, 28, 29], their molecular basis is still not well understood. In this study we evaluated the role of miRNAs in the molecular pathways associated with HPV-positive tumours, considering the high frequency (63.3%) of high-risk HPV genotypes (HPV16) in the PeCa patients. This is an unusual high frequency compared to the literature, which have described the HPV16 genotype as the most prevalent in PeCa, although in frequencies less than 50% [3, 27].

The aCGH analysis revealed 89.5% of the altered cytobands are mapped in HPV integration sites in the host genome. These findings support that HPV integration sites are generally located on or nearby regions of fragile sites, and that CNAs in these regions can be potentially triggered by the virus insertion [30]. It is also relevant to point out, that all the 15 cytobands that showed an association with clinical and histopathological parameters are regions of HPV integration sites. For example, the 14q12 region showed an association with the HPV16 genotype and multiple infections. CNAs on chromosomes 2p, 9p, 11q and 15q were associated to early stages of the tumour, while CNAs on 4p, 7p, 7q, 14q and 22q to intermediate tumour stages. These findings, together with our previous study [4], reinforce that CNAs play an important role in HPV-related PeCa, and can indicate candidate genes within these affected regions, that can be used as molecular markers for these patients.

Herein, we identified potential molecular markers affected by CNAs using prediction analysis to determine the interaction of genes and miRNAs mapped at these regions with CNAs and their corresponding gene targets. To our knowledge this is the first study conducted in PeCa patients, where this integration approach was performed. This integrative analysis, followed by functional enrichment pathways analysis, using both tarbase and microT-CDs algorithms, revealed five common pathways: Hippo signaling, lysine degradation, mucin type O-Glycan biosynthesis, prion diseases, and proteoglycans in cancer. Hippo signaling, which is a pathway associated with cell proliferation and survival [31]. The miRNAs that were mapped in the CNAs, miR-548d-5p, miR-548k, miR-548n, and miR-4517, were previously demonstrated to regulate several genes involved in these cellular process [32, 33, 34, 35]. Interestingly, this pathway presents a feedback relationship with the TGF-B signaling pathway [36], in which the downstream genes BMPR2 and SMAD4 are involved. These genes are also targets of miR-548d-5p (BMPR2) and miR-4517 (BMPR2 and SMAD4), which loss was previously associated with aggressive tumour phenotypes, such as tumours with larger sizes, metastasis development, and chemotherapy resistance [37, 38]. Therefore, we suggest that the negative regulation of these genes by miR-548d-5p, and miR-4517 can be one of the molecular mechanisms relevant to the penile tumorigenesis.

The role of miRNAs in the lysine’s degradation pathway is not completely understood. MiR-548d-3p, mapped at 8q24.13 and affected in our cases by gain of copy number, is predicted to regulate several gene targets involved in this pathway, such as KMT2 A, KMT2D, SUV420H1, WHSC1L1, and SETD7. The altered expression of these genes has been described in several cancers [39, 40, 41, 42, 43], however the mechanisms that lead to their expression alterations is not clear.

The third pathway affected by the miRNAs located in regions with CNAs of this study, is the biosynthesis of mucin’s type O-glycan. O-glycans-type mucins are constituents of the main compounds of mucous, forming physical-chemical barriers [44]. The O-glycosylation process involve several genes from the GALNT or GalNac-T (N-acetyl-galactosaminyltransferases) family, which abnormal expression have been associated with cancer progression [44, 45]. In our study, the GALNT1, GALNT2, GALNT3 and GANLT7 genes were identified as potential targets for miR-30b-5p and miR-30d-5p, mapped at 8q24.22 region, affected by gain of copy number. Interestingly, Gaziel-Sovran et al. [46] observed an increase in the metastatic potential in melanoma by the repression of GALNT1 and GALNT7 due to the high expression of miR-30b/30d. In gastric carcinoma, negative regulation of GALNT2 is associated with more advanced stages of the tumour, lymph node metastasis and reduced disease-free survival [45]. It was also observed that the downregulation of GALNT2 modifies the O-glycosylation of EGFR and affects the EGFR-AKT signaling [44]. Additionally, the lower expression of GALNT3 in pancreatic cancer has been associated with altered glycosylation of the ErbB family receptors, and with increased aggressiveness [47]. Altogether, these findings support the role of miR-30b/30d in the O-glycan mucin pathway via regulation of the GALNT gene family, and their association with aggressive tumour types.

Our analysis also revealed that miR-20d-3d and miR-148a-3p target the PRNP gene. This gene is involved in Prion diseases, where it regulates pathways of oxidative stress, proliferation, adhesion, and cell survival. Although most of the studies on PRNP are focused on the nervous system, recent evidence indicates that this gene also regulates these process in cancer cells [48]. The role of miRNAs in the regulation of PRNP is still uncertain, since the physiological functions of the PrP ${}^{\rm c}$ protein, encoded by the PRNP gene, remain ambiguous [49]. However, it is suggested that changes in PrP ${}^{\rm c}$ expression affect NF- $\kappa$ B activation, TNF $\alpha$ production [49], and EGFR signaling [50], pathways and mechanisms closely related to cancer development. Moreover, our integrated analysis showed that miR-100-5p, miR-148a-3p, and miR-217 negatively regulate genes involved in the Proteoglycans pathway. Proteoglycans are important constituents of the extracellular matrix and regulate the bioavailability of hormones, growth factors, receptors, cytokines, and gene expression [51].

Considering that in this study all the tumours were HPV positive, we also investigated the HPV-associated carcinogenesis pathway. The tarbase algorithm revealed that both miRNAs and HPV oncoproteins can negatively regulate genes related to proliferation, cell cycle and apoptosis. For example, although endosomal acidification is required for the initial stages of HPV infection [52], inhibition of acidification by the HPV oncoprotein E5 has already been reported [53, 54, 55]. Moreover, the ATP6V0D1 gene has been identified as a target for miR-148a-3p. This gene encodes one of the subunits of the vacuolar-ATPase enzyme and is responsible for endosomal acidification and consequent degradation of some cellular components, such as the EGFR [56]. Thus, miR-148a-3p and HPV-E5 may both inhibit the endosomal acidification and, consequently, increase the availability of EGFR. The involvement of EGFR in PeCa was previously reported in other studies, including our own [4], where EGFR was observed with gain of copy number and gene overexpression.

Other genes that are usually targeted by HPV oncoproteins have also been identified as targeted by miRNAs in the HPV-mediated carcinogenesis pathway. In the cases of this study, miR-148a-3p and miR-30d-3p are mapped at regions with gains (at 7p15.2 and 8q24.22, respectively), which were among the miRNAs predicted to regulate RB1 and CDKN1B, both genes involved in cell cycle and apoptosis [57]. Previously, we demonstrated the reduced expression of RB1 in PeCa and suggested the occurrence of other mechanisms for RB repression than the canonical HPV/TP53/RB1 signaling pathway [4]. Therefore, based on this present study, upregulation of these miRNAs’ expression, could be another mechanism that leads to the negative regulation of RB1 in HPV ${}^{+}$ PeCa.

Additional integrated analysis considering the 78 genes affected by CNAs and regulated by CNA-miRNAs highlighted the involvement of other genes (CCNE2, E2F5, EPAS1, INHBA, LYN, PDE7A, PPP1CB, PPP3R1, RAC1, SEMA3D, SOS1, YWHAQ, and YWHAZ) in the previously discussed pathways, such as P53 signaling, hippo and TGF- $\beta$ signaling, proteoglycans in cancer and viral carcinogenesis. Thus, we emphasize the critical role of these pathways in penile tumorigenesis, however, functional assays are essential to confirm their biological role in PeCa.

The conventional treatment for PeCa patients is usually based on 5-Fluorouracil, Doxorubicin, Cisplatin, Ifosfamide, Methotrexate, and Paclitaxel, mostly in combination with radiotherapy. This treatment regimen is effective, however, patients’ tumours eventually become resistant, which can lead to disease recurrence and death [6, 9].

The drug prediction analysis conducted in our study, based on the main pathways described above indicated several molecular markers which inhibitory drugs could be repurposed for the treatment of HPV-related PeCa. The EGFR gene emerged as one of these known druggable targets, which we consider one of the most promising markers that can be used for PeCa treatment. This is based on our findings showing its gain of copy number and higher expression in most of the PeCa cases analyzed [4], and its involvement in the HPV-mediated carcinogenesis pathway, where its expression is affected by the endosomal acidification triggered by miR-148a-3p and the HPV-E5 oncoprotein. Thus, Cetuximab, a recombinant monoclonal antibody directed to EGFR, could be considered an alternative therapy for patients with PeCa that express this receptor. Interestingly, a recent study [58] has demonstrated that the growth of one PeCa-derived cell line overexpressing EGFR, was inhibited by EGFR inhibitors, corroborating with our findings. In addition, combined with chemotherapy, Cetuximab is the first line therapy for metastatic colorectal and squamous cell carcinomas [59, 60, 61], where it has been shown satisfactory responses [62, 63].

In addition to EGFR inhibitors, the drug prediction analysis performed revealed the Celecoxib, a selective inhibitor of COX2. Interestingly, we previously reported in the PeCa cases of this study, the association of the overexpression of EGFR and COX2 [4], suggesting that their inhibitors could be used for PeCa treatment. Furthermore, in vitro studies have demonstrated the efficacy of the combination of COX2 inhibitors with other drugs in cancer cells. For example, Mishra et al. [64] demonstrated that N-acetyl-cysteine and Celecoxib inhibit the initiation of cutaneous tumours. Velmurugan et al. [65] showed that the combined treatment of Celecoxib and Calyculin-A inhibited the epithelial-mesenchymal transition in oral cancer cells. Other studies showed that treatment with Celecoxib prevents Doxorubicin-induced drug resistance in canine and mouse lymphoma cell lines [66]. These findings reinforce that COX2 and EGFR are promising druggable targets and could be considered as alternative therapeutic strategies to treat PeCa patients.

5. Conclusion

In conclusion, our analysis showed that CNAs present a significant impact on the regulation of miRNAs and their corresponding gene targets. The genes, miRNAs, and mRNA targets mapped in chromosome regions with CNAs, can affect signaling pathways with relevance to the HPV-mediated carcinogenesis, and point out to druggable targets that could be potentially “repurposed” to treat PeCa patients with HPV infection.

Author contributions

Conception: J.S. and S.P.

Interpretation or analysis of data: J.S., L.C. and S.P.

Preparation of the manuscript: J.M., L.C. and S.P.

Revision for important intellectual content: J.M., L.N., R.C. A.D., A.K., R.M., A.P.S., L.C. and S.P.

Supplementary data

The supplementary files are available to download from http://dx.doi.org/10.3233/CBM-210035.

sj-pdf-1-cbm-10.3233_CBM-210035.pdf - Supplemental material

Supplemental material, sj-pdf-1-cbm-10.3233_CBM-210035.pdf

Footnotes

Acknowledgments

The authors would like to thank the personnel from Tissue Culture and Cytogenetics Laboratory of Institute of Evandro Chagas (Belém, PA, Brazil) for performing the scanning of the array slides. The authors thank Coordenação de Aperfeiçoamento Pessoal de Nível Superior (CAPES) for providing scholarship for J.S (Master Postgraduate Program in Health Science - Federal University of Maranhão, Brazil).

Conflict of interest

The authors declare that they have no competing interest.

References

Coelho

R.W.P.

et al., Penile cancer in Maranho, Northeast Brazil: The highest incidence globally? BMC Urology 18 (2018), 1–7. doi: 10.1186/s12894-018-0365-0.

Douglawi

and Masterson

T.A.

, Penile cancer epidemiology and risk factors: A contemporary review, Current Opinion in Urology 29 (2019), 145–149. doi: 10.1097/MOU.0000000000000581.

Sakamoto

et al., Etiological role of human papillomavirus infection in the development of penile cancer, International Journal of Infectious Diseases 78 (2019), 148–154. doi: 10.1016/j.ijid.2018.11.003.

Macedo

et al., Genomic profiling reveals the pivotal role of hrHPV driving copy number and gene expression alterations, including mRNA downregulation of TP53 and RB1 in penile cancer, Molecular Carcinogenesis 59 (2020), 604–617. doi: 10.1002/mc.23185.

O’Brien

J.S.

et al., Penile cancer: Contemporary lymph node management, Journal of Urology 197 (2017), 1387–1395. doi: 10.1016/j.juro.2017.01.059.

Necchi

et al., Clinical outcomes of patients with locally advanced penile squamous-cell carcinoma: Results of a multicenter analysis, Clinical Genitourinary Cancer 15 (2017), 548–555.e3. doi: 10.1016/j.clgc.2017.02.002.

O’Brien

J.S.

et al., Penile cancer: Contemporary lymph node management, Journal of Urology 197 (2017), 1387–1395. doi: 10.1016/j.juro.2017.01.059.

Hakenberg

O.W.

et al., The diagnosis and treatment of penile cancer, Deutsches Aerzteblatt Online 115 (2018), 646–652. doi: 10.3238/arztebl.2018.0646.

Dickstein

R.J.

et al., Prognostic factors influencing survival from regionally advanced squamous cell carcinoma of the penis after preoperative chemotherapy, BJU International 117 (2016), 118–125. doi: 10.1111/bju.12946.

10.

et al., Neoadjuvant docetaxel, cisplatin and ifosfamide (ITP) combination chemotherapy for treating penile squamous cell carcinoma patients with terminal lymph node metastasis, BMC Cancer 19 (2019), 1–8. doi: 10.1186/s12885-019-5847-2.

11.

Adhami

et al., Candidate miRNAs in human breast cancer biomarkers: A systematic review, Breast Cancer 25 (2018), 198–205. doi: 10.1007/s12282-017-0814-8.

12.

et al., MiRNA-based therapeutics for lung cancer, Current Pharmaceutical Design 23 (2018), 5989–5996. doi: 10.2174/1381612823666170714151715.

13.

Khan

et al., Role of miRNA-regulated cancer stem cells in the pathogenesis of human malignancies, Cells 8 (2019), 840. doi: 10.3390/cells8080840.

14.

Calin

G.A.

et al., Human microRNA genes are frequently located at fragile sites and genomic regions involved in cancers, Proceedings of the National Academy of Sciences 101 (2004), 2999–3004. doi: 10.1073/pnas.0307323101.

15.

Zhang

et al., microRNAs exhibit high frequency genomic alterations in human cancer, Proceedings of the National Academy of Sciences 103 (2006), 9136–9141. doi: 10.1073/pnas.0508889103.

16.

et al., An integrated approach to reveal miRNAs’ impacts on the functional consequence of copy number alterations in cancer, Scientific Reports 5 (2015). doi: 10.1038/srep11567.

17.

Peng

and Croce

C.M.

, The role of microRNAs in human cancer, Signal Transduction and Targeted Therapy 1 (2016). doi: 10.1038/sigtrans.2015.4.

18.

Soh

et al., Identification and characterization of microRNAs associated with somatic copy number alterations in cancer, Cancers 10 (2018), 1–18. doi: 10.3390/cancers10120475.

19.

Kozomara

Birgaoanu

and Griffiths-Jones

, MiRBase: From microRNA sequences to function, Nucleic Acids Research 47 (2019), 155–162. doi: 10.1093/nar/gky1141.

20.

Kumar Gupta

and Kumar

, HPVbase – A knowledgebase of viral integrations, methylation patterns and microRNAs aberrant expression: As potential biomarkers for Human papillomaviruses mediated carcinomas, Scientific Reports 5 (2015), 1–3. doi: 10.1038/srep12522.

21.

Karagkouni

et al., DIANA-TarBase v8: A decade-long collection of experimentally supported miRNA-gene interactions, Nucleic Acids Research 46 (2018), 239–245. doi: 10.1093/nar/gkx1141.

22.

Vlachos

I.S.

et al., DIANA-miRPath v3.0: Deciphering microRNA function with experimental support, Nucleic Acids Research 43 (2015), 460–466. doi: 10.1093/nar/gkv403.

23.

Rukov

J.L.

et al., Pharmaco-miR: Linking microRNAs and drug effects, Briefings in Bioinformatics 15 (2013), 648–659. doi: 10.1093/bib/bbs082.

24.

Kutmon

et al., PathVisio 3: An extendable pathway analysis toolbox, PLoS Computational Biology 11 (2015), 1–13. doi: 10.1371/journal.pcbi.1004085.

25.

Busso-Lopes

A.F.

et al., Genomic profiling of human penile carcinoma predicts worse prognosis and survival, Cancer Prevention Research 8 (2014), 149–156. doi: 10.1158/1940-6207.CAPR-14-0284.

26.

La-Touche

et al., DNA copy number aberrations, and human papillomavirus status in penile Carcinoma. Clinico-pathological correlations and potential driver genes, PLoS ONE 11 (2016), 1–18. doi: 10.1371/journal.pone.0146740.

27.

Kuasne

et al., Integrative miRNA and mRNA analysis in penile carcinomas reveals markers and pathways with potential clinical impact, Oncotarget 8 (2017), 15294–15306.

28.

Wang

et al., Mutational landscape of penile squamous cell carcinoma in a Chinese population, International Journal of Cancer 145 (2019), 1280–1290. doi: 10.1002/ijc.32373.

29.

Feber

et al., CSN1 somatic mutations in penile squamous cell carcinoma, Cancer Research 76 (2016), 4720–4747. doi: 10.1158/0008-5472.CAN-15-3134.

30.

Tuna

and Amos

C.I.

, Next generation sequencing and its applications in HPV-associated cancers, Oncotarget 8 (2017), 8877–8889. doi: 10.18632/oncotarget.12830.

31.

Kovar

Bierbaumer

and Radic-Sarikas

, The YAP/TAZ pathway in osteogenesis and bone sarcoma pathogenesis, Cells 9 (2020), 972. doi: 10.3390/cells9040972.

32.

et al., MicroRNA expression profiles identify biomarkers for predicting the response to chemoradiotherapy in rectal cancer, Molecular Medicine Reports 18 (2018), 1909–1916. doi: 10.3892/mmr.2018.9215.

33.

Yang

et al., Screening of exosomal miRNAs derived from subcutaneous and visceral adipose tissues: Determination of targets for the treatment of obesity and associated metabolic disorders, Molecular Medicine Reports 18 (2018), 3314–3324. doi: 10.3892/mmr.2018.9312.

34.

Zhang

et al., The chromosome 11q13.3 amplification associated lymph node metastasis is driven by miR-548k through modulating tumor microenvironment, Molecular Cancer 17 (2018), 1–18. doi: 10.1186/s12943-018-0871-4.

35.

Tripathi

et al., Unravelling the role of long non-coding RNA – LINC01087 in breast cancer, Non-Coding RNA Research 5 (2020), 1–10. doi: 10.1016/j.ncrna.2019.12.002.

36.

Cheng

J.C.

Auersperg

and Leung

P.C.K.

, TGF-beta induces serous borderline ovarian tumor cell invasion by activating EMT but triggers apoptosis in low-grade serous ovarian carcinoma cells, PLoS ONE 7 (2012), 1–9. doi: 10.1371/journal.pone.0042436.

37.

Wasserman

et al., SMAD4 loss in colorectal cancer patients correlates with recurrence, loss of immune infiltrate, and chemoresistance, Clinical Cancer Research 25 (2018), 1948–1956. doi: 10.1158/1078-0432.CCR-18-1726.

38.

Pickup

M.W.

et al., BMPR2 loss in fibroblasts promotes mammary carcinoma metastasis via increased inflammation, Molecular Oncology 9 (2015), 179–191. doi: 10.1016/j.molonc.2014.08.004.

39.

Huang

Y.C.

et al., Epigenetic regulation of NOTCH1 and NOTCH3 by KMT2A inhibits glioma proliferation, Oncotarget 8 (2017), 63110–63120. doi: 10.18632/oncotarget.18668.

40.

Huang

et al., SETD7 is a prognosis predicting factor of breast cancer and regulates redox homeostasis, Oncotarget 8 (2017), 94080–94090. doi: 10.18632/oncotarget.21583.

41.

Vougiouklakis

et al., SUV420H1 enhances the phosphorylation and transcription of ERK1 in cancer cells, Oncotarget 6 (2015), 43162–43171. doi: 10.18632/oncotarget.6351.

42.

et al., Downregulation of NSD3 (WHSC1L1) inhibits cell proliferation and migration via ERK1/2 deactivation and decreasing CAPG expression in colorectal cancer cells, OncoTargets and Therapy 12 (2019), 3933–3943. doi: 10.2147/OTT.S191732.

43.

Xiong

et al., Downregulation of KMT2D suppresses proliferation and induces apoptosis of gastric cancer, Biochemical and Biophysical Research Communications 504 (2018), 129–136. doi: 10.1016/j.bbrc.2018.08.143.

44.

W.T.

et al., The O-glycosylating enzyme GALNT2 suppresses the malignancy of gastric adenocarcinoma by reducing EGFR activities, American Journal of Cancer Research 8 (2018), 1739–1751.

45.

Liu

S.Y.

et al., Mucin glycosylating enzyme GALNT2 suppresses malignancy in gastric adenocarcinoma by reducing MET phosphorylation, Oncotarget 7 (2016), 11251–11262. doi: 10.18632/oncotarget.7081.

46.

Gaziel-Sovran

et al., miR-30b/30d Regulation of GalNAc transferases enhances invasion and immunosuppression during metastasis, Cancer Cell 20 (2011), 104–118. doi: 10.1016/j.ccr.2011.05.027.

47.

Chugh

et al., Loss of N-acetylgalactosaminyltransferase 3 in poorly differentiated pancreatic cancer: Augmented aggressiveness and aberrant ErbB family glycosylation, British Journal of Cancer 114 (2016), 1376–1386. doi: 10.1038/bjc.2016.116.

48.

Lee

J.H.

Yun

C.W.

and Lee

S.H.

, Cellular prion protein enhances drug resistance of colorectal cancer cells via regulation of a survival signal pathway, Biomolecules and Therapeutics 26 (2018), 313–321. doi: 10.4062/biomolther.2017.033.

49.

G.R.

, Prion protein is required for tumor necrosis factor α (TNFα)-triggered nuclear factor κb (NF-κB) signaling and cytokine production, Journal of Biological Chemistry 292 (2017), 18747–18759. doi: 10.1074/jbc.M117.787283.

50.

Atkinson

C.J.

et al., EGFR and Prion protein promote signaling via FOXO3a-KLF5 resulting in clinical resistance to platinum agents in colorectal cancer, Molecular Oncology 13 (2019), 725–737. doi: 10.1002/1878-0261.12411.

51.

Nikitovic

et al., Proteoglycans-Biomarkers and targets in cancer therapy, Frontiers in Endocrinology 9 (2018), 1–8. doi: 10.3389/fendo.2018.00069.

52.

Mller

K.H.

et al., Inhibition by cellular vacuolar atpase impairs human papillomavirus uncoating and infection, Antimicrobial Agents and Chemotherapy 58 (2014), 2905–2911. doi: 10.1128/AAC.02284-13.

53.

Hemmat

and Baghi

H.B.

, Human papillomavirus E5 protein, the undercover culprit of tumorigenesis, Infectious Agents and Cancer 13 (2018), 4–5. doi: 10.1186/s13027-018-0208-3.

54.

DiMaio

and Petti

L.M.

, The E5 proteins, Virology 445 (2013), 99–114. doi: 10.1016/j.virol.2013.05.006.

55.

Disbrow

G.L.

Hanover

J.A.

and Schlegel

, Endoplasmic reticulum-localized human papillomavirus type 16 E5 protein alters endosomal pH but not trans-Golgi pH, Journal of Virology 79 (2005), 5839–5846. doi: 10.1128/jvi.79.9.5839-5846.2005.

56.

Miles

A.L.

et al., The vacuolar-ATPase complex and assembly factors, TMEM199 and CCDC115, control HIF1α prolyl hydroxylation by regulating cellular Iron levels, ELife 6 (2017), 1–28. doi: 10.7554/eLife.22693.

57.

Fischer

et al., Human papilloma virus E7 oncoprotein abrogates the p53-p21-DREAM pathway, Scientific Reports 7 (2017), 2603. doi: 10.1038/s41598-017-02831-9.

58.

Kuasne

et al., Penile cancer-derived cells molecularly characterized as models to guide targeted therapies, Cells 10 (2021), 814. doi: 10.3390/cells10040814.

59.

Caetano-Pinto

et al., Cetuximab prevents methotrexate-induced cytotoxicity in vitro through epidermal growth factor dependent regulation of renal drug transporters, Molecular Pharmaceutics 14 (2017), 2147–2157. doi: 10.1021/acs.molpharmaceut.7b00308.

60.

Meads

, Cetuximab for the first-line treatment of metastatic colorectal cancer, Health Technology Assessment 1 (2010), 1–8. doi: 10.3310/hta14suppl1/01.

61.

Bokemeyer

et al., Addition of cetuximab to chemotherapy as first-line treatment for KRAS wild-type metastatic colorectal cancer: Pooled analysis of the CRYSTAL and OPUS randomised clinical trials, European Journal of Cancer 48 (2012), 1466–1475. doi: 10.1016/j.ejca.2012.02.057.

62.

Koyama

et al., Cetuximab and paclitaxel combination therapy for recurrent basaloid squamous cell carcinoma in the ethmoid sinus, Auris Nasus Larynx (2020). doi: 10.1016/j.anl.2020.07.002.

63.

Qin

et al., Efficacy and tolerability of first-line cetuximab plus leucovorin, fluorouracil, and oxaliplatin (FOLFOX-4) versus FOLFOX-4 in patientswith RASwild-typemetastatic colorectal cancer: The open-label, randomized, phase III TAILOR trial, Journal of Clinical Oncology 36 (2018), 3031–3039. doi: 10.1200/JCO.2018.78.3183.

64.

Mishra

et al., N-acetyl-cysteine in combination with celecoxib inhibits Deoxynivalenol induced skin tumor initiation via induction of autophagic pathways in swiss mice. Free Radical Biology and Medicine 156 (2020), 70–82. doi: 10.1016/j.freeradbiomed.2020.06.001.

65.

Velmurugan

B.K.

et al., Combination of celecoxib and calyculin-A inhibits epithelial-mesenchymal transition in human oral cancer cells, Biotechnic and Histochemistry 95 (2020), 341–348. doi: 10.1080/10520295.2019.1700429.

66.

Karai

et al., Celecoxib prevents doxorubicin-induced multidrug resistance in canine and mouse lymphoma cell lines, Cancers 12 (2020), 1117. doi: 10.3390/cancers12051117.