Skip to main content
Download PDF

Coordinated action of human papillomavirus type 16 E6 and E7 oncoproteins on competitive endogenous RNA (ceRNA) network members in primary human keratinocytes

BMC Cancer volume 21, Article number: 673 (2021) Cite this article

Abstract

Background

miRNAs and lncRNAs can regulate cellular biological processes both under physiological and pathological conditions including tumour initiation and progression. Interactions between differentially expressed diverse RNA species, as a part of a complex intracellular regulatory network (ceRNA network), may contribute also to the pathogenesis of HPV-associated cancer. The purpose of this study was to investigate the global expression changes of miRNAs, lncRNAs and mRNAs driven by the E6 and E7 oncoproteins of HPV16, and construct a corresponding ceRNA regulatory network of coding and non-coding genes to suggest a regulatory network associated with high-risk HPV16 infections. Furthermore, additional GO and KEGG analyses were performed to understand the consequences of mRNA expression alterations on biological processes.

Methods

Small and large RNA deep sequencing were performed to detect expression changes of miRNAs, lncRNAs and mRNAs in primary human keratinocytes expressing HPV16 E6, E7 or both oncoproteins. The relationships between lncRNAs, miRNAs and mRNAs were predicted by using StarBase v2.0, DianaTools-LncBase v.2 and miRTarBase. The lncRNA-miRNA-mRNA regulatory network was visualized with Cytoscape v3.4.0. GO and KEEG pathway enrichment analysis was performed using DAVID v6.8.

Results

We revealed that 85 miRNAs in 21 genomic clusters and 41 lncRNAs were abnormally expressed in HPV E6/E7 expressing cells compared with controls. We constructed a ceRNA network with members of 15 lncRNAs – 43 miRNAs – 358 mRNAs with significantly altered expressions. GO and KEGG functional enrichment analyses identified numerous cancer related genes, furthermore we recognized common miRNAs as key regulatory elements in biological pathways associated with tumorigenesis driven by HPV16.

Conclusions

The multiple molecular changes driven by E6 and E7 oncoproteins resulting in the malignant transformation of HPV16 host cells occur, at least in part, due to the abnormal alteration in expression and function of non-coding RNA molecules through their intracellular competing network.

Peer Review reports

Background

Human papillomaviruses (HPVs) are small, double-stranded circular DNA viruses with significant oncogenic potential. HPVs infect epithelial cells of skin and mucous membranes mainly in the anogenital and head and neck region. High-risk HPV-types, such as HPV16, 18, 31, 33, 45, 51, 58 etc., are main etiologic factors of numerous different carcinomas including cervical, vulvar, penile, anal and head and neck cancers. Nearly 100% of cervical carcinomas are positive for high-risk HPVs, HPV16 being the most prevalent type detected in these malignancies [1, 2]. E6 and E7 early viral proteins are the major oncoproteins of high-risk HPVs. The transforming activity of E6 and E7 is mainly mediated through binding and inactivation of cellular tumour suppressor proteins p53 and pRb, respectively, although they are able to directly or indiretly interact with numerous cellular factors involved in regulation of cell cycle, cell death or cellular differentiation, thereby promote malignant transformation [3,4,5].

High-throughput RNA-seq studies have revealed the great complexity of human transcriptome by demonstrating that only approximately 2% of the human genome encodes protein-coding sequences. A significant portion of the human genes is actively transcribed but not translated into functional proteins suggesting their crucial regulatory role in biological processes under different conditions. Among these transcripts, several classes of non-coding RNA molecules have been identified including microRNAs (miRNAs), small nucleolar RNAs (snoRNAs), PIWI-interacting RNAs (piRNAs) and long non-coding RNAs (lncRNAs) [6, 7]. The approximately 18–25 nt miRNAs are the most intensively studied and best-known subtype of non-coding RNAs. miRNAs downregulate the expression of their targets post-transcriptionally by binding to the 3′-UTR of mRNAs leading to the translational inhibition or degradation of their target. Expression of miRNAs may vary in different tissues, developmental stages or pathological conditions [8,9,10,11,12]. lncRNAs are a diverse group of non-coding transcripts longer than 200 nt which are thought to have important function in gene expression regulation. Their ability to interact with DNA, RNA or proteins suggests their multiple functions in diverse biological processes. They have been linked to epigenetic mechanisms, splicing, translational regulation and protein stability, although the accurate mechanism of their action remains largely unknown [9, 13,14,15,16].

Emerging evidence shows the aberrant function of non-coding RNAs especially miRNAs and lncRNAs in various diseases and tumours including HPV-induced cervical cancers. The non-coding transcripts involved in cell cycle, differentiation, proliferation, migration and apoptosis tend to have altered expression during the initiation and progression of tumorigenesis [9, 15, 17,18,19,20,21,22,23,24,25,26,27]. Although several details of the complex pathogenesis of HPV-induced carcinogenesis have been explored primarily by focusing on oncogenic activity of E6 and E7 proteins of high-risk HPV types, the complete mechanism is not fully understood [3, 28]. An increasing number of studies have reported both abnormal upregulation and downregulation of numerous miRNAs in cervical cancer cell lines and tissue samples, strongly suggesting the oncogenic or tumour suppressor function of certain miRNAs in HPV-associated carcinogenesis [29,30,31,32,33,34,35,36,37,38]. As for the lncRNAs, despite the growing body of evidence regarding their importance in oncogenesis, only a limited number of studies have directly linked certain lncRNAs to HPV-associated cervical carcinoma [15, 21, 39,40,41,42,43]. Thus, the modulation of a wide range of host lncRNAs by E6 and E7 oncoproteins of high-risk HPVs remains to be extensively investigated.

According to the newly established competing endogenous RNA (ceRNA) hypothesis, all types of RNAs, both protein coding and non-coding transcripts are able to “talk to each other” and influence the functions of the others as a part of a complex intracellular regulatory network. This communication is thought to be mediated by miRNAs through binding to microRNA response elements (MREs) on various RNA transcripts [44]. Recent studies have proved that lncRNAs have the ability to regulate miRNA activity by acting as decoys for them, leading to the alteration of protein turnover. In this context, the altered interaction between lncRNAs and miRNAs by perturbation in their levels may contribute to the pathogenesis of cancer [9, 45,46,47,48,49].

In the current study, we performed high-throughput new generation RNA sequencing to detect cellular lncRNA and miRNA expression changes in primary human keratinocytes expressing HPV16 E6, E7 and E6/E7 oncoproteins. In addition to non-coding RNAs, mRNA expression changes were measured and we constructed a ceRNA network of coding and non-coding genes to suggest a regulatory network associated with HPV16 infection. We also refined the mRNAs of this ceRNA network by key biological pathways involved in cancer development as determined by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.

Methods

Cell culture and RNA isolation

