Identification of Deregulated Pathways, Key Regulators, and Novel miRNA-mRNA Interactions in HPV-Mediated Transformation

Next to a persistent infection with high-risk human papillomavirus (HPV), molecular changes are required for the development of cervical cancer. To identify which molecular alterations drive carcinogenesis, we performed a comprehensive and longitudinal molecular characterization of HPV-transformed keratinocyte cell lines. Comparative genomic hybridization, mRNA, and miRNA expression analysis of four HPV-containing keratinocyte cell lines at eight different time points was performed. Data was analyzed using unsupervised hierarchical clustering, integrated longitudinal expression analysis, and pathway enrichment analysis. Biological relevance of identified key regulatory genes was evaluated in vitro and dual-luciferase assays were used to confirm predicted miRNA-mRNA interactions. We show that the acquisition of anchorage independence of HPV-containing keratinocyte cell lines is particularly associated with copy number alterations. Approximately one third of differentially expressed mRNAs and miRNAs was directly attributable to copy number alterations. Focal adhesion, TGF-beta signaling, and mTOR signaling pathways were enriched among these genes. PITX2 was identified as key regulator of TGF-beta signaling and inhibited cell growth in vitro, most likely by inducing cell cycle arrest and apoptosis. Predicted miRNA-mRNA interactions miR-221-3p_BRWD3, miR-221-3p_FOS, and miR-138-5p_PLXNB2 were confirmed in vitro. Integrated longitudinal analysis of our HPV-induced carcinogenesis model pinpointed relevant interconnected molecular changes and crucial signaling pathways in HPV-mediated transformation.


Introduction
Development of cervical cancer is a multi-step process initiated by a persistent infection with a high-risk type of the human papillomavirus (HPV) [1]. Following infection of the basal epithelial cells by HPV, a productive infection is established, which is characterized by the production of new viral particles [2,3]. The expression of viral proteins is tightly linked to host cell differentiation and is necessary for viral genome replication in differentiated cells. Aberrant expression of viral oncogenes E6 and E7 in proliferating cells, however, leads to the abortion of the viral life cycle and triggers malignant transformation. While E6 and E7 initiate and maintain transforming infections, additional molecular changes in the host cell genome are required for the development of invasive cancer.
Molecular alterations associated with cervical carcinogenesis are of both genetic and epigenetic nature and ultimately lead to aberrant expression of oncoproteins and tumor suppressors. Genome-wide analyzes of cervical tissue specimens have led to the identification of numerous chromosomal aberrations and differentially expressed coding and non-coding genes in cervical (pre)cancer [4][5][6]. As the sequential order of molecular changes as well as their causative relevance for cancer development cannot be extrapolated from cross-sectional data, it has proven problematic to distil promising disease markers and potential therapeutic targets from these observations.
In vitro cell line models of HPV-mediated transformation, in which primary keratinocytes are transfected with HPV16 or HPV18, have been shown to faithfully mimic cervical cancer development [7][8][9]. This offers the unique opportunity to study (epi)genetic alterations during carcinogenesis. Integrative analysis of longitudinal data obtained from multiple molecular levels allows for reconstruction of the sequential order in which molecular changes occur and is likely to result in the identification of crucial molecular alterations that drive the carcinogenic process. Importantly, HPV-transformed keratinocytes have previously been shown to provide useful model systems to study chromosomal instability in cancer in general [10], indicating that insights in HPV-induced transformation might not only provide a better understanding of cervical cancer development, but could be applicable to other cancers, too.
Previous studies have demonstrated that HPV-induced transformation can be divided into four stages ( Figure 1a) [11,12]: an extended lifespan (1) is acquired as a result of E6 and E7 mediated inhibition of tumor suppressor genes TP53 and RB1 [13]. Genetic instability induced by E6 and E7 subsequently generates an immortal phenotype (2) usually associated with the activation of the telomerase enzyme resulting from upregulation of TERT. During prolonged culturing of HPV-immortalized keratinocytes additional (epi)genetic alterations accumulate that can eventually lead to anchorage-independent growth (3), considered as proof for complete transformation in vitro [14,15], and tumorigenicity (4).
Here, we present the first comprehensive molecular profiling of HPV-induced carcinogenesis over time. DNA copy number changes, mRNA, and miRNA expression were determined in four individual HPV-transformed keratinocyte cell lines at eight different passages representing different stages of HPV-induced transformation (hereafter referred to as time points). Integrative temporal analysis of this unique longitudinal multi-level dataset allowed for the identification of relevant pathways and associated key regulators as well as the prediction of miRNA-mRNA interactions in HPV-induced transformation. The obtained results will greatly accelerate the identification of those molecular alterations that are crucial to HPV-induced carcinogenesis and as such will likely provide highly specific and sensitive disease markers. These are urgently needed to improve current cervical cancer screening programs for the coming era, especially with the implementation of HPV-based screening and HPV vaccination in more and more countries.  [16]. MiRNA microarrays of cell line FK18A at T4 and T5 did not pass quality control and were therefore excluded from miRNA analysis. Unsupervised hierarchical cluster results based on (b) DNA copy number, (c) overall mRNA expression, and (d) overall miRNA [16] expression are shown for FK16A, FK16B, FK18A, and FK18B.

