Host proteome linked to HPV E7-mediated specific gene hypermethylation in cancer pathways

Background Human papillomavirus (HPV) infection causes around 90% of cervical cancer cases, and cervical cancer is a leading cause of female mortality worldwide. HPV-derived oncoprotein E7 participates in cervical carcinogenesis by inducing aberrant host DNA methylation. However, the targeting specificity of E7 methylation of host genes is not fully understood but is important in the down-regulation of crucial proteins of the hallmark cancer pathways. In this study, we aim to link E7-driven aberrations in the host proteome to corresponding gene promoter hypermethylation events in the hope of providing novel therapeutic targets and biomarkers to indicate the progression of cervical cancer. Methods HEK293 cells were transfected with pcDNA3.1-E7 plasmid and empty vector and subjected to mass spectrometry-based proteomic analysis. Down-regulated proteins (where relative abundance was determined significant by paired T-test) relevant to cancer pathways were selected as gene candidates for mRNA transcript abundance measurement by qPCR and expression compared with that in SiHa cells (HPV type 16 positive). Methylation Specific PCR was used to determine promoter hypermethylation in genes downregulated in both SiHa and transfected HEK293 cell lines. The FunRich and STRING databases were used for identification of potential regulatory transcription factors and the proteins interacting with transcription factor gene candidates, respectively. Results Approximately 400 proteins totally were identified in proteomics analysis. The transcripts of six genes involved in the host immune response and cell proliferation (PTMS, C1QBP, BCAP31, CDKN2A, ZMYM6 and HIST1H1D) were down-regulated, corresponding to proteomic results. Methylation assays showed four gene promoters (PTMS, C1QBP, BCAP31 and CDKN2A) were hypermethylated with 61, 55.5, 70 and 78% increased methylation, respectively. Those four genes can be regulated by the GA-binding protein alpha chain, specificity protein 1 and ETS-like protein-1 transcription factors, as identified from FunRich database predictions. Conclusions HPV E7 altered the HEK293 proteome, particularly with respect to proteins involved in cell proliferation and host immunity. Down-regulation of these proteins appears to be partly mediated via host DNA methylation. E7 possibly complexes with the transcription factors of its targeting genes and DNMT1, allowing methylation of specific target gene promoters.


Background
Cervical cancer is the fourth leading cause of cancerrelated death for women worldwide and a particular problem in developing countries. Over 200,000 people, from approximately 550,000 new cases (50.4%), died in 2017 [1]. Human papillomavirus (HPV) is a causative agent of cervical cancer, and almost 90% of overall cervical cancer cases are linked with high-risk HPV infection [2][3][4][5]. About 80-90% of infections are symptomless and removed by the host within a couple of years. However, 10-20% are persistent and can cause carcinoma development [6,7]. There are more than 100 types of HPV, of which 14 types are cancer-causing and defined as high-risk. From this subset, approximately 70% of all cervical cancer cases are linked to HPV type 16 infection [8,9]. E6 and E7 are well known HPV oncogenes that induce cervical carcinogenesis [2]. Despite the implementation of screening and vaccination for HPV, HPV-caused cervical cancer remains an issue among patients who have been infected by high-risk HPV. The elucidation of E6/E7 interactions with host proteins and their effect on host physiology could be helpful in defining HPV elimination strategies [10]. The regulation of host physiology by HPV, particularly with respect to E7, has been extensively studied. It has been found that E7 competes with E2F in binding to retinoblastoma protein (pRb), which promotes cell cycle progression and carcinogenesis by increasing degradation of pRb. E7 is also the most prominent HPV protein contributing to host immune dysregulation and evasion [11]. E7 regulates host physiology not only through interaction with host signaling proteins but also by modulating host epigenetics, including DNA methylation. E7 binds host DNA-methyltransferase 1 (DNMT1) and enhances the methylase activity of DNMT1 [12]. DNMT1 is responsible for cytosine methylation in mammals, which downregulates protein expression, playing an important role in gene silencing. E6 also increases the abundance of DNMT1 and p53 contributing to aberrant host gene promoter methylation; another HPV-associated carcinogenesis mechanism [7,8]. As an oncogenic protein, E7 is likely to methylate a number of genes involved in the hallmark cancer processes, such as cell proliferation and immune evasion. However, the targeting mechanism of E7-directed methylation to certain gene promoters is not fully understood and E7 has no specific binding motif for a host gene promoter. Therefore, we postulated that E7 might bind both to host transcription factors specific to E7 target genes and to DNMT1.
The cell proteome is more closely related to phenotype than readouts from genomics, transcriptomics and epigenetics. Additionally, the relationship between mRNA and protein levels is highly variable between tissue types and cancers and cannot be used for reliable predictions. In this study, we use both proteomics and epigenetics to search for E7-mediated gene methylation, guided by the identification of down-regulated proteins. Our hypothesis; that E7-mediated methylation alters the expression of specific tumor suppressor genes which are important in cervical cancer pathogenesis, is supported by previous studies [13,14]. One study has mentioned that the hypermethylation of p16 and CCNA1 genes stimulate cervical neoplastic progression and contributes to a decrease in cell adhesion molecule 1 (CADM1), which functions in epithelial cell adhesion and is involved in metastasis [14][15][16]. However, this study did not link this activity to E6 or E7. Numerous studies suggest that DNA methylation occurs at the early stages of cervical cancer and in precancerous lesions [17][18][19]. HPV persistence alone is insufficient to predict progression of cervical cancer because additional factors participate in tumorigenesis. Therefore, host DNA methylation analysis combined with HPV testing could be a promising option for predicting progression from precancerous to invasive cancer in HPV-positive women [3,20,21]. This study was designed to find aberrant E7mediated DNA methylation events related to cancer pathways to clarify its influence in cervical cancer progression. We hope this study will provide preliminary data for host DNA methylation states in clinical samples, which may identify useful biomarkers.