Primary human foreskin keratinocytes (HFK) transduced with LXSN-based retroviral vectors encoding HPV16 oncogenes were obtained from a previous study of our research group [50]. HPV16 E6, HPV16 E7 and HPV16 E6/E7 expressing cells and LXSN control cells were cultured in Defined Keratinocyte-Serum Free Medium (DK-SFM, Invitrogen) as described previously [51]. Two independent, passage-matched HFK populations were used in our experiments. Total RNA was isolated from each transduced keratinocyte cultures using TriReagent (Sigma-Aldrich) according to the manufacturer’s instructions. HPV16 E6 and E7 specific reverse transcription PCR was performed to confirm the expression of viral oncogenes in all HFK samples [52].

RNA sequencing and data analysis

RNA integrity was checked on an Agilent BioAnalyzer 2100 instrument using RNA 6000 Nano kit (Agilent Technologies) according to manufacturer’s protocol and samples with RNA integrity number > 7.0 were accepted to good quality and were used in the further experiments. A NanoDrop ND-1000 was used to determine RNA concentration. Global transcriptome analysis was performed on Illumina sequencing platform. Sequencing libraries for small RNA-Seq were generated from 1 μg total RNA using TruSeq Small RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. Fragment size distribution and molarity of libraries were checked on Agilent BioAnalyzer DNA1000 chip (Agilent Technologies, Santa Clara, CA, USA). Then single read 50 bp sequencing run was performed on Illumina HiScan SQ instrument (Illumina, San Diego, CA, USA). cDNA library for RNA-Seq was generated from 1 μg total RNA using TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. Briefly, poly-A tailed RNAs were purified by oligodT conjugated magnetic beads and fragmented on 94 C degree for 8 min, then 1st strand cDNA was transcribed using random primers and SuperScript II reverse transcriptase (Lifetechnologies, Carslbad, CA, USA). Following this step second strand cDNA synthesized, double stranded cDNA end repaired and 3′ ends adenylated then Illumina index adapters were ligated. After adapter ligation enrichment PCR was performed to amplify adapter ligated cDNA fragments. Fragment size distribution and molarity of libraries were checked on Agilent BioAnalyzer DNA1000 chip (Agilent Technologies, Santa Clara, CA, USA). Single read 50 bp sequencing run was performed on Illumina HiScan SQ instrument (Illumina, San Diego, CA, USA). Quality control and alignment of raw reads to the reference human genome hg19 was done using FASTQC package. Quantification and Fold Change analysis (FC) were performed using StrandNGS software. Expressed genes were identified by filtering out the lowest 20 percentile of genes based on raw signal intensity. RNA sequencing and data analysis were carried out by Genomic Medicine and Bioinformatic Core Facility, Deparment of Biochemistry and Molecular Biology, Faculty of Medicine, University of Debrecen (Debrecen, Hungary).

Quantitative real-time RT-PCR

For miRNA RT-qPCR, cDNA was synthesised using TaqMan MicroRNA Reverse Transcription Kit (Applied Biosystems) with miRNA-specific stem-loop primers. Mature miRNAs were detected on 7500 Real Time PCR System (Applied Biosystems) using TaqMan Universal Master Mix and MicroRNA Assays (Applied Biosystems) according to manufacturer’s protocol. PCR reactions were performed in triplicate, small RNA RNU44 was used as an endogenous control.

For large RNA RT-qPCR, the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems) was used to prepare cDNA. qPCR was performed using TaqMan Gene Expression Master Mix and Gene Expression Assays (Applied Biosystems) according to the manufacturer’s instructions. Each reaction was performed in triplicate and GAPDH (glyceraldehyde 3-phosphate dehydrogenase) was used as endogenous control. In data analysis, comparative Ct method was used to obtain the Relative Quantification (RQ) values with standard deviation and confidence intervals (7500 System SDS Software, version 2.0.6).

Construction of ceRNA network, GO and KEGG pathway analysis

The construction of a regulatory network was a multistep process. At first, we identified miRNAs, lncRNAs and mRNAs which were differentially expressed in HFK-E6/E7 cells. ≥10 reads in raw data and at least ±1.5 fold change were applied as treshold cutoffs in data analysis. For each transcript, z-score and p value were calculated with Benjamini-Hochberg FDR correction of 0.10 to include the transcript into downstream analyses. lncRNA-miRNA interactions were identified using StarBase v2.0 and DianaTools-LncBase v.2 databases which provides experimentally supported interaction data [53, 54]. In the following step, mRNAs targeted by miRNAs were explored using the advanced search function of miRTarBase [55]. To establish the lncRNA-miRNA-mRNA network, we integrated lncRNA-miRNA and miRNA-mRNA interactions and the regulatory network was visualized with Cytoscape v3.4.0 [56]. In RNA-RNA interaction search, we used only the experimentally validated search module with default settings so only miRNA-lncRNA and miRNA-mRNA pairs were included in further analyses, which were previously experimentally validated and were differentially expressed in our experimental system, too. The lists of StarBase v2.0, DianaTools-LncBase v.2 and miRTarBase search results are shown in Additional file 1.

The coding genes involved in ceRNA network were input into the DAVID v6.8 (Database for Annotation, Visualization and Integrated Discovery) for GO (Gene Ontology) and KEEG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis [57].

Results

Altered miRNA profile in HPV16 E6/E7 expressing keratinocytes

Next generation microRNA sequencing was used to investigate the alteration of miRNA expression levels in response to high risk HPV16 E6/E7 in undifferentiated primary human keratinocytes. Two independent passage-matched populations of human foreskin kerationcytes expressing HPV16 E6 and/or E7 or control vector (HFK-E6, HFK-E7, HFK-E6/E7 and HFK-LXSN) were used in our experiments. Cells were cultured in low-calcium, serum-free medium and harvested for RNA extraction at maximum 80% confluency after up to 5 passages to keep them undifferentiated and to avoid acquiring a transformed phenotype. In addition, E6 and E7 specific RT-PCR in this study (Additional file 2) and p53 and Rb Western blot analyses in our previous study verified the presence of functionally active oncoproteins in the transduced cells [58]. NGS analysis of cellular mRNAs - that was performed in this study - provides further evidence for the presence of functionally active E6 and E7 oncoproteins in the transduced cells. The effects of the HPV16 E6 and E7 on the level of well-known transcriptional targets of high-risk HPV oncoproteins (TERT, CDKN1A, CDKN2A, MDM2, etc.) are shown in Additional file 3. As a consequence, expression changes in RNA levels most probably represent the effect of HPV oncogenes rather than the effect of differentiation or transformation. Only miRNAs which matched the treshold criteria and showed consistent expression changes in the two independent experiments were involved in downstream analyses.