Time Point
To obtain an unbiased overview of the observed molecular changes over time in the four cell lines, we performed unsupervised clustering analysis of genomic as well as mRNA and miRNA expression profiles. Chromosomal profiles of all four cell lines were separated almost perfectly based on their ability to grow anchorage-independently (Figure 1b), which supports our previous observation based on miRNA expression profiles ( Figure 1d) [16]. Revealing a slightly less clear-cut separation, a sub-cluster enrichment for either anchorage-dependent or anchorage-independent time points was still noticeable for mRNA profiles (Figure 1c).

Approximately One Third of Differentially Expressed Genes is Associated with Copy Number Changes
To identify genes relevant to HPV-induced carcinogenesis, we performed integrative longitudinal differential gene expression analysis, which also takes the genomic copy number into account [17]. In total, expression of 3642 mRNA genes and 106 mature miRNAs (corresponding to 118 miRNA genes) was found to either increase or decrease consistently over time in at least three out of four cell lines analyzed (Tables S1 and S2). An association with copy number was observed in 33.9% of significant mRNAs (1233 genes, Table S1) and 32.2% of miRNA genes (38 genes, Table S2).
A technical validation of our microarray results was performed by qRT-PCR for three mRNAs (DEK, DKK3, SLC25A36, Figure S1a) and five miRNAs (miR-100-5p, miR-103a-3p, miR-125b-5p, miR-15b-5p, miR-221-5p, Figure S1b). Moreover, in vivo relevance of our findings was investigated for those differentially expressed mRNAs and miRNAs for which in-house microarray data of cervical tissue samples was available (mRNAs: 2661 out of 3642, miRNA: 39 out of 106). A concordant significant pattern of up-or downregulation was observed from normal to CIN to SCC by Spearman correlation for 24% of mRNAs and 49% of miRNAs, whereas 15% of mRNAs and 18% of miRNAs exhibited an opposite expression pattern in tissues and the cell line model (Table S1 and S2).

Pinpointing Key Regulators in Enriched Pathways
Based on the unsupervised hierarchical cluster results we hypothesized that a substantial proportion of the above described DNA copy number induced gene expression changes are involved in the acquisition of anchorage independence. Indeed, three of the most significantly overrepresented pathways among all DNA copy number associated genes (Table S3), namely focal adhesion (KEGG hsa04510), TGF-beta (KEGG hsa04350), and mTOR signaling (KEGG hsa04150), are implicated in processes underlying anchorage independence (i.e., anoikis resistance and induction of EMT) [18]. To further investigate this notion, we translated our longitudinal data into pathway-based networks for focal adhesion, TGF-beta signaling, and mTOR signaling (Table S3 and Materials and Methods section). Data-driven longitudinal networks were built for all genes in the respective pathways taking copy number into account and key regulators were identified (Table 2, Figure 2, Table S4).

