- Research article
- Open Access
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)
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.
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.
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.
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.
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 . 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.
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 . 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 . 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 .
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 . 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 . 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 .
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 . 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).
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.
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 . 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).
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).
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.
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).
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).
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.
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 . 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 . 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 . 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 .
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.
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.
long non-coding RNA
competitive endogenous RNA
microRNA response element
Human foreskin keratinocyte
Kyoto Encyclopedia of Genes and Genomes
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.
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.
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.
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.
Moody CA, Laimins LA. Human papillomavirus oncoproteins: pathways to transformation. Nat Rev Cancer. 2010;10(8):550–60. https://doi.org/10.1038/nrc2886.
Mattick JS. The genetic signatures of noncoding RNAs. PLoS Genet. 2009;5(4):e1000459. https://doi.org/10.1371/journal.pgen.1000459.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The lists of StarBase v2.0, DianaTools-LncBase v.2 and miRTarBase search results.
Result of HPV16 E6 and E7 specific RT-PCR.
The effects of the HPV16 E6 and E7 on the level of well-known transcriptional targets of high-risk HPV oncoproteins.
Influence of individual HPV16 oncoproteins on miRNA and lncRNA expression changes.
The lists of genes significantly enriched GO and KEGG pathways.
About this article
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
- Long non-coding RNA
- Competing endogenous RNA (ceRNA)