85 differentially expressed miRNAs were detected in E6/E7 expressing HFKs compared to control vector cells (E6/E7 vs. LXSN). Of these dysregulated miRNAs, 44 miRNAs were downregulated (e.g. miR-34a-5p, miR-127-3p, miR-138-5p, miR-370-3p, miR-377-5p, miR-487b-3p, miR-654-3p) and 41 miRNAs (e.g. let-7c-5p, miR-10-5p, miR − 20b-5p, miR − 132-3p, miR-195-5p, miR-675-5p, miR-3605-5p) were upregulated (Fig. 1).

Fig. 1
figure 1

Differentially expressed miRNAs in HPV16 E6/E7 expressing keratinocytes compared to control vector cells. Average Log2 transformed FC values of up - and downregulated miRNAs are shown

Full size image

To validate the miRNA-seq results, altered expression of several miRNAs was confirmed using TaqMan real-time assays in three independent experiments. The RNA preparations used as template in 2 of 3 TaqMan assays were also used in the RNA sequencing experiments to verify the validity of RNA sequencing results. RNA preparation from an additional HFK population was used in the third TaqMan assay to further strengthen our observations in RNA sequencing. The downregulation of miR-34a-5p, miR-138-5p and miR-370-3p and the upregulation of miR-132-3p and miR-675-5p matched with the expression tendency seen in the sequencing results and their expression changes were statistically significant. The integrated results of sequencing and TaqMan analyses of selected miRNAs are shown in Fig. 2.

Fig. 2
figure 2

Validation of selected miRNAs identified by miRNA sequencing. The integrated results of sequencing and TaqMan analyses of miRNAs are shown. Green lines represent FC values obtained in miRNA sequencing. Altered expression of selected miRNAs was confirmed using TaqMan real-time assays in HPV16 E6, HPV16 E7, HPV16 E6/E7 and LXSN control cells in three independent experiments. In data analysis, comparative Ct method was used to obtain the Relative Quantification (RQ) values with standard deviation and confidence intervals (95%). The RQ value of LXSN control was set to 1. Representative graphs are shown. *p < 0.01

Full size image

We investigated the changes in miRNA expression levels in HFK cells expressing only the E6 or E7 alone to understand the importance of individual oncoproteins in miRNA profile alteration. We considered the miRNA expression changes as a result of E6 oncoprotein if the expression changes were predominant in E6 vs LXSN and E6/E7 vs E7 comparisons in addition to E6/E7 vs LXSN results in both experiments. Similarly, E7 effect was supposed if the altered miRNA expression was observed mainly in E7 vs LXSN and E6/E7 vs E6 comparisons. In 56 of 85 dysregulated miRNAs, the modulation of expression was considered primarily the individual effect of HPV16 E6 oncoprotein. In contrast, miRNA expression alterations were observed to be driven predominantly by E7 only in 5 cases while up- or downregulation of 24 miRNAs was driven both by E6 and E7 together. Hence, majority of miRNA expression alterations in E6/E7 expressing keratinocytes occured as a result of the independent expression of HPV16 E6 or the cumulative effect of both E6 and E7 (Additional file 4).

Effect of HPV16 E6/E7 and human miRNA clusters

We analyzed the effect HPV16 E6 and E7 oncoproteins by human miRNA clusters. The 85 miRNAs with altered expression were located to known miRNA clusters on the basis of miRBase, and 55 were identified as belonging to genomic clusters [59]. Out of the 55 clustered miRNAs, 20 were upregulated and 35 were downregulated. In total, the expression of 21 miRNA clusters were modulated by HPV16 E6/E7 oncoproteins. 21 differentially expressed miRNAs were members of small clusters containing only 2–4 miRNAs, while 34 miRNAs were members of larger clusters containing 6 to 22 miRNAs. According to the genomic location of miRNA clusters, chromosome X and chromosome 14 seem to have an importance consisting 5–5 clusters with altogether 38 altered miRNA members. The directions of miRNA expressional changes were either consistently up or consistently down within all but one clusters containing multiple miRNAs with altered expression (Table 1).

Table 1 miRNA clusters modulated in HPV16 E6/E7 expressing HFKs
Full size table

Long non-coding RNAs modulated by HPV oncoproteins

In addition to miRNA analysis, polyA-tailed long RNA (≥ 200 nt) sequencing was performed from the same HFK populations in two independent experiments to investigate the modulatory effect of HPV16 oncoproteins on long non-coding RNA expressions. The same HPV oncoprotein expressing keratinocytes, as well as the same validation criteria were applied in the analyses of our data. Only lncRNAs which matched the treshold criteria and showed consistent expression results in the two independent experiments were involved in downstream analyses.

41 differentially expressed lncRNAs were detected in E6/E7 expressing HFKs compared to control vector cells (E6/E7 vs. LXSN). Of these, 14 lncRNAs were downregulated (e.g. MEG3, MEG9, LINC00520, LINC00923, FBXL19-AS1) and 27 lncRNAs (e.g. H19, CECR7, LINC00087, LINC00925, DLX6-AS1, HOXA-AS2) were upregulated (Fig. 3).

Fig. 3
figure 3

Differentially expressed lncRNAs in HPV16 E6/E7 expressing keratinocytes compared to control vector cells. Average Log2 transformed FC values of up - and downregulated lncRNAs are shown

Full size image

TaqMan assays of selected lncRNAs also confirmed the results of sequencing. Significant downregulation of MEG3 and upregulation of CECR7 and H19 are represented on Fig. 4.

Fig. 4
figure 4

Validation of selected lncRNAs identified by large RNA seqencing. The integrated results of sequencing and TaqMan analyses of lncRNAs are shown. Green lines represent FC values obtained in large RNA sequencing. Altered expression of selected lncRNAs was confirmed using TaqMan real-time assays in HPV16 E6, HPV16 E7, HPV16 E6/E7 and LXSN control cells in three independent experiments. In data analysis, comparative Ct method was used to obtain the Relative Quantification (RQ) values with standard deviation and confidence intervals (95%). The RQ value of LXSN control was set to 1. Representative graphs are shown. *p < 0.01

Full size image

We applied the same criteria that were used in the case of miRNAs when we investigated the influence of individual oncoproteins on lncRNA expression changes. Like in the case of miRNAs, majority of the modifications in lncRNA expressions was either the dominant effect of E6 oncoprotein alone or the outcome of the presence of E6 and E7 together. The importance of E7 individually was observed only in two cases, while expressions of 21 lncRNAs were altered predominantly by E6 and expression changes of 18 lncRNAs were driven by the two oncoproteins together (Additional file 4).

Coordinated action of HPV16 E6/E7 oncoproteins on ceRNA members: construction of regulatory network, GO and KEGG analyses

