A systems biology analysis of the changes in gene expression via silencing of HPV-18 E1 expression in HeLa cells

Previous studies have reported the detection of a truncated E1 mRNA generated from HPV-18 in HeLa cells. Although it is unclear whether a truncated E1 protein could function as a replicative helicase for viral replication, it would still retain binding sites for potential interactions with different host cell proteins. Furthermore, in this study, we found evidence in support of expression of full-length HPV-18 E1 mRNA in HeLa cells. To determine whether interactions between E1 and cellular proteins play an important role in cellular processes other than viral replication, genome-wide expression profiles of HPV-18 positive HeLa cells were compared before and after the siRNA knockdown of E1 expression. Differential expression and gene set enrichment analysis uncovered four functionally related sets of genes implicated in host defence mechanisms against viral infection. These included the toll-like receptor, interferon and apoptosis pathways, along with the antiviral interferon-stimulated gene set. In addition, we found that the transcriptional coactivator E1A-binding protein p300 (EP300) was downregulated, which is interesting given that EP300 is thought to be required for the transcription of HPV-18 genes in HeLa cells. The observed changes in gene expression produced via the silencing of HPV-18 E1 expression in HeLa cells indicate that in addition to its well-known role in viral replication, the E1 protein may also play an important role in mitigating the host's ability to defend against viral infection.


Summary
Previous studies have reported the detection of a truncated E1 mRNA generated from HPV-18 in HeLa cells. Although it is unclear whether a truncated E1 protein could function as a replicative helicase for viral replication, it would still retain binding sites for potential interactions with different host cell proteins. Furthermore, in this study, we found evidence in support of expression of full-length HPV-18 E1 mRNA in HeLa cells. To determine whether interactions between E1 and cellular proteins play an important role in cellular processes other than viral replication, genome-wide expression profiles of HPV-18 positive HeLa cells were compared before and after the siRNA knockdown of E1 expression. Differential expression and gene set enrichment analysis uncovered four functionally related sets of genes implicated in host defence mechanisms against viral infection. These included the tolllike receptor, interferon and apoptosis pathways, along with the antiviral interferon-stimulated gene set. In addition, we found that the transcriptional coactivator E1A-binding protein p300 (EP300) was downregulated, which is interesting given that EP300 is thought to be required for the transcription of HPV-18 genes in HeLa cells. The observed changes in gene expression produced via the silencing of HPV-18 E1 expression in HeLa cells indicate that in addition to its well-known role in viral replication, the E1 protein may also play an important role in mitigating the host's ability to defend against viral infection.

Introduction
The International Agency for Research on Cancer has indicated that there is convincing evidence that infection with human papilloma virus type 16 (HPV-16) and HPV-18 can result in cervical cancer, as well as in cancers of the vulva, vagina, penis and anus [1]. HPV DNA has also been found in cancers of the oral cavity, oropharynx, larynx, oesophagus and lung, and thus its association with cancers of the upper aero-digestive tract has also been suspected [2][3][4][5]. In addition, HPV DNA has been detected in cell lines derived from cervical cancers: HPV-16 DNA in the cell lines CaSki and SiHa, and HPV-18 DNA in HeLa cells [6][7][8][9]. Owing to the difficulty of establishing tissue culture systems, which are susceptible to transformation by HPV, these cell lines provide unique systems to study the expression of HPV- 16  In HeLa cells, a cervical adenocarcinoma-derived cell line containing multiple copies of integrated HPV-18 DNA [10][11][12][13][14], HPV-18 early mRNAs are transcribed as polycistronic RNA [10], which give rise to three differentially spliced mRNA species [11,12]. These early mRNAs contain information for the translation of three potential HPV-18 proteins: E1, E6 and E7. Seedorf et al. [11] detected two early proteins-E7 (12 kDa) and a truncated E1 (70 kDa)-in HeLa cells, as assessed using polyacrylamide gels, following hybrid selection and Western blotting analyses, in which the sizes were consistent with the sizes predicted from the HPV-18 DNA sequence.
Abundant information is available on the E7 protein and its role in both the virus infection cycle and carcinogenesis. The viral protein E7 can induce cell immortalization by binding to pRb, a tumour suppressor protein, which binds to and inactivates the E2F transcription factor. E2F is released from pRb, which results in the transcription of genes involved in DNA replication and cell division [15]. In addition, E1 is a replicative helicase protein, which binds to the host DNA polymerase alpha-primase [16], replication protein A [17] and topoisomerase I to perform viral replication [18]. Moreover, E1 also interacts with the molecular chaperones, Hsp40 and Hsp70 [19], cellular WD-repeat protein 80, which maintains the viral genome in keratinocytes [20], E/cdk2 cyclin complex [21], and with cellular proteins involved in chromatin remodelling and co-transcription activation, such as histone H1 [22] and Ini1/hSNF5, a subunit of the SWI/SNF chromatin remodelling complex [23]. However, it is still unknown whether these interactions play a role in cellular processes other than viral replication.
To determine if these potential interactions play an important role in cellular processes other than viral replication, HPV-18 E1 mRNA, which is endogenously expressed in HeLa cells, was silenced using short interfering RNA sequences (siRNAs). Subsequently, we examined changes in the expression of cellular genes using microarray analysis. Differentially expressed genes were analysed using gene set enrichment in order to look for functionally related sets of genes that may be coordinately changed upon the E1 mRNA knockdown.

