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

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.

Keywords: Competing endogenous RNA (ceRNA); E6; E7; HPV16; Long non-coding RNA; microRNA.

Fig. 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

Fig. 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 C_t 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

Fig. 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

Fig. 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 C_t 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

Fig. 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

Fig. 6

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

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