Predicted pairwise lncRNA-miRNA interactions were identified between the 85 differentially expressed miRNAs and the 41 differentially expressed lncRNAs using StarBase v2.0 and DianaTools-LncBase v.2 experiment - supported databases. As a result, 43 of 85 specific miRNAs and 15 of 41 specific lncRNAs were identified to interact in this experimental system. A total of 77 interactions between 15 lncRNAs and 43 miRNAs were identified based on combined results of database search. MEG3 lncRNA was found the most potent interaction partner with 29 miRNAs, followed by H19 lncRNA with 17 different miRNA interactions. Among miRNAs, miR-20b-5p, miR-143-3p, miR − 497-5p and miR-34a-5p, miR-195-5p were the most common molecules with interactions with 5 and 4 different lncRNAs, respectively (Table 2).

Table 2 lncRNA – miRNA interactions with potential functional significance identified in HPV16 E6/E7 expressing HFKs
Full size table

To construct the whole lncRNA-miRNA-mRNA regulatory network, we searched mRNAs targeted by miRNAs using the miRTarBase. We mapped miRNA - mRNA pairs with experimental evidence using the advanced search function of miRTarBase. We investigated the 43 miRNAs that are in interactions with lncRNAs and mRNAs which showed differential expressions by HPV16 E6/E7. A total of 358 mRNAs targeted by 43 miRNAs were identified with expressions changed by E6/E7 oncoproteins (at least ±1.5 fold change in E6/E7 HFKs similarly to miRNAs and lncRNAs) and a total of 522 miRNA-mRNA pairs were observed (Table 3). After combination of lncRNA-miRNA and miRNA-mRNA interactions, we constructed a ceRNA network with members of 15 lncRNAs – 43 miRNAs – 358 mRNAs with significantly altered expressions driven by HPV16 oncoproteins. The regulatory network was visualized with Cytoscape v3.4.0. (Fig. 5).

Table 3 miRNA-mRNA target pairs with potential functional significance identified in HPV16 E6/E7 expressing HFKs
Full size table
Fig. 5
figure 5

Construction of ceRNA network of differentially expressed lncRNAs-miRNAs-mRNAs in HPV16 E6/E7 expressing HFKs. Interactions of 15 lncRNAs – 43 miRNAs – 358 mRNAs are represented. The regulatory network was visualized with Cytoscape v3.4.0

Full size image

In order to understand the biological processes affected by ceRNA regulatory network, the mRNAs were submitted to GO and KEGG analyses using DAVID v6.8 database. Functional enrichment analysis revealed 95 of 358 genes as cancer – associated genes, such as E2F2, MDM2, CDC25C, CDKN1A, CDKN1C, CDKN2A, GAS1, IGF2, NOTCH1, TERT etc.. According to KEGG analysis the genes in the regulatory network are involved in pathways such as pathways in cancer, signaling pathways regulating pluripotency of stem cells, proteoglycans in cancer, microRNAs in cancer, TNF signaling pathway, Hippo signaling pathway, viral carcinogenesis, Wnt signaling pathway. In the GO analysis (GO, level BP3), key biological processes which may play important roles in HPV-induced carcinogenesis were significantly enriched such as signal transduction, regulation of cell communication, regulation of cell cycle, cell cycle arrest, cell death, regulation of cell death, cell migration, regulation of cell proliferation, epithelial cell proliferation, etc.. The lists of genes significantly enriched GO and KEGG pathways are reported as additional file (Additional file 5).

We examined the effect of HPV16 E6/7 on the expression of ceRNA members involved in specific biological processes important in cancer development to see if there are „hot points” in the biological activity of oncoproteins. Nine common miRNAs targeting mRNAs which are components of key canonical pathways were identified. miR − 20b-5p, miR-34a-5p, miR-106a-5p, miR-125b-5p, miR-143-3p, miR-181c-5p, miR-335-5p, miR-363-3p and miR-432-5p were observed as key regulatory elements in regulation of cell cycle and cell death, regulation of cell proliferation, regulation of epithelial cell proliferation, in such biological processes which have an impact in cancer development (Table 4, Fig. 6). The integrated ceRNA network of above mentioned biological processes consists of 11 lncRNAs and 82 mRNAs targeted by common miRNAs in the regulatory modules (Fig. 7). These results support the existence of regulatory „hot points” in key cancer related pathways driven by the coordinatied action of HPV16 E6 and E7 oncoproteins on ceRNA members.

Table 4 Common miRNAs with functional importance in four cancer related regulatory pathways
Full size table
Fig. 6
figure 6

Venn diagram of dysregulated miRNAs involved in particular biological processes with significance in cancer development

Full size image
Fig. 7
figure 7

The integrated ceRNA network of cancer related biological processes in HPV16 E6/E7 expressing cells. The network consists of 11 lncRNAs and 82 mRNAs targeted by nine common miRNAs in the regulatory modules. Diamonds represent lncRNAs, rectangles represent miRNAs and ovals represent mRNAs. RNAs with red were downregulated and RNAs with blue were upregulated. The regulatory network was visualized with Cytoscape v3.4.0

Full size image

Discussion

E6 and E7 oncoproteins of high-risk HPVs cause large-scale gene expression changes in the HPV infected epithelial cells. They interact with numerous signaling molecules, thereby they are able to reprogram multiple cellular regulatory pathways involved in cell cycle, cell death, cell differentiation and migration playing a crucial role in HPV-induced abnormal cell proliferation and malignant transformation [5, 58, 60]. In addition to protein coding genes, E6 and E7 can also alter the expression of non-coding RNAs, including lncRNAs and miRNAs [29, 39].