TGF-Beta Pathway
TGFB1 has previously been shown to act as tumor suppressor during early stages of cervical carcinogenesis, although it takes on a tumor promoting function at a later stage [19]. Here, we identified PITX2 as key regulator of TGF-beta signaling (Table 2c). PITX2 was downregulated in all four keratinocyte cell lines ( Figure 3a) and this was concordant with our in-house tissue data (normal to CIN2/3 to SCC, Figure 3b).
To investigate the functional relevance of PITX2, a PITX2 expression vector was transduced in late FK18B cells. Cell growth analysis demonstrated a decrease in cell viability after three days of PITX2 overexpression compared to empty vector (GFP only) transduced FK18B cells (Figure 3c,d). PITX2 has previously been shown to cause G0/G1 arrest in HPV18-positive HeLa cells through accumulation of TP53 and CDKN1A (p21) [20,21]. In line with this, we found elevated TP53 and CDKN1A protein levels in PITX2 overexpressing FK18B cells (Figure 3d). To elucidate whether the observed growth reduction in PITX2 overexpressing cells is caused by TP53/CDKN1A-mediated cell cycle arrest or apoptosis, we further analyzed known marker proteins CDKN1B (p27), PCNA, and cleaved CASP3. Even though there was no change in PCNA levels upon PITX2 overexpression, increased levels of CDKN1B and cleaved CASP3 suggested that reduced cell viability in PITX2 overexpressing cells is associated to both cell cycle arrest and apoptosis ( Figure 3d). (a) PITX2 expression in HPV-transformed keratinocyte cell lines. Using our tigaR method, a model without (red continuous line) and with copy number effect (blue discontinuous line) was fitted [17]. Similarity between the fitted models indicates that differential PITX2 expression is not copy number-driven.