Plasmid transformation
The pcDNA3.1-E7 (E7), pcDNA3.1 empty vector (EV) plasmid and DH5α competent cell were thawed on ice for 5 min. The 100 μL of DH5α competent cell was aliquoted into each 1.5 micro-centrifuge tube. Then 5 μL of E7, EV plasmid were added separately into each tube which containing competent cell followed by gently mixed. Then, tubes were incubated on ice for 5 min. Heat shock method was performed, the tubes were placed in 42°C thermomixer well (Eppendorf, USA) for 45 s, immediately the tubes were placed on ice for 2 min. Next, 900 μL of SOC medium was added (Biolabs, USA) into each tube and gently incubated the tubes in thermomixer machine for 45 min at 37°C, 400 rpm. Bacterial cells were collected by spinning down at 8000 rpm for 5 min. The 900 μL of supernatant was discarded and the cell pellet 100 μL was resuspended by pipetting. The 100 μL competent cell contained plasmid was spread on Luria-Bertani (LB) agar (Titan Biotech, India) contained ampicillin antibiotic final concentration 0.1 mg/mL (Merck, Germany). The agar plates were incubated at 37°C in 5% CO 2 incubator.

Plasmid extraction and purification
The positive colonies on Luria-Bertani (LB) agar (Titan Biotech, India) contained ampicillin antibiotic final concentration 0.1 mg/mL (Merck, Germany) were selected and continued cultured in 10 mL Luria-Bertani (LB) broth with 10 μL ampicillin (final concentration 0.1 mg/ mL) at 37°C, 250 rpm shaking incubator overnight. Then, 1000 mL of LB broth with 1 mL of ampicillin (final concentration 0.1 mg/mL) was prepared, 5 mL of overnight cultured bacteria was added then continued cultured at 37°C, 250 rpm shaking incubator overnight. Plasmid extraction (E7, EV) was done afterward using Maxi Plasmid Kit Endotoxin Free according to the manufacturer's instructions (Geneaid, Taiwan).

Cell culture medium preparation
The cell culture medium was performed. The preparation of 1 l of MEM supplemented with 10% FBS required 9.5 g MEM (GE Healthcare, USA), 2.2 g NaHCO 3 (Merck, Germany), 100 mL fetal bovine serum (Gibco, New Zealand) and 900 mL deionized water. The MEM was gently mixed in 1000 mL beaker and steriled filtered through 0.22 μm filter upper cups (Jet Biofil, China) into 1000 mL duran bottle (Duran, Germany). The culture medium was kept at 4°C, and warmed at 37°C before use. The preparation of 1 l of DMEM supplemented with 10% FBS required 13.4 g DMEM (GE Healthcare, USA), 3.7 g NaHCO 3 (Merck, Germany), 100 mL fetal bovine serum (Gibco, New Zealand) and 900 mL deionized water. The DMEM was gently mixed in 1000 mL beaker and steriled filtered through 0.22 μm filter upper cups (Jet Biofil, China) into 1000 mL duran bottle (Duran, Germany). The culture medium was kept at 4°C, and warmed at 37°C before use.

Mammalian cell transfection
The human cervical carcinoma cell line SiHa (HPV type 16-positive), used as E7 expression positive control and C33A (HPV type 16-negative) cell line, used as negative control, were grown and maintained in DMEM supplemented with 10% FBS at 37°C in an atmosphere of 5% CO 2 , HEK293 (Human Embryonic Kidney 293) cell was grown and maintained in MEM supplemented with 10% FBS at 37°C in an atmosphere of 5% CO2. In each 6well plates, SiHa, C33A and HEK293 cells were seeded separately at 2.5× 10 5 cells/mL in 2 mL of growth medium 24 h before transfection. Then 4 μg of pcDNA3.1-E7 (E7), pcDNA3.1 empty vector (EV) plasmid were diluted in 400 μL of Opti-MEM™ serum-free growth medium (Gibco, New Zealand). TurboFect™ transfection reagent (Invitrogen, USA) was briefly vortexed then the reagent 6 μL was added to the diluted DNA. The mixture was mixed immediately by pipetting followed by incubating 20 min at room temperature. The mixture was evenly distributed at the bottom of the well of 6-well plates containing adherent cells. The plates were gently rocked and incubated at 37°C in a CO2 incubator. Protein expression analysis was done after 48 h.

Cell lysis buffer preparation
The cell lysis buffer solution (4% SDS, 0.1 M DTT) were prepared by adding 0.4 g of high-purity SDS (Merck, Germany) and 0.0154 g of DTT (Merck, Germany) to 10 mL of 0.1 M Tris-HCl (Biorad, USA) stock solution at pH 7.6.