Host cell’s non-coding RNAs, especially miRNAs and lncRNAs function as key regulators of cellular biological processes both under physiological and pathological conditions, including tumour initiation and progression by multiple types of cancers including HPV-associated malignancies. However, the majority of the studies have focused on high-risk HPV associated expressional changes and/or interactions of either miRNAs or lncRNAs. Several studies indicated that miRNAs and lncRNAs may interact with and regulate cellular processes involved in cancer development through modulation of the expression of their downstream mRNA target [47, 48, 61]. Nevertheless, non-coding RNA molecules, both miRNAs and lncRNAs, have the potential to interact many other RNAs and modulate numerous mRNA targets and vice versa, a single mRNA may be targeted by hundreds of miRNAs. Therefore, the aberrant expression of these RNA molecules will influence multiple signaling pathways in HPV infected cells. Thus, analyzing the global context of coding and non-coding RNAs provides valuable information on the complexity of the structure and function of endogenous RNA network regarding high-risk HPV infection. Therefore our goal was to investigate the global expression changes of miRNAs, lncRNAs and mRNAs driven by HPV16 E6/E7 oncoproteins. We described the network of abnormally expressed non-coding RNAs in stably transduced undifferentiated human primary keratinocytes using small RNA (≤ 200 nt) and large RNA (≥ 200 nt) deep sequencing. To understand the potential impact of expression alterations of miRNAs and lncRNAs modulated by E6/E7 oncoproteins on cellular mRNA abundance and functions, we have paired the expression data of the various RNA populations. Furthermore, additional GO and KEGG analyses were performed to understand the consequences of mRNA expression alterations on biological processes. To our best knowledge, to date this is the first study which describes such an integrated expression data of non-coding RNAs and mRNAs and construct a ceRNA network of miRNAs, lncRNAs and mRNAs in association with high-risk HPV infection. It is important to highlight, that all of our functional analyses were performed on the basis of E6/E7 expressing HFK cells, because this scenario represents best the in vivo HPV16 infection in epithelial cells. Beyond that, we examined the individual effect of E6 and E7 oncogenes on non-coding RNA expression, too. It is well documented, that in high-risk HPV-types, including HPV16, E6 and E7 expression is controlled by alternative splicing resulting in three different E6 isoforms (E6*I, E6*II and E6^E7) besides the unspliced full-length E6. The expression of different HPV16 isoforms may vary among cells and tissues, nevertheless E6*I is reported to be the most abundant isoform in cervical tumors and cells [62]. Two different E6 products were detected in HPV16 E6 specific RT-PCR performed in this study which matched the full-length E6 and an E6 splice isoform. Only the full length E6 transcript is responsible for coding the transforming E6 oncoprotein while different E6 isoforms may have an antagonistic effect in tumorigenesis. However, both the regulation of these isoforms and their effect on non-coding RNA expression is poorly understood, there are limited number of studies reporting relationship between miRNAs and E6 isoforms [63]. Although, we were not able to differentiate the effect of full-length E6 and E6 splice isoform on the expression of different RNA populations in our study, our observations likely show the cumulative effect of different E6 isoforms and E7 oncoprotein on RNA expression in our experimental system which may be significant in tumorigenesis. At the same time, further research needs to be performed to clarify the exact role of E6 isoforms on non-coding RNA expression in HPV16 host cells.

In our study, we revealed that 85 miRNAs and 41 lncRNAs were abnormally expressed in HPV E6/E7 expressing cells compared with controls. When we examined the impact of oncoproteins on RNA alteration profile individually, we found that the vast majority of RNA expression changes was driven primarly by E6 oncoprotein, which was not a surprise as E6 has a strong effect on host cell processes [64]. In total, the expression of 21 miRNA clusters were modulated by HPV16 E6/E7 oncoproteins, affecting 55 differentially expressed miRNAs. According to chromosome location, 38 miRNA members of miRNA clusters on chromosome X and chromosome 14 may have functional importance in HPV16 infected cells, seeing the clustered miRNAs often coordinately regulate certain biological processes. Harden and coworkers also described modulation of multiple miRNA clusters by HPV16 E6/E7 oncoproteins and they found several miRNAs with expression changes on chromosome X and chromosome 14 [29].

This study also revealed that several differentially expressed non-coding RNAs tend to have regulatory interaction with multiple RNA species also altered in the presence of HPV16 oncogenes. We constructed a ceRNA network with members of 15 lncRNAs – 43 miRNAs – 358 mRNAs with significantly altered expressions driven by HPV16 oncoproteins. Among non-coding RNAs, lncRNAs such as MEG3 and H19, and miRNAs such as miR − 20b-5p, miR-143-3p, miR-497-5p, miR-34a-5p and miR-195-5p had the highest number of relationships in the interaction network. In accordance with our results, downregulation of MEG3 and overexpression of H19 lncRNAs were found by several research groups in association with various neoplasms and the link of these lncRNAs with p53 and pRb pathways underlines their significance in E6 and E7 mediated transformation [19, 65, 66]. The relationships between lncRNAs and miRNAs, and miRNAs and mRNAs were predicted by using StarBase v2.0, DianaTools-LncBase v.2 and miRTarBase, which are experiment - supported databases corroborating that the RNA-RNA interactions identified in our experiment would occur not only in silico situation [53,54,55]. These supporting data justified the construction of a ceRNA network from the experimental data of this study.

The functional enrichment analyses identified numerous ceRNA member coding genes as cancer related genes and the majority of mRNAs are involved in GO biological processes and KEGG pathways associated with cancer initiation and progression. Furthermore, our comprehensive examination recognized nine miRNAs as key regulatory elements in regulation of cell cycle and cell death, regulation of cell proliferation and regulation of epithelial cell proliferation. Although, all of these miRNAs have been reported in previous studies associated with HPV infection, their functions in such a complex ceRNA network were so far unknown and they are worth to be further investigated as potential biomarkers in disease detection and progression.

Conclusions

In this study, our main goal was to highlight the complexity of regulatory network of cellular non-coding RNAs in HPV16 infected cells. In the last few years, meta-analyses of sequencing data from different databases or analyzing the global context of coding and non-coding RNAs using integrative methods have become a new and widespread direction of current cancer research. Nevertheless, carrying further experiments focusing on real biological impact of certain RNA-RNA interactions is very useful and essential to explore the exact functions and regulations in depth of non-coding RNAs in cancer development.

We hypothesize, that the multiple molecular changes driven by E6 and E7 oncoproteins resulting in the malignant transformation of HPV16 infected host cells may occur, at least in part, due to the abnormal alteration in expression and function of non-coding RNA molecules through their intracellular competing network. Although several lncRNAs and miRNAs in this study have been previously reported regarding both to HPV-associated and non-HPV cancers, to our best knowlege, this is the first study which describes a ceRNA network of lncRNAs-miRNAs-mRNAs in HPV16 infected cells. Functional enrichment analyses carried out in this study assume important regulatory roles of certain miRNAs and lncRNAs in key biological processes involved in tumour development. We believe, that our comprehensive results will provide a basis for further detailed experiments. Based on these expressional data, we plan to perform further experiments in the future to investigate the relationships of certain non-coding RNAs which funtion is so far poorly understood. Furthermore, miRNAs and lncRNAs can serve as promising biomarkers for molecular diagnosis and progression of cervical cancers.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Abbreviations

HPV:

Human papillomavirus

miRNA:

microRNA

lncRNA:

long non-coding RNA

ceRNA:

competitive endogenous RNA

MRE:

microRNA response element

HFK:

Human foreskin keratinocyte

GO:

Gene Ontology

KEEG:

Kyoto Encyclopedia of Genes and Genomes