Identification of Potential miRNA-mRNA Target Interactions
Using our developed tigaR framework we next investigated the relation between miRNA and mRNA expression in a longitudinal fashion. For this analysis, we only included miRNAs for which altered expression was not contradictory to the pattern observed in the tissue specimens (n = 86). When restricting the analysis to miRNA-mRNA pairs showing significant expression changes over time in the opposite direction, 634 interactions between 27 miRNAs and 372 mRNAs were identified (Table S5).
To validate our miRNA target prediction, the four miRNA-mRNA interactions that were predicted by all three databases, i.e., miR-138-5p_PLXNB2, miR-221-3p_BRWD3, miR-221-3p_FOS, and miR-30a-3p_PECR, were selected. For all four miRNA-mRNA pairs, mRNA expression increased and miRNA expression decreased with increasing passage number in at least three out of four cell lines. Ectopic overexpression of the respective miRNA in late passage FK18B cells resulted in reduced expression levels of the predicted mRNA targets BRWD3, PECR, and PLXNB2 (Figure 4a-d). Ectopic overexpression of miR-221-3p surprisingly led to 40% increase of FOS expression in FK18B. Dual-luciferase reporter assay to confirm direct interaction between the miRNA and its predicted target mRNA. HEK293 cells were transiently transfected with a negative control (cntr#1) or the respective miRNA mimic in combination with either an empty pmiRGLO vector, a pmiRGLO construct containing the predicted binding site (pmiRGLO-predicted target-UTR), or a pmiRGLO construct containing a mutated binding site (pmiRGLO-predicted target-UTR_mut). MiRNA, UTR, and binding sites are indicated in the lower panel. Mean and standard deviation of three technical replicates are shown. * p < 0.05, ** p < 0.005., according to the two-sided Student's t-test.
Next, dual-luciferase assays were performed to investigate whether the effects of miRNAs on their target mRNAs are the result of direct interactions. As shown in Figure 4a, co-transfection of pmiRGLO-BRWD3-UTR and ectopic miR-221-3p in HEK293 cells decreased luciferase activity compared to either co-transfection of pmiRGLO-BRWD3-UTR with a non-targeting control sequence or the empty pmiRGLO vector with ectopic miR-221-3p. Directed mutagenesis of the predicted miR-221-3p seed sequence (pmiRGLO-BRWD3-UTR_mut) abolished the reduction in luciferase activity observed with the wild type vector (Figure 4a). Likewise, a direct interaction between miR-221-3p and the cFOS 3'UTR could be demonstrated (Figure 4b). Luciferase activity of a pmiRGLO-PECR-UTR construct was comparable to that of the empty pmiRGLO vector in the presence of ectopic miR-30a-3p or a non-targeting control sequence, suggesting that the observed reduction in PECR mRNA levels upon ectopic overexpression of miR-30a-3p is an indirect effect (Figure 4c). Interestingly, PLXNB2 has two predicted miR-138-5p binding sites in the 3'UTR. Luciferase activity of the pmiRGLO-PLXNB2-UTR vector was reduced in the presence of miR-138-5p compared to the empty vector (Figure 4d). Mutation of the first binding site significantly restored luciferase activity, while mutation of the second binding site alone had almost no effect. The largest and most significant restoration of luciferase activity was obtained when both binding sites were mutated.
Interestingly, ectopic expression of miR-30a-3p and miR-138-5p decreased anchorage-independent growth of late passage FK18B cells as measured on ultra-low attachment plates ( Figure 5). No reduction in cell viability was observed on normal, adherent plates, suggesting involvement of these miRNAs specifically in anchorage independence.

Discussion
We here present one of the most extensive studies characterizing HPV-induced transformation in a longitudinal fashion. Unsupervised hierarchical clustering of our copy number data separated nearly perfectly between anchorage-dependent and -independent passages of HPV-transformed cells, demonstrating the importance of chromosomal alterations in the acquisition of anchorage independence. Moreover, we show that differential expression of a considerable fraction of mRNAs and miRNAs is directly attributable to copy number changes of their respective genes. Pathway enrichment analysis on copy number affected mRNAs identified focal adhesion, mTOR signaling, and TGF-beta signaling as altered pathways in HPV-induced carcinogenesis. The biological significance of these findings is supported by functional evidence for PITX2, a key regulator of TGF-beta signaling, and novel miRNA-mRNA interactions. Anchorage-independent cell growth is a hallmark of transformation in vitro [14,15]. While loss of appropriate cell-cell and cell-matrix interactions leads to aberrant integrin signaling and eventually detachment-induced cell death (anoikis) in healthy epithelial cells, cancerous epithelial cells overcome anoikis by epithelial-to-mesenchymal transition (EMT), activation of survival and proliferation pathways, or temporary dormancy [25]. Strikingly, the focal adhesion, mTOR signaling, and TGF-beta signaling pathways identified in present study are implicated in anoikis resistance and induction of EMT [18,26].
TGFB1 has been reported to exert a tumor inhibiting function in early stages of cervical carcinogenesis, while it plays an oncogenic role at a later stage by promotion of EMT and metastasis, induction of angiogenesis, and escape from immune surveillance [19]. In line with this, expression of TGFB1 was found to be downregulated in cervical precancerous lesions compared to normal epithelium although it is increased in cervical cancers [27,28]. This TGF-beta paradox is not unique to cervical carcinogenesis but has been recognized in many other cancers as well [29]. An analysis of 178 cervical carcinomas by The Cancer Genome Atlas (TCGA) network found recurrent deletions of TGF-beta receptor 2 (TGFBR2, 3p24.1) and SMAD4 (18q21.2) in 36% and 28% of analyzed samples, respectively [30]. In accordance, downregulation of both TGFBR2 and SMAD4 was related to copy number loss in our cell line data (Table S1). Moreover, the TCGA network identified common somatic mutations in the TGF-beta pathway for TGFBR2, SMAD4, CREBBP, and EP300 [30]. These mutations are of an inactivating nature and are believed to deactivate growth-suppressive and pro-apoptotic functions of TGF-beta, leading to a more tumorigenic phenotype.
In a previous study, we reported that the presence of TGFB1 reduced the growth rates of mortal and early immortal FK16A and FK18B cells, but not of late immortal FK16A and FK18B cells, in a SMAD2/SMAD4-independent fashion [31]. Similarly, Creek and colleagues showed that HPV16-transformed cells become insensitive to TGFB1 with increasing passaging and demonstrated that TGFB1 promotes EMT and cell migration [32][33][34]. A recent integrative screen comparing transcriptome and proteome data of HPV16 E6E7 expressing keratinocytes to non-HPV expressing keratinocytes confirmed that TGFB1 plays a regulatory role in HPV-induced carcinogenesis [35]. Here, our extensive longitudinal pathway analysis identified PITX2 (pituitary homeobox 2) as key regulator of TGF-beta signaling during cervical carcinogenesis. We show that elevated PITX2 expression induces apoptosis and cell cycle arrest, indicating that downregulation of PITX2 is important for the survival of HPV-transformed keratinocytes. In line with our results, PITX2 has been shown to cause cell cycle arrest through accumulation of TP53 and CDKN1A in HeLa cells [20,21]. This is achieved through binding of the HPV E6 protein by PITX2 and the subsequent inhibition of TP53 degradation [20]. Moreover, PITX2 has been described as transcriptional activator of CDKN1A in other tissues such as the dental epithelium and neural stem cells [36,37]. Besides an effect on cell growth, others have observed decreased migration upon PITX2 overexpression [20]. Wei and Adelstein showed that PITX2 activates the Rho GTPases RAC1 and RHOA, causing changes to the actin-myosin cytoskeleton thereby reducing cell adhesion and increasing cell motility [20]. Interestingly, our pathway reconstruction also identified RHOA as PITX2 target, as well as ROCK2 and THBS1, two other genes known to be involved in cell adhesion and motility. Future studies on the involvement of PITX2 in growth inhibition, cell cycle, migration, and invasion during HPV-induced transformation are therefore warranted.
The tumor suppressive role of PITX2 in cervical cancer development is further supported by our previously published whole methylome data showing that the PITX2 promoter is highly methylated in high passages of our four cell lines but not in primary donor keratinocytes [16]. PITX2 hypermethylation has been described as a promising prognostic biomarker in cervical cancer as well as breast cancer, head and neck squamous cell carcinoma, and prostate cancer [38][39][40][41].
Besides studying the effect of copy number aberrations on gene expression, the tigaR framework also proved suitable for the identification of temporal miRNA-mRNA interactions [17]. Opposed to BRWD3 and PLXNB2 being downregulated upon overexpression of their validated miRNA interaction partner, FOS expression increased upon miR-221-3p overexpression. This observation is potentially caused by the high ectopic levels of miR-221-3p in our experiments, which may result in an indirect effect of miR-221-3p on FOS expression via other miR-221-3p targets. Under physiological conditions, we do not observe this indirect effect as illustrated by the opposing expression directions of FOS (upregulated) and miR-221-3p (downregulated) over time in our cell line model as well as during progression in cervical tissue specimens. In our cell line model, all investigated target mRNAs were upregulated with increasing passages, suggesting oncogenic roles for these genes. While little is known about the functional involvement of BRWD3 and PECR in human cancers, an enrichment of FOS as part of a heterodimer with JUN has previously been linked to anchorage-independent cell growth in our HPV-transformed cell lines as well as HeLa [42,43]. Moreover, PLXNB2 has been described to act as oncogene in ovarian cancer cells, where a knock-down of PLXNB2 led to decreased cell viability and invasion [44]. We demonstrate that miR-30a-3p and miR-138-5p reduce anchorage-independent cell growth of late passage FK18B cells, indicating their functional role in anchorage independence. Additional studies investigating the exact functional roles of the here identified miRNA-mRNA pairs in cervical carcinogenesis might provide valuable therapeutic targets in the future.

Cell Lines and Clinical Specimens
Establishment and culture of the HPV16-(FK16A and FK16B) and HPV18-(FK18A and FK18B) immortalized keratinocyte cell lines has been described previously [7]. Renal epithelium cell line HEK293 was authenticated by STR testing using the Powerplex16 System (Promega, Leiden, The Netherlands) and cultured as described previously [45]. From all four HPV-immortalized keratinocyte cell lines eight passages (Table 1), including both anchorage-dependent (no shading) and anchorage-independent (grey shading) cells, were selected.

RNA and DNA Isolation
Total RNA was isolated using TRIzol Reagent according to the manufacturer's instructions (Thermo Fisher Scientific, Bleiswijk, The Netherlands). RNA integrity was determined by gel electrophoresis. Total DNA was isolated by standard proteinase K digestion followed by phenol-chloroform purification [46].

mRNA Arrays
Global mRNA expression profiles were generated using whole human genome oligo microarrays (G4112F, mRNA 4 × 44K; Agilent Technologies) following the manufacturer's instructions. High-resolution mRNA expression data are available from GEO through series accession number GSE138079.

miRNA Arrays
Global miRNA expression profiles were determined using human miRNA microarrays (Sureprint G3 human v16 miRNA 8 × 60K; Agilent Technologies) according to the manufacturer's instructions. These arrays contain in situ synthesized 60-mer oligonucleotides representing 1205 human miRNAs based on the Sanger miRBase release 16. Microarrays for two passages of FK18A (passage 60 and 92, i.e., T4 and T5) failed subsequent quality control and were excluded from analysis. Microarray data have previously been published and are available from GEO through series accession number GSE78279 [16].

Expression Profiling on Tissue Specimens
To confirm in vivo relevance of our cell line data, miRNA, and mRNA expression was analyzed in microarray data obtained from cervical tissue specimens. Global miRNA profiles were obtained from normal HPV-positive cervical epithelium (n = 10), high-grade precancerous lesions (CIN2/3, n = 18), and squamous cell carcinomas (SCC, n = 10). Data is available from GEO through series accession number GSE30656 [47]. The same tissue specimens were used to determine global mRNA expression profiles using whole human genome oligo microarrays (G4112F, mRNA 4 × 44K; Agilent Technologies). Three CIN2/3 did not pass quality control and were excluded from analysis. High-resolution mRNA expression data obtained from cervical tissue samples are available from GEO through series accession number GSE138080. Tissue samples were used in an anonymous fashion in accordance with the Code for Proper Secondary Use of Human Tissues in the Netherlands as formulated by the Dutch Federation of Medical Scientific Organizations (www.federa.org).

Pre-Processing and Matching
DNA copy number data were pre-processed employing median normalization and segmented using the circular binary segmentation method [48]. mRNA gene expression data were background corrected and between array normalization was performed using the robust quantile method [49,50]. Resulting gene expression intensities were transformed using variance stabilizing transformation [51]. MiRNA gene expression data were processed similarly but without background correction.
Prior to integrative downstream analysis, mRNA and miRNA gene transcripts were matched to DNA copy number data based on the genes' respective chromosomal locations (genomic build GRCh37/hg19) using the overlapPlus matching procedure described in Van Wieringen et al. [52].

Cluster Analysis
Unsupervised hierarchical clustering was performed per cell line on complete mRNA/miRNA expression profiles to explore the overall similarities and differences in expression patterns using maximum as distance measure. The R-package WECCA was used for unsupervised hierarchical clustering of samples based on their DNA copy number profiles per cell line [53]. All dendrograms were constructed using Ward's linkage as it yields compact and well-separated clusters.

Differential Expression Analysis
To identify m(i)RNA genes with temporal differential expression and to study the association between gene expression and underlying copy number changes, the methodology presented in Miok et al. was employed [17]. Temporal gene expression analysis was performed with a common spline for the cell lines, to identify consistently altered genes over time at a 5% false-discovery rate (FDR). The same methodology was used to investigate associations between mRNA and miRNA gene expression data over time [17]. From this analysis, significantly associated miRNA-mRNA gene pairs with a negative regression parameter were selected. Identified miRNA-mRNA interactions were verified by three independent databases: RNA22 v2.0, miRDB v5, and TargetScan v7 [22][23][24]. RNA22 v2.0 uses conserved sequence features among miRNAs, so-called patterns, to identify binding site hotspots in the 3' untranslated region (3'UTR) of mRNAs [22]. The identified putative binding sites are then associated to individual miRNAs independent of their seed sequences. miRDB and TargetScan are both based on target prediction features such as miRNA seed conservation or the mRNA sequence flanking the binding site [23,24]. While the MirTarget2 algorithm employed by miRDB uses a support vector machine framework, TargetScan employs a linear regression model. In both cases, prediction features were selected and weighted based on large experimental datasets.