Protein extraction
The transfected HEK293 cell, SiHa and C33A cell line were harvested and washed 3 times with sterile 1X phosphate-buffered saline (1XPBS), by diluting 10X PBS (Merck, Germany) ratio 1:9 to deionized water. The cells were centrifuged at 1000×g for 5 min and removed the supernatant. The cell pellet was snap freezing in a microcentrifuge tube using liquid nitrogen (take extreme caution when doing this). The cell pellet was stored at − 80°C until cell lysis was performed. The transfected cells and control cell corresponding to 2 × 10 6 cells (approximately 1 mg pellet of cells) were lysed with 1 mL of lysis buffer. The cell pellet was suspended in the lysis buffer and mixed the solution well using a pipette. Further, the cells were disrupted using a sonicating probe with 20 amplitudes, 10 s pulse on, 10 s pulse off, total 1 min while maintaining the cell lysate at 4°C and centrifuged at 16,000 x g for 20 min at 4°C to collect the supernatant. The lysate was heated for 20 min at 56°C using a heating block to denature the protein. The protein concentration was determined using a BCA assay kit (Thermo Fisher Scientific, USA) as manufacturer's instructions. The extracted proteins have been kept in − 80°C for further experiments.

2D − Clean up
The non-specific reagent and detergent were certainly removed from protein samples before sample preparation by using 2D − Clean up kit (GE Healthcare, USA) according to manufacturer's instructions.

Protein samples preparation for mass spectrometry
The preparation of 100 mM ammonium bicarbonate buffer was done as follow, 0.158 g ABC (Sigma-Aldrich, Germany) dissolved in 20 mL HPLC grade water (Merck, Germany) then 50 mM ammonium bicarbonate buffer was prepared by diluting 5 mL 100 mM ammonium bicarbonate buffer with 5 mL HPLC grade water. The preparation of 100 mM DTT was done as follow, 0.0154 g DTT (Sigma-Aldrich, Germany) in 1 mL 50 mM ammonium bicarbonate buffer. The preparation of 100 mM IAA was done as follow, 0.0185 g IAA (Sigma-Aldrich, Germany) in 1 mL 50 mM ammonium bicarbonate. The preparation of trypsin from porcine pancreas (Sigma-Aldrich, Germany) digestion enzyme was done by adding 200 μL of 50 mM ABC into 20 μg of trypsin to generate the concentration of 0.1 μg/μL. An aliquot of 100 μg of protein was taken and 50 mM ammonium bicarbonate was added to the protein to fill the total volume up to100 μL. Next, 10 μL of 100 mM DTT (final concentration 10 mM DTT) was added in the tube for 40 min at 56°C. Protein Alkylation was done by adding 2 μL of 100 mM IAA (final concentration 20 mM IAA) for 60 min in the dark at room temperature. Protein samples were digested later by adding trypsin (ratio trypsin: sample = 1:50) then mixture was incubated at 37°C overnight. The samples were continued proceeded desalting procedure.

Desalting by oasis column
The preparation of 50% acetonitrile was done by diluting 5 mL of 100% acetonitrile (Merck, Germany) stock solution with 5 mL of HPLC grade water then 0.1% TFA was prepared by diluting 10% TFA stock solution (Merck, Germany) with 49.5 mL of HPLC grade water (Merck, Germany). Finally, elution buffer, 75% ACN / 0.1% TFA was prepared by adding 7.5 mL of 100% acetonitrile stock solution, 0.1 mL of 10% TFA stock solution and 2.4 mL of HPLC grade water. Oasis columns were placed in 15 mL tube, 1 mL of 50% ACN was added and let it flow, 1 mL of 0.1% TFA was added to equilibrate column, let it flow. The equilibration step was repeated 3 times. The protein samples were added into column and slowly let it flow by gravity. Then 1 mL of 0.1% TFA was added and let it flow. Then, oasis columns were taken and placed in the new clean 15 mL tubes, peptide was eluted by adding 1 mL elution buffer, 75% ACN / 0.1% TFA, let it flow by gravity. The samples were taken to evaporate at speed vacumm centrifuge cincentrator machine (Meditop, Thailand), then kept samples in − 80°C until further analysis.

Mass spectrometry analysis
The peptides were identified using an UltiMate® 3000 RSLCnano system (Thermo Fisher Scientific, USA) coupled with Q Exactive Hybrid Quadrupole-Orbitrap mass spectrometer (Thermo Fisher Scientific, USA) through EASY-Spray nano-electrospray ion source (Thermo Fisher Scientific, USA). The samples were loaded with 5-7% Acetonitrile (ACN) in 5 min, 7-45% in 60 min, 45-50% in 5 min, and 50-97% in 5 min, followed by washing at 100% at 300 nL/min flow rate for 90 min. Full MS scan was carried with mass ranges of m/z 200-2000. Precursor ions with + 1 and greater than + 8 charge state were excluded. Fragmentation of precursor ions was performed using Higher-energy collisional dissociation and data acquisition was performed by Thermo Xcalibur 2.2 (Thermo Fisher Scientific, USA).