References

  1. Faridi R, Zahra A, Khan K, Idrees M. Oncogenic potential of human papillomavirus (HPV) and its relation with cervical cancer. Virol J. 2011;8(1):269. https://doi.org/10.1186/1743-422X-8-269.

    Article  PubMed  PubMed Central  Google Scholar 

  2. zur Hausen H. Papillomaviruses and cancer: from basic studies to clinical application. Nat Rev Cancer 2002;2(5):342–350, DOI: https://doi.org/10.1038/nrc798.

  3. de Freitas AC, Coimbra EC, Leitao MC. Molecular targets of HPV oncoproteins: potential biomarkers for cervical carcinogenesis. Biochim Biophys Acta. 2014;1845(2):91–103. https://doi.org/10.1016/j.bbcan.2013.12.004.

    Article  CAS  PubMed  Google Scholar 

  4. Korzeniewski N, Spardy N, Duensing A, Duensing S. Genomic instability and cancer: lessons learned from human papillomaviruses. Cancer Lett. 2011;305(2):113–22. https://doi.org/10.1016/j.canlet.2010.10.013.

    Article  CAS  PubMed  Google Scholar 

  5. Moody CA, Laimins LA. Human papillomavirus oncoproteins: pathways to transformation. Nat Rev Cancer. 2010;10(8):550–60. https://doi.org/10.1038/nrc2886.

    Article  CAS  PubMed  Google Scholar 

  6. Mattick JS. The genetic signatures of noncoding RNAs. PLoS Genet. 2009;5(4):e1000459. https://doi.org/10.1371/journal.pgen.1000459.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Ponting CP, Belgard TG. Transcribed dark matter: meaning or myth? Hum Mol Genet. 2010;19(R2):R162–8. https://doi.org/10.1093/hmg/ddq362.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005;120(1):15–20. https://doi.org/10.1016/j.cell.2004.12.035.

    Article  CAS  PubMed  Google Scholar 

  9. Liz J, Esteller M. lncRNAs and microRNAs with a role in cancer development. Biochim Biophys Acta. 2016;1859(1):169–76. https://doi.org/10.1016/j.bbagrm.2015.06.015.

    Article  CAS  PubMed  Google Scholar 

  10. Melo SA, Esteller M. Dysregulation of microRNAs in cancer: playing with fire. FEBS Lett. 2011;585(13):2087–99. https://doi.org/10.1016/j.febslet.2010.08.009.

    Article  CAS  PubMed  Google Scholar 

  11. Ribeiro J, Sousa H. MicroRNAs as biomarkers of cervical cancer development: a literature review on miR-125b and miR-34a. Mol Biol Rep. 2014;41(3):1525–31. https://doi.org/10.1007/s11033-013-2998-0.

    Article  CAS  PubMed  Google Scholar 

  12. Romero-Cordoba SL, Salido-Guadarrama I, Rodriguez-Dorantes M, Hidalgo-Miranda A. miRNA biogenesis: biological impact in the development of cancer. Cancer biology & therapy. 2014;15(11):1444–55. https://doi.org/10.4161/15384047.2014.955442.

    Article  CAS  Google Scholar 

  13. Hosseini ES, Meryet-Figuiere M, Sabzalipoor H, Kashani HH, Nikzad H, Asemi Z. Dysregulated expression of long noncoding RNAs in gynecologic cancers. Mol Cancer. 2017;16(1):107. https://doi.org/10.1186/s12943-017-0671-2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Hu W, Alvarez-Dominguez JR, Lodish HF. Regulation of mammalian cell differentiation by long non-coding RNAs. EMBO Rep. 2012;13(11):971–83. https://doi.org/10.1038/embor.2012.145.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Peng L, Yuan X, Jiang B, Tang Z, Li GC. LncRNAs: key players and novel insights into cervical cancer. Tumour Biol. 2016;37(3):2779–88. https://doi.org/10.1007/s13277-015-4663-9.

    Article  CAS  PubMed  Google Scholar 

  16. Rinn JL, Chang HY. Genome regulation by long noncoding RNAs. Annu Rev Biochem. 2012;81(1):145–66. https://doi.org/10.1146/annurev-biochem-051410-092902.

    Article  CAS  PubMed  Google Scholar 

  17. Fu Y, Biglia N, Wang Z, Shen Y, Risch HA, Lu L, et al. Long non-coding RNAs, ASAP1-IT1, FAM215A, and LINC00472, in epithelial ovarian cancer. Gynecol Oncol. 2016;143(3):642–9. https://doi.org/10.1016/j.ygyno.2016.09.021.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Jones M, Lal A. MicroRNAs, wild-type and mutant p53: more questions than answers. RNA Biol. 2012;9(6):781–91. https://doi.org/10.4161/rna.20146.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Kitagawa M, Kitagawa K, Kotake Y, Niida H, Ohhata T. Cell cycle regulation by long non-coding RNAs. Cell Mol Life Sci. 2013;70(24):4785–94. https://doi.org/10.1007/s00018-013-1423-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Li DS, Ainiwaer JL, Sheyhiding I, Zhang Z, Zhang LW. Identification of key long non-coding RNAs as competing endogenous RNAs for miRNA-mRNA in lung adenocarcinoma. Eur Rev Med Pharmacol Sci. 2016;20(11):2285–95.

    PubMed  Google Scholar 

  21. Qin R, Chen Z, Ding Y, Hao J, Hu J, Guo F. Long non-coding RNA MEG3 inhibits the proliferation of cervical carcinoma cells through the induction of cell cycle arrest and apoptosis. Neoplasma. 2013;60(5):486–92. https://doi.org/10.4149/neo_2013_063.

    Article  CAS  PubMed  Google Scholar 

  22. Wilting SM, Verlaat W, Jaspers A, Makazaji NA, Agami R, Meijer CJ, et al. Methylation-mediated transcriptional repression of microRNAs during cervical carcinogenesis. Epigenetics. 2013;8(2):220–8. https://doi.org/10.4161/epi.23605.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Yan X, Hu Z, Feng Y, Hu X, Yuan J, Zhao SD, et al. Comprehensive genomic characterization of long non-coding RNAs across human cancers. Cancer Cell. 2015;28(4):529–40. https://doi.org/10.1016/j.ccell.2015.09.006.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Zhang J, Fan D, Jian Z, Chen GG, Lai PB. Cancer specific long noncoding RNAs show differential expression patterns and competing endogenous RNA potential in hepatocellular carcinoma. PLoS One. 2015;10(10):e0141042. https://doi.org/10.1371/journal.pone.0141042.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Vojtechova Z, Tachezy R. The Role of miRNAs in Virus-Mediated Oncogenesis. International Journal of Molecular Sciences 2018;19:1217. https://doi.org/10.3390/ijms19041217.

  26. Kolenda T, Kopczynska M, Guglas K, Teresiak A, Blizniak R, Lasinska I, et al. EGOT lncRNA in head and neck squamous cell carcinomas. Polish J Pathol. 2018;69(4):356–65. https://doi.org/10.5114/pjp.2018.81695.

    Article  Google Scholar 

  27. Sailer V, Charpentier A, Dietrich J, Vogt TJ, Franzen A, Bootz F, et al. Intragenic DNA methylation of PITX1 and the adjacent long non-coding RNA C5orf66-AS1 are prognostic biomarkers in patients with head and neck squamous cell carcinomas. 2018;13(2):e0192742.

    Google Scholar 

  28. Harden ME, Munger K. Human papillomavirus molecular biology. Mutation Res Rev Mutation Res. 2017;772:3–12. https://doi.org/10.1016/j.mrrev.2016.07.002.

    Article  CAS  Google Scholar 

  29. Harden ME, Prasad N, Griffiths A, Munger K. Modulation of microRNA-mRNA Target Pairs by Human Papillomavirus 16 Oncoproteins. mBio. 2017;8(1):e02170-16. https://doi.org/10.1128/mBio.02170-16.

  30. Sharma S, Hussain S, Soni K, Singhal P, Tripathi R, Ramachandran VG, et al. Novel MicroRNA signatures in HPV-mediated cervical carcinogenesis in Indian women. Tumour Biol. 2016;37(4):4585–95. https://doi.org/10.1007/s13277-015-4248-7.

    Article  CAS  PubMed  Google Scholar 

  31. Wang X, Wang HK, Li Y, Hafner M, Banerjee NS, Tang S, et al. microRNAs are biomarkers of oncogenic human papillomavirus infections. Proc Natl Acad Sci U S A. 2014;111(11):4262–7. https://doi.org/10.1073/pnas.1401430111.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Wilting SM, Snijders PJ, Verlaat W, Jaspers A, van de Wiel MA, van Wieringen WN, et al. Altered microRNA expression associated with chromosomal changes contributes to cervical carcinogenesis. Oncogene. 2013;32(1):106–16. https://doi.org/10.1038/onc.2012.20.

    Article  CAS  PubMed  Google Scholar 

  33. Yablonska S, Hoskins EE, Wells SI, Khan SA. Identification of miRNAs dysregulated in human foreskin keratinocytes (HFKs) expressing the human papillomavirus (HPV) Type 16 E6 and E7 oncoproteins. MicroRNA (Shariqah, United Arab Emirates). 2013;2(1):2–13.

    CAS  Google Scholar 

  34. Zheng ZM, Wang X. Regulation of cellular miRNA expression by human papillomaviruses. Biochim Biophys Acta. 2011;1809(11–12):668–77. https://doi.org/10.1016/j.bbagrm.2011.05.005.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Deftereos G, Corrie SR, Feng Q, Morihara J, Stern J, Hawes SE, et al. Expression of mir−21 and mir-143 in cervical specimens ranging from histologically normal through to invasive cervical cancer. PLoS One. 2011;6(12):e28423. https://doi.org/10.1371/journal.pone.0028423.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Fang J, Zhang H, Jin S. Epigenetics and cervical cancer: from pathogenesis to therapy. Tumour Biol. 2014;35(6):5083–93. https://doi.org/10.1007/s13277-014-1737-z.

    Article  CAS  PubMed  Google Scholar 

  37. Peta E, Sinigaglia A, Masi G, Di Camillo B, Grassi A, Trevisan M, et al. HPV16 E6 and E7 upregulate the histone lysine demethylase KDM2B through the c-MYC/miR-146a-5p axys. 2018;37(12):1654–1668.

    Google Scholar 

  38. Zamani S, Sohrabi A, Hosseini SM, Rahnamaye-Farzami M, Akbari A. Deregulation of miR−21 and miR-29a in Cervical Cancer Related to HPV Infection. MicroRNA (Shariqah, United Arab Emirates). 2019;8(2):110–5.

    CAS  Google Scholar 

  39. Yang L, Yi K, Wang H, Zhao Y, Xi M. Comprehensive analysis of lncRNAs microarray profile and mRNA-lncRNA co-expression in oncogenic HPV-positive cervical cancer cell lines. Oncotarget. 2016;7(31):49917–29. https://doi.org/10.18632/oncotarget.10232.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Zhao M, Qiu Y, Yang B, Sun L, Hei K, Du X, et al. Long non-coding RNAs involved in gynecological cancer. Int J Gynecol Cancer. 2014;24(7):1140–5. https://doi.org/10.1097/IGC.0000000000000212.

    Article  PubMed  Google Scholar 

  41. Barr JA, Hayes KE, Brownmiller T, Harold AD, Jagannathan R, Lockman PR, et al. Long non-coding RNA FAM83H-AS1 is regulated by human papillomavirus 16 E6 independently of p53 in cervical cancer cells. Sci Rep. 2019;9(1):3662. https://doi.org/10.1038/s41598-019-40094-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Sharma S, Munger K. Expression of the cervical carcinoma expressed PCNA regulatory (CCEPR) long noncoding RNA is driven by the human papillomavirus E6 protein and modulates cell proliferation independent of PCNA. Virology. 2018;518:8–13. https://doi.org/10.1016/j.virol.2018.01.031.

    Article  CAS  PubMed  Google Scholar 

  43. Iancu IV, Anton G, Botezatu A, Huica I, Nastase A, Socolov DG, et al. LINC01101 and LINC00277 expression levels as novel factors in HPV-induced cervical neoplasia. 2017;21(12):3787–3794.

    Google Scholar 

  44. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell. 2011;146(3):353–8. https://doi.org/10.1016/j.cell.2011.07.014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Braconi C, Kogure T, Valeri N, Huang N, Nuovo G, Costinean S, et al. microRNA-29 can regulate expression of the long non-coding RNA gene MEG3 in hepatocellular cancer. Oncogene. 2011;30(47):4750–6. https://doi.org/10.1038/onc.2011.193.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Liang WC, Fu WM, Wong CW, Wang Y, Wang WM, Hu GX, et al. The lncRNA H19 promotes epithelial to mesenchymal transition by functioning as miRNA sponges in colorectal cancer. Oncotarget. 2015;6(26):22513–25. https://doi.org/10.18632/oncotarget.4154.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Liu S, Song L, Zeng S, Zhang L. MALAT1-miR-124-RBG2 axis is involved in growth and invasion of HR-HPV-positive cervical cancer cells. Tumour Biol. 2016;37(1):633–40. https://doi.org/10.1007/s13277-015-3732-4.

    Article  CAS  PubMed  Google Scholar 

  48. Zhang J, Yao T, Wang Y, Yu J, Liu Y, Lin Z. Long noncoding RNA MEG3 is downregulated in cervical cancer and affects cell proliferation and apoptosis by regulating miR-21. Cancer biology & therapy. 2016;17(1):104–13. https://doi.org/10.1080/15384047.2015.1108496.

    Article  CAS  Google Scholar 

  49. Zhou X, Ji G, Ke X, Gu H, Jin W, Zhang G. MiR−141 inhibits gastric Cancer proliferation by interacting with long noncoding RNA MEG3 and Down-regulating E2F3 expression. Dig Dis Sci. 2015;60(11):3271–82. https://doi.org/10.1007/s10620-015-3782-x.

    Article  CAS  PubMed  Google Scholar 

  50. Gyongyosi E, Szalmas A, Ferenczi A, Konya J, Gergely L, Veress G. Effects of human papillomavirus (HPV) type 16 oncoproteins on the expression of involucrin in human keratinocytes. Virol J. 2012;9(1):36. https://doi.org/10.1186/1743-422X-9-36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Gyongyosi E, Szalmas A, Ferenczi A, Poliska S, Konya J, Veress G. Transcriptional regulation of genes involved in keratinocyte differentiation by human papillomavirus 16 oncoproteins. Arch Virol. 2015;160(2):389–98. https://doi.org/10.1007/s00705-014-2305-y.

    Article  CAS  PubMed  Google Scholar 

  52. Borbély AA, Murvai M, Kónya J, Beck Z, Gergely L, Li F, et al. Effects of human papillomavirus type 16 oncoproteins on survivin gene expression. J Gen Virol. 2006;87(Pt 2):287–94. https://doi.org/10.1099/vir.0.81067-0.

    Article  CAS  PubMed  Google Scholar 

  53. Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 2014;42(Database issue):D92–7. https://doi.org/10.1093/nar/gkt1248.

    Article  CAS  PubMed  Google Scholar 

  54. Paraskevopoulou MD, Vlachos IS, Karagkouni D, Georgakilas G, Kanellos I, Vergoulis T, et al. DIANA-LncBase v2: indexing microRNA targets on non-coding transcripts. Nucleic Acids Res. 2016;44(D1):D231–8. https://doi.org/10.1093/nar/gkv1270.

    Article  CAS  PubMed  Google Scholar 

  55. Chou CH, Chang NW, Shrestha S, Hsu SD, Lin YL, Lee WH, et al. miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 2016;44(D1):D239–47. https://doi.org/10.1093/nar/gkv1258.

    Article  CAS  PubMed  Google Scholar 

  56. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. https://doi.org/10.1101/gr.1239303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    Article  PubMed  Google Scholar 

  58. Szalmas A, Gyongyosi E, Ferenczi A, Laszlo B, Karosi T, Csomor P, et al. Activation of Src, Fyn and yes non-receptor tyrosine kinases in keratinocytes expressing human papillomavirus (HPV) type 16 E7 oncoprotein. Virol J. 2013;10(1):79. https://doi.org/10.1186/1743-422X-10-79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47(D1):D155–d62. https://doi.org/10.1093/nar/gky1141.

    Article  CAS  PubMed  Google Scholar 

  60. McLaughlin-Drubin ME, Munger K. Oncogenic activities of human papillomaviruses. Virus Res. 2009;143(2):195–208. https://doi.org/10.1016/j.virusres.2009.06.008.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Ou L, Wang D, Zhang H, Yu Q, Hua F. Decreased expression of miR−138-5p by lncRNA H19 in cervical Cancer promotes tumor proliferation. Oncol Res. 2018;26(3):401–10. https://doi.org/10.3727/096504017X15017209042610.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Graham SV, Faizo AAA. Control of human papillomavirus gene expression by alternative splicing. Virus Res. 2017;231:83–95. https://doi.org/10.1016/j.virusres.2016.11.016.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Li Y, Cai Q, Lin L, Xu C. MiR-875 and miR-3144 switch the human papillomavirus 16 E6/E6* mRNA ratio through the EGFR pathway and a direct targeting effect. Gene. 2018;679:389–97. https://doi.org/10.1016/j.gene.2018.09.015.

    Article  CAS  PubMed  Google Scholar 

  64. Vande Pol SB, Klingelhutz AJ. Papillomavirus E6 oncoproteins. Virology. 2013;445(1–2):115–37. https://doi.org/10.1016/j.virol.2013.04.026.

    Article  CAS  PubMed  Google Scholar 

  65. Benetatos L, Vartholomatos G, Hatzimichael E. MEG3 imprinted gene contribution in tumorigenesis. Int J Cancer. 2011;129(4):773–9. https://doi.org/10.1002/ijc.26052.

    Article  CAS  PubMed  Google Scholar 

  66. Baldassarre A, Masotti A. Long non-coding RNAs and p53 regulation. Int J Mol Sci. 2012;13(12):16708–17. https://doi.org/10.3390/ijms131216708.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