Short interfering RNA sequences
Six siRNAs targeting HPV-18 E1 mRNAs (siE1.1 -6) were selected using an algorithm based on guidelines developed by Ui-Tei et al. [24] and the web-based online software system SIDIRECT [25]. The siRNA molecule sequences were designed under the following conditions: A/U at the 5 0 end of the antisense strand; G/C at the 5 0 end of the sense strand; at least five A/U residues in one-third from the 5 0 end of the antisense strand; and the absence of more than nine GC alignments (table 1). These six siRNAs and nontargeting control siRNA (siControl: 5 0 -UAG CGA CUA AAC ACA UCA A-3 0 ) were synthesized by Sigma GenosyssiRNA service (Sigma-Aldrich Japan K.K., Tokyo, Japan).

Quantitative real-time RT-PCR
Total RNA was extracted using an Isogen-LS RNA extraction kit (Nippon Gene, Tokyo, Japan). Genomic DNA was removed from RNA preparations using the DNase I recombinant RNase-free (Roche, Germany) prior to RT-PCR. After quantification of RNAs using the NanoDrop ND-1000 Spectrophotometer (Nano-Drop Technologies), RNA samples were stored at 2208C, until needed for no more than 2 days. Single-strand cDNAs were synthesized from mRNAs and quantified using the Quanti-TectTM SYBR Green RT-PCR Kit (Qiagen) and ABI Prism 7000 sequence detection system (Applied Biosystems). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) gene was used as an internal control to normalize the cDNA expression levels. GAPDH primer sequences for forward and reverse were 5 0 -TGA TGA CAT CAA GAA GGT GGT GAA G-3 0 and 5 0 -TCC TTG GAG GCC ATG TGG GCC AT-3 0 , respectively. For each target gene, 1 ml of cDNA was amplified in a total volume of 20 ml containing 2 Â QuantiTect SYBR Green RT-PCR master mixes supplemented with 300 nM of each primer. The data were analysed using ABIPRISM 7000 SDS software (Applied Biosystems). For all the samples, crossing points (CP) normalized expression of the target gene versus GAPDH for each treatment or control sample were calculated [26] using the following formula: 2 2(DCP treatment sample -DCP control sample) . Each sample was assayed three times.

RNA-seq analysis
RNA-seq data from a transcriptome analysis of HPV-18þ HeLa cells [27] were downloaded from the European Nucleotide Archive (ENA; accession ERP000959). Sequence reads from this dataset were mapped to the HPV-18 complete genome reference sequence (NCBI Refseq accession NC_001357) using the program BOWTIE2 [28]. HPV-18 gene expression levels were quantified with log 10 -transformed sequence coverage values (i.e. number of mapped tags per position), and mapped reads were visualized along the HPV-18 genome sequence using the Integrative Genomics Viewer [29,30].