Protein identification and peptide quantification
The MS Data were analyzed by Proteome discover version 2.2 (Thermo Fisher Scientific). Spectra were matched with protein database using Homo sapiens (Proteome ID: UP000005640) database from Uniprot. Parameters were set to optimized quantification method of label free labeling including; processing workflow were input data of enzyme cleavage as trypsin; Max. Missed cleavage sites were 2; Fragment tolerance and Precursor tolerance were 0.02 Da and 10 ppm, respectively; Modification: Oxidation as dynamic modification and static modification included carbamidomethyl. Consensus workflow was set to result filter for high confidence peptides which PSM was filtered as SEQUEST: XCorr. Chromatographic Alignments included 10 ppm mass tolerance and 10 min was maximum RT shift. The protein results were cut off by false discovery rate lower than 0.5%. Protein abundance was calculated by PD program between E7 transfected and control. In addition, the MS data were proceeding to protein identification and label free quantitation (LFQ) by Maxquant 16.5. In addition, the MS data were proceeding to protein identification and label free quantitation (LFQ) by Maxquant 16.5. For Maxquant database search setting, the criteria were setup as following; peptide tolerance: 20 ppm, MS/MS tolerance: 0.5 Da, Isotope tolerance: 0.1 Da and FDR: 0.05.

Gene ontology
Protein classification was carried out by their biological processes and molecular functions by Proteome discover version 2.2 (Thermo Fisher Scientific, USA) and FunRich software version 3.1.3.

Hierarchical clustering of proteins and bioinformatic analysis
The Hierarchical clustering heatmap for identified proteins by Maxquant was generated by Perseus software version 1.6.1.1 (http://www.perseus-framework.org). The matrix generation contained data from protein identification and quantification which analyzed with Maxquant and Proteome Discoverer 2.2 (Demo version). LFQ intensity were transformed to log2 and missing values were imputed with normal distribution (width:3 and down: 1.8) then transformed values were further clustered and performed density estimation for enrichment expression ratio. Scatterplot was also built by Perseus with use transformed LFQ intensity of EV against E7. Thirteen of down-regulated proteins, were found from Maxquant and PD were selected for further measure of mRNA expression and DNA methylation based on gene ontology related cancer. Furthermore, protein annotation was performed with FunRich software version 3.1.3. Custom databases were generated by upload Homo sapiens accession downloaded from Uniprot proteome (ID: UP000005640) and protein annotation database downloaded through Perseus annotation (annotations.perseusframework.org). Dataset of protein down-regulation in each condition was uploaded to compared gene enrichment analysis. Heatmap was generated based on fold enrichment of biological pathways. Putative regulatory transcription factors of DNA methylation genes were predicted and generated heatmap by FunRich database.

RNA extraction
The plasmid transfected HEK293, SiHa and C33A cell line were performed conventional RNA extraction. Cell culture media was removed from T25 flask and washed with 1XPBS. Then, 1 mL of trizol reagent (Thermo Fisher Scientific, USA) was added to adherent cell and all cells were scraped. Cell lysate was passed up and down through pipette then collected to1.5 micro-centrifuge tube incubated for 10 min then 0.2 mL of chloroform (Thermo Fisher Scientific, USA) was added to 1 mL trizol, the sample was vortexed for 15 s then incubated at room temperature for 3 min. The sample was centrifuged at 12,000 x g 15 min at 4°C. The upper aqueous phase was collected into the next 1.5 micro-centrifuge tube. RNA precipitation was done by adding 250 μL isopropanol (Merck, Germany) then the RNA briefly mixed and incubated at room temperature for 20 min. The sample was centrifuged at 12,000 x g 10 min at 4°C. The supernatant was completely removed and the RNA pellet was washed with 1 mL 75% ethanol (Merck, Germany), mixed by vortex. The mixture was centrifuged at 7500 x g 5 min at 4°C. The washing step was repeated one more time then the supernatant was discarded. The RNA pellet was air dried for 20 min then The RNA pellet was dissolved with 15 μL DEPC-treated water (Thermo Fisher Scientific, USA). The RNA concentration measurement was done by using Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, USA).

RT-qPCR primer design
The target genes were selected and further investigated gene expression analysis in RNA level. Oligonucleotides primers for RT-qPCR were designed manually and sent to the company to be further synthesized by Integrated Device Technology (IDT), Inc., USA (Table 1).

Quantitative reverse transcriptase polymerase chain reaction
The target genes were investigated for gene expression by one-step RT-qPCR technique using KAPA SYBR® FAST One-Step RT-qPCR Master Mix (2X) Kit (Kapabiosystems, USA) according to manufacturer's instructions for Bio-Rad iCycler™ real-time machine (Biorad, USA). The total volume of PCR reaction was 20 μL; KAPA SYBR FAST qPCR Master Mix (2X) 10 μL, 10 mM dUTP 0.4 μL, 10 μM forward primer 0.4 μL, 10 μM reverse primer 0.4 μL, 50X KAPA RT Mix 0.4 μL, Template DNA as required and PCR-grade water (Bioline, UK) up to 20 μL. The qPCR parameters were adjusted properly according to the manufacturer's protocol (reverse transcription 42°C for 5 min, enzyme activation 95°C for 3 min, denaturation 95°C for 3 s and annealing/extension/data acquisition 60°C for 20 s for total 40 cycles, dissociation step is depended on instrument guidelines). The melting curves were carried out in each RT-qPCR to verify singleproduct amplification. The relative level of gene expression was calculated with the Livak method (2 −(ΔΔCt) ). The genes protein Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the reference genes. The measurements were recorded from three technical and three biological replicates for each experimental condition.