The study was supported by a grant from the Faculty of Medicine, University of Debrecen. AS was supported by the János Bolyai Research Scholarship awarded by the Hungarian Academy of Sciences.

Author information

Authors and Affiliations

  1. Department of Medical Microbiology, Faculty of Medicine, University of Debrecen, Nagyerdei krt. 98, Debrecen, H-4032, Hungary

    Brigitta László, Eszter Gyöngyösi, Anita Szalmás, György Veress & József Kónya

  2. Department of Hydrobiology, University of Debrecen, Egyetem tér 1, Debrecen, H-4032, Hungary

    László Antal

  3. Genomic Medicine and Bioinformatics Core Facility, Department of Biochemistry and Molecular Biology, Faculty of Medicine, University of Debrecen, Nagyerdei krt. 98, Debrecen, H-4032, Hungary

    Szilárd Póliska

Authors
  1. Brigitta László
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. László Antal
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Eszter Gyöngyösi
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Anita Szalmás
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Szilárd Póliska
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. György Veress
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. József Kónya
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

BL, JK and GV conceived and designed the study. BL, EG and SP carried out the experiments. BL, EG, SP and JK performed data analysis and interpretations. BL and LA prepared the first draft of the manuscript. BL, LA, EG, AS, GV and JK drafted and revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Brigitta László.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

The lists of StarBase v2.0, DianaTools-LncBase v.2 and miRTarBase search results.

Additional file 2.

Result of HPV16 E6 and E7 specific RT-PCR.

Additional file 3.

The effects of the HPV16 E6 and E7 on the level of well-known transcriptional targets of high-risk HPV oncoproteins.

Additional file 4.

Influence of individual HPV16 oncoproteins on miRNA and lncRNA expression changes.

Additional file 5.

The lists of genes significantly enriched GO and KEGG pathways.

Rights and permissions

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

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

László, B., Antal, L., Gyöngyösi, E. et al. Coordinated action of human papillomavirus type 16 E6 and E7 oncoproteins on competitive endogenous RNA (ceRNA) network members in primary human keratinocytes. BMC Cancer 21, 673 (2021). https://doi.org/10.1186/s12885-021-08361-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-021-08361-y

Keywords

BMC Cancer

ISSN: 1471-2407

Contact us