Microarray analysis
RNA qualities were examined by capillary electrophoresis using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), and their qualities were confirmed (28S/18S rRNA, E260/E280 ratio . 1.8). Microarray analysis (including labelling, hybridization, image scanning and sample conditions) was performed by Takara Bio Dragon Genomix Center (Mie, Japan). Briefly, fluorescent-labelled cDNAs were generated from 200 ng of total RNA in each reaction using the Agilent Quick Amp labelling kit in two colours. The cyanine 5-labelled cDNAs from the siE1.6-treated HeLa cells were mixed with the same amount of reverse-colour cyanine 3-labelled cDNAs from the non-treated control cells and then applied to a whole human genome oligonucleotide microarray Kit (Agilent Technologies), which contained 41 000 unique human genes and transcripts. These transcripts all contained public domain annotations (http://www.chem.agilent. com/CAG/bsp/gene_lists.asp). Upon hybridization and washing, the Agilent dual-laser DNA microarray scanner scanned the arrays. The data were then extracted from images using the Agilent FEATURE EXTRAC-TION v. 9.5 software, and candidate genes were selected using GENESPRINGGX software (Agilent Technologies). An intensitydependent normalization (LOWESS: locally weighted regression) was applied to correct artefacts caused by nonlinear rates of dye incorporation as well as inconsistencies of the relative fluorescence intensity between the red and green dyes. Differentially expressed genes were then identified by comparing the log 2 ratio of expression levels between conditions with the overall log 2 expression signals (LogRatio versus Log Processed Signal; electronic supplementary material, figure S1).

Gene set enrichment analysis
The potential functional significance of the differentially expressed genes was analysed using gene set enrichment analysis by combining the Ingenuity Pathway Analysis tool (http://www.ingenuity.com/products/ipa) with analyses of custom gene sets taken from the literature (electronic supplementary material, figure S2) [31]. To do this, the enrichment of differentially expressed genes in pathways, or functionally coherent sets of genes, was determined using the Fisher's exact test with the p-value determined as where N is the total number of genes on the microarray, n is the total number of differentially expressed genes, m is the total number of genes in the pathway and k is the total number of differentially expressed genes in the pathway.

Validation of microarray data
Validation of the microarray data was performed using a realtime RT-PCR assay for EP300 and cyclin-dependent kinase inhibitor 1A (CDKN1A) mRNAs using real-time RT-PCR. The primer sequences for EP300 forward and reverse were 5 0 -CTC CAA CTT TCT GCC ACA ACA-3 0 and 5 0 -CCA GAA TCC CTT CCC TTT CG-3 0 , respectively. For CDKN1A, the primer sequences forward and reverse were 5 0 -GGA CCT GGA GAC TCT CA-3 and 5 0 -CCT CTT GGA GAA GAT CAG-3, respectively. The transcript level of each gene was determined using the QuantiTectTM SYBR Green RT-PCR Kit (Qiagen) and ABIPRISM 7000 SDS. GAPDH mRNA was also used as an internal control to normalize the mRNA expression levels.

Expression of a full-length E1 mRNA and protein
Previous reports have shown evidence suggesting that the HPV-18 early genes E6, E7 and E1 are transcribed as a polycistron that terminates within the E1 ORF [10,12]. This rsob.royalsocietypublishing.org Open Biol. 4: 130119 has been taken to indicate the production of a truncated E1 protein with attenuated viral replication activity [14]. However, Western blot analysis of HeLa cell protein extracts revealed a 70 kDa E1 protein [11], which is consistent with the molecular weight calculated from the full-length E1 protein sequence (73.75 kDa) and substantially larger than the calculated molecular weight of the protein corresponding to the shorter E1 transcript (58.55 kDa). To address this apparent contradiction, we re-evaluated HPV-18 gene expression patterns using RNA-seq data from a more recently published transcriptome analysis of HeLa cells [27]. These RNA-seq data support the expression of a full-length E1 mRNA transcript, with the caveat that there may indeed be a higherabundance E6-E7-E1 polycistron that contains a shorter E1 sequence (electronic supplementary material, figure S3). This shorter E1 sequence may be responsible for the previously observed sequestration of cofactors needed to initiate HPV-18 replication in HeLa cells [14]. Nevertheless, in light of our new results on E1 mRNA expression and the previous results on E1 protein expression [11], we conclude that HeLa cells are in fact capable of producing full-length E1 protein product.

Expression data analysis identified 2669 differentially expressed cellular genes after silencing E1 endogenous mRNAs in HeLa cells
The cellular gene expression profile of HeLa cells was compared before and after siE1.6 transfection using the Agilent human oligonucleotide microarray analysis, which contained 41 000 unique human gene transcripts. The microarray data were normalized and the data quality was assessed. Globally, 2669 differentially expressed genes with log 2 ! 1 were identified, comprising 1718 upregulated and 951 downregulated genes. The potential biological significance of these differentially expressed genes was evaluated using gene set enrichment analysis. Upon HPV-18 E1 silencing, we found significant enrichment ( p , 0.05, Fisher's exact test) of differentially expressed genes in four closely related gene sets, including three canonical pathways and one functional group: the toll-like receptor (TLR) signalling pathway, the interferon (IFN) signalling pathway, the antiviral interferon-stimulated genes (ISG) and the apoptosis pathway ( figure 2 and table 2). In addition, the majority of genes in these sets were upregulated after E1 silencing (figures 2 and 3).

Discussion
This study revealed that silencing of the HPV-18 E1 endogenous mRNA expressions in HeLa by siRNA molecules induced changes in the expression of cellular genes, specifically with high expression in three canonical pathways and one functional group: the TLR signalling pathway, the IFN signalling pathway, the antiviral ISG and the apoptotic pathway (figures 2 and 3).
Targeting the TLR signalling pathway in elucidating the cellular and molecular mechanisms of cervical cancer has gained tremendous importance [32]. TLR-4 was overexpressed in cervix cancer, and its activation by LPS promotes proliferation and anti-apoptosis in HeLa cells in vitro [33]. However, members of the interferon regulatory factor (IRF) family of DNA-binding transcription factors have played a role in Table 2. Functional role and statistics on differentially expressed genes of specific pathways and gene sets. growth regulation, antiviral responses and transcriptional induction of IFN-activated early response genes [34] and critical mediators of the induction of early viral transcription and replication in several mucosal HPVs, which require IRF-1 binding to a conserved interferon response element [35]. HPV infections disrupt cytokine expression and signalling with E6 and E7 oncoproteins, and particularly target the type I IFN pathway [36]. However, the precise mechanism in which the IFN complex acts on tumours remains unknown, although its use in clinical practice has been widely recommended, particularly in tumours that are resistant to conventional treatments [37]. The TLRs belong to the T L R s i g n a l l i n g i n t e r f e r o n s i g n a l l i n g a n t i v i r a l I S   . EP300 downregulation after silencing of HPV-18 E1. Confirmation of micro-array data using real-time RT-PCR. Transcript levels of EP300 and CDKN1A were quantified using RT-PCR. EP300 and CDKN1A gene expressions were affected after siE1.6 and siE1.2 transfections, at 100 nM for 48 h, but not after siControl transfection. Decreased EP300 and increased CDKN1A expression were consistent with the results obtained from the microarray analysis. The relative logarithmic expressions were normalized against the non-treatment control (NTC). GAPDH was used as a reference gene control.  figure S3). This shorter HeLa cell-derived HPV-18 E1 protein is truncated at the carboxyl terminus, but still retains carboxyl terminal protein-binding domains that can associate with host cellular factors such as cyclin E/CDK2 [14] and Ini1/hSNF5, a subunit of the SWI/SNF chromatin remodelling complex [23]. Recruitment of the SWI/SNF complex by E1 can alter the local arrangement of transcription factors, replication factors and histones. This rearrangement may reprogram the transcriptional and replication state of DNA. We also found a downregulation of the transcription and translation of EP300 (figure 4), a histone acetyltransferase that epigenetically regulates chromatin remodelling via histone acetylation, which enables access to several transcription factors, thereby activating gene expression [38]. However, the regulatory mechanism of EP300 expression is not well understood. According to a report by Yu et al. [39], early growth response 1 (Egr-1), which binds to GC-rich elements in the promoters, acts upstream of p300/CBP. Egr-1, a zinc finger transcriptional factor, accumulates in the cell nucleus upon stimulation by mitogens and a variety of cytokines, as well as extracellular effects, including hypoxia and 17-b-oestradiol [40,41]. Recently, Saegusa et al. [42] reported that EP300 gene expression is upregulated by Egr-1 binding to the promoter region of the EP300 gene in endometrial carcinoma cells. Although there was no evidence of Egr-1 regulation by HPV-encoded proteins, the mRNA levels of Egr-1 in cervical cancer were significantly higher compared with normal tissue [43]. In this study, we observed a decrease of Egr-1 expression after the silencing of HPV-18 E1 mRNA, suggesting that Egr-1 can be stimulated by E1. However, negative feedback of EP300 to Egr-1 expression was suggested [39][40][41][42]. Moreover, Bouallaga et al. [44] reported that EP300 is recruited by the HPV-18 enhanceosome and is required for HPV18 transcription in HeLa cells. Also, the cell confluence ceiling during the siRNA experiments prior to harvesting was at least 50%, and in control cultures was almost 90%. The above is important because it proves that change in the expression of cell cycle genes such as RB1, CDC25C, CDC2 and CDKN1A has the consequence that cell proliferation is arrested.  An important question is whether the siRNA targeting of E1 also affects the stability (and thus expression) of the E1-E6-E7 polycistronic mRNAs. To address this point, we conducted two additional qPCR analyses. We carefully evaluated both E6 and E7 mRNA expression after siRNA treatments in HeLa cells, and we did not find any significant effect on the expression of E6 and E7 after treatment with siRNAs. Our explanation for these results is based on the fact that the process of maturation of E1-E6-E7 polycistronic mRNAs occurs in the nucleus, and then each one of them (E1, E6 and E7 mRNAs) are exported to the cytosol. On the other hand, many studies have shown that the process of siRNA inhibition is performed mainly in the cytosol. Therefore, we consider there to be little chance that the siRNA targeting of E1 can affect the stability and expression of the E1-E6-E7 polycistronic mRNAs, since they may be in different places in the cell.

Conclusion
Taken together, our data indicate a consistent suppression pattern at the transcriptional level induced by the HPV-18 E1 protein across four intensively collaborating components of the host cell immune machinery (figure 5), suggesting that in addition to its well-understood role in viral replication, the viral E1 protein may also play an important role in mitigating the host's ability to defend against infection.