Network Modeling
For network enrichment analysis, 118 potentially relevant pathways stored in the KEGG repository were considered (Table S1) [54]. Enrichment analysis showed significant overrepresentation of 21 pathways within the copy number related differentially expressed mRNAs (chi-square test with FDR correction). From these 21 pathways, focal adhesion (KEGG hsa04510), TGF-beta signaling (KEGG hsa04350), and mTOR signaling (KEGG hsa04150) were selected for subsequent pathway reconstruction analysis based on the following criteria: (1) total number of genes in pathway is between 50 and 200 (requirement for downstream pathway reconstruction), (2) more than 5% of genes within the pathway showed copy number induced differential expression. Although the B-cell receptor signaling pathway (hsa04662) and cell cycle (hsa04110) also met these criteria, they were not investigated further. This analysis aimed to identify hub genes, which represent possible key regulators of the pathway in the given context. For this purpose, mRNA gene expression data from the experiment was mapped to the pathways as defined by the KEGG repository and the vector autoregressive model with time-varying DNA copy number was employed [54,55]. This model allows to identify both temporal (signals propagated over time) and contemporaneous (signal propagated among the genes) interactions among the genes. This integrative methodology allows identification of the 'within-level' relationship in the mRNA gene levels as well as relations between different molecular levels [55,56].