DNA extraction
The plasmid transfected HEK293, SiHa and C33A cell line, DNA isolation from cultured cells were performed by using ISOLATE II Genomic DNA Kit (Bioline, UK) followed the manufacturer's instructions. The 10 7 cells were resuspended up in 200 μL lysis buffer GL, 25 μL proteinase K solution and 200 μL of lysis buffer G3 were added then incubated at 70°C for 15 min. DNA binding condition was adjusted by vortex. Then, the 210 μL of 100% ethanol was added to samples and vigorously vortexed. The ISOLATE II Genomic DNA spin columns were placed in a 2 mL collection tube. The samples were loaded to columns and centrifuged for 1 min at 11,000 x g. The flow-through was discarded and collection tubes were reused. The silica membrane was washed by adding 500 μL wash buffer GW1 and centrifuged for 1 min at 11,000 x g. The flow-through was discarded and collection tubes were reused. Then, 600 μL wash buffer GW2 was added and centrifuged for 1 min at 11,000 x g. The flow-through was discarded and collection tubes were reused. The silica membrane was dried by centrifugation for 1 min at 11,000 x g, to remove residual ethanol then ISOLATE II Genomic DNA spin columns were placed in a 1.5 mL collection tube. DNA elution was done by adding 100 μL preheated Elution Buffer G (70°C) onto center of silica membrane and incubated at room temperature 1 min then the DNA collection was done by centrifuged the tubes 1 min at 11,000 x g. The DNA samples have been kept at − 20°C for further purpose.

Bisulfite conversion
The extracted DNA input of 200-500 ng for each sample was subjected to bisulfite treatment using the EpiJET Bisulfite Conversion Kit (Thermo Fisher Scientific, USA) according to the protocol provided by the manufacturer. The 20 μL of DNA sample containing 200-500 ng of purified genomic DNA was added into a PCR tube. The 120 μL of prepared modification reagent solution was added to 20 μL of DNA sample in a PCR tube. The sample was mixed by pipetting up and down, then centrifuged the liquid to the bottom of the tube. The PCR tubes were placed into a T100™ Thermal Cycler (Biorad, USA) and proceeded performing denaturation and bisulfate conversion of DNA: 98°C for 10 min, 60°C for 150 min. Then the next step was proceeded immediately by adding 400 μL of binding buffer to the DNA purification micro columns. The DNA purification micro columns were placed into the collection tube later. The converted DNA sample was loaded into the binding buffer in the columns and mixed completely by pipetting. The micro columns were placed into the collection tubes and centrifuged at 12,000 rpm for 30 s then flow-through was discarded. The micro columns were placed back into the same collection tube. The 200 μL of wash buffer was added into the micro columns and further centrifuged at 12,000 rpm for 30 s, subsequently, the flow-through was discarded. The micro-columns were placed into the same collection tube. The 200 μL of desulfonation buffer, prepared with ethanol, was added into the micro columns and left in the columns at room temperature for 20 min. The micro columns were placed into the collection tubes and then centrifuged at 12,000 rpm for 30 s. The flow-through was discarded. The micro columns were placed into the same collection tubes. The 200 μL of wash buffer, prepared with ethanol was added to the micro columns and centrifuged at 12,000 rpm for 30 s. Then, the micro columns were placed into the same collection tubes. The 200 μL of wash buffer, prepared with ethanol was added to the micro columns and centrifuged at 12,000 rpm for 60 s. Then, the columns were placed into a clean 1.5 mL microcentrifuge tubes, 10 μL of elution buffer was added to the micro columns and centrifuged at 12,000 rpm for 60 s. The converted DNA was eluted and has been ready for downstream analysis. The bisulfite treated DNA were stored at − 20°C.

Design of methylation and unmethylation primer and methylation-specific PCR
The selected target genes were determined which promoter were methylated. The oligonucleotides primers for MSP were designed manually and were sent to the company to be further synthesized by Integrated Device Technology (IDT), Inc., USA. (Table 2 and 3). The target genes were investigated for methylation analysis. The bisulfite treated DNA was consequently used to carry out methylation-specific PCR by using methylated and unmethylated specific primers. MyTaq™ HS Mix (Bioline, UK) was used according to manufacturer's instructions, following standard MyTaq™ HS Mix protocol. The total volume of PCR reaction was 25 μL (MyTaq HS Mix, 2 × 12.5 μL, 10 μM forward primer 0.5 μL, 10 μM reverse primer 0.5 μL, template DNA 200 ng and PCR-grade water (Bioline, UK) up to 25 μL) The qPCR parameters were adjusted properly according to the manufacturer's protocol (initial denaturation 95°C 1 min, denaturation 95°C 30 s, annealing temperature as described for each gene primer set 45 s and extension 72°C 30 s for total 40 cycles, dissociation step is depended on instrument guidelines). EpiTect PCR Control DNA Set (Qiagen, USA) contained methylated human control DNA (bisulfite 100 μL converted) 10 ng/μL, unmethylated human control DNA (bisulfite 100 μL converted)10 ng/μL and unmethylated human control DNA 10 ng/μL were used as internal control. PCR products were observed by 1% agarose (Bio Basic, Canada) by gel electrophoresis using GelDoc™ XR + (Biorad, USA). The methylated and unmethylated band intensities of each sample were visualized and measured using GelQuant.NET Software (Biochem Lab Solutions, CA, USA).

