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

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.


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.

Background
Human papillomaviruses (HPVs) are small, doublestranded 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 posttranscriptionally 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.

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 C t 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 passagematched 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, serumfree 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 highrisk 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.
To validate the miRNA-seq results, altered expression of several miRNAs was confirmed using TaqMan realtime 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 Taq-Man 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 [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).

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).      (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 cancerassociated 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 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.
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 noncoding 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 noncoding 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 Fig. 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 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.