mRNA qRT-PCR
To determine expression levels of BRWD3, DEK, DKK3, FOS, PECR, PLXNB2, SNRPA, and SLC25A36, cDNA was synthesized from 200 ng total RNA template using specific reverse primers (Table S6a). All primer pairs were intron-flanking with an exception for mono-exonic FOS. For analysis of FOS expression, RNA samples were therefore treated with RQ1 DNase (Promega) prior to reverse transcription. The resulting cDNA was used for SYBR Green qPCR on the ABI7500 Fast Real-Time PCR System (Thermo Fisher Scientific). Specificity of the PCR reaction was determined generating melting curves for each reaction. Each 25 µL PCR reaction contained 12.5 µL 2× SYBR Green master mix (Thermo Fisher Scientific), 0.5 µM forward and reverse primers, and 2.5 µL cDNA. Cycle conditions used were according to the manufacturer's instructions. Data were normalized so SNRPA using the 2 −∆Ct method [57].

Viral Transduction and Flow Cytometry
The coding sequence of PITX2 transcript variant 2, amplified by Phusion High-Fidelity PCR (New England Biolabs (NEB), Ipswich, MA, USA), was cloned into the lentiviral vector LeGO-iG2 [58] (Addgene, Teddington, UK) via the pcDNA3.1(+) Mammalian Expression Vector (Thermo Fisher Scientific) using HindIII and NotI restriction sites. For primer sequences please refer to Table S6b. Restriction enzymes were purchased from NEB. For virus production, HEK293 cells were transfected with lentiviral packaging plasmids pRSV-rev, pCMV-VSV-G, and pMDLg/RRE together with LeGO-iG2 empty vector or LeGO-iG2-PITX2 using polyethylenimine (PEI; Sigma-Aldrich, St. Louis, MO, USA). Late passage FK18B cells (ca. passage 190) were transduced with lentivirus using polybrene (Thermo Fisher Scientific) for 24 h. Cells were harvested for RNA 72 h after transduction. Analysis of GFP expression was performed on a FACS-Calibur flow cytometer (BD Biosciences, Breda, The Netherlands) equipped with the CellQuest Pro data acquisition and analysis software. Two independent experiments were performed.
The predicted 3'UTR binding sites of BRWD3, FOS, PECR, and PLXNB2 amplified by Phusion High-Fidelity PCR (NEB) were cloned into the pmirGLO Dual-Luciferase miRNA Target Expression Vector (Promega) using the SacI (NEB) and XhoI (NEB) restriction sites. Predicted miRNA binding sites were mutated using the Q5 Site-Directed Mutagenesis Kit (NEB) according to the manufacturer's instructions (Figure 4a-d). Primer sequences are listed in Table S6c-d. HEK293 were seeded in triplicate in 96-well plates (7500 cells/well) and co-transfected with 6.67 nM mimics and 20 ng pmiRGLO vector the following day. Firefly luciferase activity was measured using the Dual-Glo Luciferase Assay (Promega) according to the manufacturer's instructions 48 h after transfection. Renilla luciferase activity was determined as internal control and used for normalization. Each experiment was carried out three times.

Conclusions
In conclusion, we present a unique longitudinal multi-level dataset on HPV-induced transformation. Using state-of-the-art analysis tools, we identified key pathway regulators and novel miRNA-mRNA interactions relevant to the carcinogenic process. Our results highlight the importance of chromosomal alterations in the acquisition of anchorage independence during carcinogenesis. Additional functional studies on the identified key pathway regulators and differentially expressed m(i)RNAs might yield potential therapeutic targets and disease markers in the future.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-6694/12/3/700/s1, Figure S1: Technical validation of microarray results by qRT-PCR, Table S1: mRNAs that consistently change over time in at least three out of four cell lines, Table S2: miRNAs that consistently change over time in at least three out of four cell lines, Table S3: Network enrichment analysis of copy number related genes in 118 KEGG pathways, Table S4: Key regulators of selected pathways and their interactions, Table S5: Predicted mRNA-miRNA interactions, Table S6: Primer sequences.