Statistical data analysis
All data with average of triplicates was reported. Paired sample t-tests was applied for two groups comparison in each experiment, p ≤ 0.05 was indicated statistical significance.

E7 protein expression in a transfected HEK293 cell line
To examine the transfection efficiency of a recombinant histidine-tagged E7 construct in HEK293 cells, the protein expression in lysates from E7-transfected, empty vector (EV) transfected and untransfected HEK293 cells was examined by immunoblot using anti-E7 antibody. As shown in results from immunoblotting in Fig. 1, E7 bands were observed at 20 kDa in E7-transfected HEK293 samples and at 17 kDa in SiHa. E7 protein bands were not observed in EV-transfected and untransfected HEK293 and C33A lysates. β-actin was used as internal control, which presented a 42 kDa band.

E7-induced alteration of HEK293 proteome
We further explored E7-driven host proteome changes, comparing E7-and EV-transfected HEK293 cell proteomes. The Maxquant database search algorithm providing a greater number (Additional file 1: Figure S1 and Additional file 2: Table S1 and Additional file 3: Table S2) compared with SEQUEST (passed FDR ≤ 0.05). Several proteins were found at lower levels in E7-transfected HEK293 than in EV-transfected HEK293 (Fig. 2). Proteome-wide expression analysis showed differential expression after hierarchical heatmap clustering (Fig. 2b). We focused on downregulated proteins (potentially resulting from gene silencing associated with aberrant host DNA methylation). According to hierarchical clustering, protein expression pattern was clustered into 10 groups. Clusters 1 and 7 contained down-regulated proteins relative to other clusters. We also analyzed down-regulated proteins in a scatter plot (Fig. 3a). Proteins with lower than two-fold down-regulated expression (as defined by Maxquant transformed label-freequantification) were selected for further function annotation. From this information, we constructed a pathway enrichment heatmap for nominating candidate genes relevant to cancer pathways.

Protein functional enrichment analysis
Identified proteins from E7-and EV-transfected cells were analyzed using FunRich protein annotation analysis, which combines gene enrichment with protein quantity data (Fig. 3b). Major differences in functional enrichment between E7-and EV-transfected cells are observed for several protein classes. In E7-transfected HEK293, proteins were increased related to; ubiquitin-specific protease activity (4.3%), transcription factor activity (2.4%), transcription factor binding (0.2%) and transcription regulator activity (4.6%). Conversely, protein functions for; the structural constituents of ribosomes (9.7%), catalytic activity (6.5%) and heat shock protein activity (1.3%) were decreased. Another 8.9% of major differences were attributed to proteins of unknown function.

Down-regulated proteins preferentially expressed in particular gene ontology (GO) annotation pathways
To identify host gene candidates for specific E7-mediated DNA promoter hypermethylation, we mapped downregulated proteins against several cancer-associated pathways. The protein annotation heatmap (Fig. 4a) indicated the function of E7-mediated down-regulated proteins, which mostly mapped to the interferon gamma pathway, the mTOR signaling pathway, the ErbB receptor signaling network and VEGF/VEGFR signaling network. In addition, there were a smaller number of proteins mapped to other major cancer-associated pathways such as immune systems and p53 pathway. We selected indicative candidate proteins for mRNA transcript level and promotor DNA methylation status analysis. Candidates were selected that were distributed in cancer associated pathways related to host immunity, the cell cycle and HPV, that have been described in the literature.     (Fig. 4b).

Prediction of putative regulatory transcription factors of E7-mediated host specific promoter hypermethylation
We conducted analysis to determine potential host transcription factors interacting with E7 to implement specific gene-promoter hypermethylation. The heatmap in Fig. 7a identified GA-binding protein alpha chain (GABPA), specificity protein 1 (SP1) and ETS Like protein-1 (ELK1) as possible transcription factors for C1QBP, BCAP31, CDKN2A and PTMS. We used the STRING database to identify potential interacting protein partners with E7 and these host transcription factors. Protein networking prediction by STRING database demonstrated GABPA interacts to SP1, which interact to E7 interacting proteins i.e. p53 and Rb1 (Fig. 7b). Thereby, apparently, GABPA, SP1 and ELK1 potentially interact a group of proteins that function in HPV associated pathways as well as transcription factors binding, in which CDKN2A also the key protein in these pathways as well. Moreover, C1QBP, BCAP31 and PTMS protein were predicted as their interactor pathways in cancer.

Discussion
The expression of E7 altered the proteome of HEK293 cells Although it is commonly known that HPV proteins E6 and E7 play crucial and interlinked roles in cancerous tissue transformation, the effects of E7 alone in tumorigenesis are not fully understood. Particularly with respect to the links between E7-driven DNA methylation of specific gene promotors at early cancer stages and host proteome aberrations. We employed HEK293 as naïve epithelial cell model, previously unexposed to HPV. Despite HEK293 not being a cervical cell type, a previous study reported that this line can be used for studying this cancer [22]. In addition, HEK293 cells generally have high transfection efficiency, and thus we could be confident that E7 transfection would reliably influence the proteome. Our proteomic analysis revealed that E7 likely influenced protein levels through transcription factor/regulator binding and modulation, while reducing ubiquitin-specific protease activity and heat shock protein activity. Gene enrichment analysis implied that E7 also modifies host transcription from an epigenetic level. DNA methylation is a commonly employed viral mechanism to suppress specific protein expression in the host proteome. Here, we chose down-regulated proteins from proteomic analysis, that play crucial roles in cancerassociated pathways, for our analysis of the corresponding hypermethylated gene candidates. Comparing the selected down-regulated proteins with the expression in other types of cancers in human atlas database, these down-regulated proteins mostly express less compared to normal tissues. Therefore, this further convinced us that E7 influenced the proteins have been identified that down-regulation in cancer tissues. In which, the transcript abundances of these proteins were measured to strengthen the link between E7mediated protein downregulation and corresponding gene hypermethylation. Expression of candidate mRNA transcripts compared with corresponding protein levels and E7-mediated DNA methylation of gene candidates We measured the relative abundance of mRNA transcripts for the 13 gene candidates by qPCR, to determine the correlation to the corresponding protein levels assessed in proteomics data. Data were compared with expression in the SiHa cell line.
Six genes out of the 13 candidates were down-regulated at the transcriptional level (C1QBP, CDKN2A, ZMYM6, BCAP31, HIST1H1D and PTMS), while four were upregulated (EIF3K, RACK1, TAPBPL and BTNL8), however, proteomic data demonstrated all their corresponding protein abundances were down-regulated, In SiHa cells, the 13 gene candidates were expressed in a corresponding pattern to the transfected HEK293 cells, suggesting that E7 drives the expression pattern of these genes candidates in same manner that HPV drives their expression in cancer. However, SiHa also expresses E6 making it inconclusive that E7 alone is able to drive tumorigenesis through dysregulation of those 13 genes. Nevertheless, it is evident that the expression of those 13 genes is influenced by E7 from proteomic and qPCR results in E7-transfected HEK293 cells. In addition, four down-regulated genes in our transformed HEK293 cells, with corresponding expression in SiHa and linked to HPV-altered signaling pathways (C1QBP, CDKN2A, BCAP31 and PTMS), were selected for methylation analysis of their promotors. Interestingly, we found all their promoters were hypermethylated in both E7-transfected HEK293 and SiHa, whereas methylation of those genes was lower in EV-transfected HEK293 and C33A. This result implies that E7 modulates Fig. 7 Prediction of putative regulatory transcription factors and STRING database protein network analysis. Heatmap of regulatory transcription factors of hypermethylated gene candidate promoters, the heatmap demonstrated that GABPA, SP1 and ELK1 were putative transcription factors that possibly participated in DNA methylation (a). STRING database protein network analysis mapping were performed for seeking interacting proteins that presumably associate with transcription factors contributing E7-mediated hypermethylation. The STRING network analysis indicated that E7-mediated hypermethylated genes and putative transcription factors were mostly relevant to cell cycle, viral carcinogen and HPV infection (b) C1QBP, CDKN2A, BACAP31 and PTMS levels through DNA methylation and that this modulation could play a role in precancerous progression, though it should be noted that methylation of those genes in the cancerous C33A cell line was low.
Here, we postulated E7 binds to host transcription factors to form a complex with DNMT1 and binds to a specific genome region via transcription factors. Subsequently, DNMT1, which has enhanced methylase activity when bound to E7, methylates the specific genome region. We constructed a predictive heatmap to ascertain putative regulatory transcription factors for C1QBP, CDKN2A, BCAP31 and PTMS. From our analysis, GABPA demonstrated the highest possibility for regulating the four gene candidates followed by ELK1 and SP1. GABPA, SP1 and ELK1 have common cis-regulatory elements and probably form a complex as a result of E7 activation [23]. A previous study has mentioned that GABPA likely interacts with SP1 and plays a part in regulating broad biological processes such as; embryonic development, cell differentiation, mitochondrial biogenesis and the cell cycle [24]. GABPA and SP1 have both been linked to HSV type 16-mediated transcription regulation, and moreover, are able to facilitate DNA methylation of the E7 upregulated genes identified here [25,26]. GABPA belongs to the EST transcription factor family, which are able to bind to several cis-regulatory elements, including C1QBP, CDKN2A, BCAP31 and PTMS. Suppression of these genes by hypermethylation, and a reduction in their protein abundance, possibly promotes carcinogenesis. Notably cancerous cells often express CDKN2A at low levels mostly, in which through gene methylation [27,28].

Biological function of hypermethylated gene promotors
The effects of HPV-E7 oncoprotein on host gene promoter methylation have been extensively studied. For example, E7 affects the methylation of E-cadherin, involved in Langerhans cell delivery through the mucosal epithelium layer, to instigate an immune response to HPV infection [29][30][31][32][33]. E7 possibly also drives hypermethylation of the cyclin A1 (CCNA1) promoter, resulting in decreased cyclin A1 mRNA expression and uses the same mechanism to repress HLA-E expression, downregulating MHC class I, the antigen presenting complex on natural killer and CD8 + T cells [14][15][16].
C1QBP encodes the multi-functional protein complement component 1Q subcomponent binding protein (C1QBP), extensively found in various tissues and cell types. C1QBP is mainly found in mitochondria but also, occasionally, in the cytosol, cell surface and nucleus of many cell types [34][35][36]. Several studies have revealed the functions of C1QBP in tumorigenesis are varied: for instance, C1QBP maintains oxidative phosphorylation for tumor cells, enhances cell chemotaxis and metastasis in breast cancer and suppresses the Y box binding protein 1 in renal cell carcinoma [37][38][39][40][41][42]. Additionally, suppression of C1QBP by HPV type 16, dramatically induced cancer cell immune evasion [42,43]. C1QBP promoter methylation could impede the host immune response by modifying interleukin 6 production and dendritic cell metabolism and maturation; activity important in cervical cancer cell progression [44,45]. However, to our knowledge, no report has mentioned that E7 mediates suppression of C1QBP. Our findings here demonstrate that C1QBP promoter methylation is significantly increased in E7-transfected HEK293 compared with control and is a subject of E7-targeted methylation.
BCAP31, coding for B-cell receptor-associated protein 31 (BCAP31), is an endoplasmic reticulum membrane protein that acts as a molecular chaperone. It is associated with apoptosis regulation through the caspase 8 pathway, protein transport and degradation, promoting the vesicular transport of transmembrane proteins (such as MHC 1, CD11b/CD18, cytochrome 450 and cellubrevin) and B-and T-cell activation [46][47][48][49][50][51][52][53]. Enhanced BCAP31 expression has been observed to decrease survival of non-small-cell lung carcinoma cells [54]. However, overexpression in gastric cancer stimulated proteasome degradation of p27kip1, leading to cell proliferation [55]. Interestingly, from a previous study, BCAP31 has been established as a key target of high-risk HPV E5 protein, modulating cancer cell differentiation in cervical cancer [56]. To our knowledge, the effect of E7 on BCAP31 had not been identified. Here, we found that BCAP31 protein and BCAP31 gene expression levels were decreased in E7-transfected HEK293 compared with control and that promoter methylation status in E7-transfected cells was markedly increased. This result suggests that E7, in addition to E5, contributes to HPV pathogenesis by targeting BCAP31 through gene silencing.
CDKN2A, encodes the p16 INK4a tumor suppressor protein of the cyclin dependent kinase inhibitor family, which functions by blocking cyclin dependent kinase 4/6 interaction with cyclin D1 in the cell cycle. This event can heavily suppress phosphorylation of pRb, prohibiting cell cycle progression from G1 to S phase [57]. Inactivation of p16 INK4a can strongly support cell proliferation in various types of carcinomas, including gastric cancer, glioma, bladder cancer, breast cancer and head and neck cancer [58,59]. Dysregulation of the CDKN2A gene has also been a prominent observation in oral squamous cell carcinoma [60][61][62][63]. Suppression of CDKN2A gene expression generally occurs by aberrant methylation of its promoter, and it also plays a crucial role during neoplastic progression in cervical cancer [64,65]. The methylation of CDKN2A promoter may occur between low-grade and high-grade cervical dysplasia and is common in invasive carcinoma which express the E7 oncoprotein [66,67]. Here, we observed that promoter methylation of CDKN2A was significantly increased in E7-transfected HEK293 compared with control. This result is consistent with protein and gene expression level in a recent study identifying CDKN2A promoter methylation as a diagnostic biomarker of E7-expressing cervical dysplasia [68].
PTMS encodes for parathymosin alpha (PTMS), an immunoregulatory protein which elicits interleukin 2 expression in human lymphocytes, and cooperates with αB-interferon in stimulating natural killer cell activity and involved in the early DNA replication and apoptosis process [69][70][71][72]. However, PTMS has not been studied in cervical cancer intensively, particularly HPV-mediated gene methylation. In addition, our results show high methylation status for the PTMS promoter in E7-transfected HEK293, compared with control, in agreement with low expression levels of protein and mRNA. These results suggest that E7-targeted methylation is responsible for PTMS activities in immune regulation and apoptosis activation enhancement by HPV. However, in summary, this study still requires further research on to validate the interaction of gene candidates, transcription factors and how E7 interacts with specific transcription factors by other experiments such as chromatin immunoprecipitation and further gene methylation in clinical samples.

Conclusion
This study has identified proteins downregulated by the HPV E7 oncoprotein and linked them to targeted promoter hypermethylation of their corresponding genes. In particular, four genes, C1QBP, BCAP31, CDKN2A and PTMS are identified. These genes may serve as potential biomarkers for HPV infection or be target genes for cervical cancer therapeutics. Our findings are a primary report, and further investigations into the roles of these genes in HPV-mediated cervical cancer progression will be important to improve cervical cancer detection and treatment.
Additional file 1 : Figure S1. Venn diagram of comparison of number of identified proteins with SEQUEST (Proteome Discoverer 2.2) and Maxquant database searching algorithm. The Venn diagram showed that Maxquant identified more proteins compared to SEQUEST but are mostly common. The commonly identified proteins by two algorithms were selected to further transcript abundance determination.
Additional file 2 : Table S1. List of total identified proteins of E7transfected HEK293 by Proteome Discover.
Additional file 3 : Table S2. List of total identified proteins of E7transfected HEK293 by Maxquant.