A transcriptomic analysis of malignant transformation of human embryonic esophageal epithelial cells by HPV18 E6E7
Original Article

A transcriptomic analysis of malignant transformation of human embryonic esophageal epithelial cells by HPV18 E6E7

Duo Tang1, Biqi Wang1, Sara Khodahemmati2, Jingtao Li1, Zhixiang Zhou1, Jingfeng Gao2, Wang Sheng1, Yi Zeng1

1College of Life Science and Bioengineering, Beijing University of Technology, Beijing 100124, China; 2College of Environmental and Energy Engineering, Beijing University of Technology, Beijing 100124, China

Contributions: (I) Conception and design: Z Zhou, D Tang; (II) Administrative support: Z Zhou, J Gao, WSheng, Y Zeng; (III) Provision of study materials or patients: J Li, Y Zeng; (IV) Collection and assembly of data: Z Zhou, D Tang, S Khodahemmati, B Wang; (V) Data analysis and interpretation: D Tang, B Wang; (VI) Manuscript writing: All authors; (VII) Final approval of manuscript: All authors.

Correspondence to: Dr. Zhixiang Zhou. College of Life Science and Bioengineering, Beijing University of Technology, Beijing 100124, China. Email: zhouzhixiang@bjut.edu.cn.

Background: Esophageal cancer is one of the most common malignant tumours in humans. A series of esophageal cancer cell lines are accompanied by human papilloma virus (HPV) infection, but the mechanism behind HPV in cancer malignancy is not clear.

Methods: This research was conducted in different generations of HPV E6E7 gene-induced human foetal esophageal epithelial immortalised cells (Shantou Human Embryonic Esophageal Epithelial cell line; SHEE); the RNA sequencing transcriptomic analysis was performed to explore the mechanism of HPV infection in these cell lines.

Results: The results showed that there are 9,990 differential genes in late-stage cells compared with HPV18 E6E7-infected early foetal esophageal epithelial immortalised cells. Among these, 4,882 genes are upregulated, and 5,108 genes are downregulated. We used bioinformatics to analyze the expression and function of aberrantly expressed lncRNA, miRNA, mRNA and construct the competing endogenous RNA (ceRNA) network and protein protein interaction (PPI) network.

Conclusions: we predicted TP53TG1 promotes to malignant transformation of SHEEs by acting as a ceRNA to competitively bind to miR-6835 and regulate IGF2 expression. We also predicted IL6 serve as prognostic biomarkers and therapy target. With these results maybe provides new insights into the mechanisms of HPV carcinogenesis in esophageal cancer.

Keywords: Human papilloma virus (HPV); esophageal cancer; RNA-Seq; qRT-PCR


Submitted Sep 19, 2019. Accepted for publication Jan 18, 2020.

doi: 10.21037/tcr.2020.02.23


Introduction

Esophageal cancer is one of the most common malignant tumours in humans. Data from 2012 show that there were 455,800 of new cases of esophageal cancer worldwide, and of these cases, 400,200 died (1). One problem is that esophageal cancer lacks the typical early symptoms or is misdiagnosed as other diseases, such as cardiac spasms. It demonstrated that most advanced esophageal cancer patients were unable to undergo surgery, meaning their survival rate is less than 10% in 5 years after diagnosis. However, once esophageal cancer is diagnosed early and treated promptly, the 10-year survival rate of patients with esophageal cancer can reach more than 95% (2). Therefore, it is particularly important to develop a new diagnostic strategy to improve the cure rate of esophageal cancer.

HPV is a double-stranded small DNA virus and has been reported to be closely related to various tumours (3). In recent decades, the works of our laboratory and other international laboratories have shown that HPV is also an important cause of esophageal cancer (4-12). More than 80% of esophageal cancer patients in China’s Shantou area are concomitant with HPV 16 and 18 infections. E6 and E7, two HPV early oncoproteins, are responsible for malignant transformation. The E6 protein can bind and degrade p53 and destroy the normal control of the cell cycle, while the E7 protein can bind with the pRb protein, hence, the pRb losses its ability to be an antioncogene. The effects of the E6 and E7 proteins lead to malignant cells to proliferate (13). Therefore, we established an immortalised esophageal epithelial cell line SHEEs with an adeno-associated virus containing HPV18 E6E7 to infect human foetal esophageal epithelial cells, these transformed cells are more likely to progress to anchorage-independent growth and tumorigenicity and may provide an in vitro model of carcinogenesis induced by HPV (14,15). However, it is still not clear if genetic and cellular mechanisms take place in the immortalisation and malignant transformation of human foetal esophageal epithelial cells induced by HPV18 E6E7.

In the current work, we passaged the SHEE cell to 26th generation (P26) and 79th generations (P79), and then, the RNA-Seq analysis was used to compare the transcriptional differences between the two generation cells. Previous studies have shown, the SHEE cell from the 10th to the 85th passage were changed gradually from preimmortal, immortal to malignantly transformed stages (14). In our work, we passaged the SHEE cell to 26th generation (P26) which is an early immortalized cell and 79th generations (P79) which is a pre-malignant transformation cell. And then, the RNA-Seq analysis was used to compare the transcriptional differences between the two generation cells.


Methods

Cell culture

SHEE cells induced by the HPV 18 E6 and E7 genes were obtained from the Institute of Virology, Chinese Academy of Preventive Medicine. The SHEE cells were cultured in Dulbecco’s modified Eagle’s media (DMEM; Gibico, USA) supplemented with 10% fetal bovine serum (FBS; Gibico, USA), 100 U/mL penicillin and 100 mg/mL streptomycin (Hyclone, USA) and then incubated at 37 °C, 5% CO2. SHEE cells were taken from 26 passages (P26) and 79 passages (P79) for further analysis.

Quantification and identification of RNA

RNA degradation and contamination were detected on a 1% agarose gel. RNA purity was detected using a NanoPhotometer® spectrophotometer (IMPLEN, CA, USA). The RNA concentration was measured using Qubit® RNA assay kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA), and the RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Transcription groups

The total amount of each sample was 3 µg of RNA used as an experimental material for RNA sample preparation. A sequencing library was generated by using the NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA), per the instructions, and an index code was added to the attribute sequence for each sample. The mRNA was purified from the total RNA by using magnetic beads attached with poly-T oligo. In the NEBNext first-strand synthesis reaction buffer (5X), divalent cations were used for fragmentation as the temperature increased. The first-strand cDNA was synthesised using a random hexamer primer and M-MuLV reverse transcriptase (RNase H). Second-strand cDNA synthesis was subsequently performed using DNA polymerase I and RNase H. The remaining overhangs were converted to blunt ends by exonuclease/polymerase activity. After being adenylated at the 3' end of the DNA fragment, a NEBNext adaptor with a hairpin loop structure was ligated to prepare for hybridisation. To select the cDNA fragments of, preferably, 150–200 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, USA). The cDNA that was ligated with a size-selected adaptor was then used at 37 °C for 15 minutes using 3 µL USER enzyme (NEB, USA) and then at 95 °C for 5 minutes. A PCR was carried out using the Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. Finally, the PCR product (AMPure XP system) was purified, and the library quality was assessed using an Agilent Bioanalyzer 2100 system.

Clustering and sequencing

Clustering of the samples was performed on the cBot Cluster Generation System using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina). After clustering generation, library preparations were sequenced on an Illumina Hiseq platform, generating paired-end reads of 125 bp/150 bp.

Quantification of gene expression levels

HTSeq V0.6.0 was used to count each gene reading. The fragments per kilobase of transcript per million (FPKM) of each gene was then calculated based on the length of the gene, and the reading count was mapped to the gene. The FPKM takes into account the effect of sequencing depth and gene length using how many reads are used and is the most commonly used method for estimating gene expression levels.

Differential expression analysis

Differential expression analysis of two conditions/groups (two biological replicates per condition) was performed using the DESeq2 R package (1.10.1). DESeq2 provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P value <0.05 found by DESeq2 were assigned as differentially expressed.

Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis of two conditions was performed using the edgeR R package (3.12.1). The P values were adjusted using the Benjamini and Hochberg method. Corrected P value of 0.05 and absolute fold change of 2 were set as the threshold for significantly differential expression.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of differentially expressed genes

A GO enrichment analysis of differentially expressed genes (DEGs) was performed using the cluster Profiler R package, in which the gene length deviation was corrected. GO items with a P value of less than 0.05 correction were considered to be significantly enriched by DEGs. The KEGG is a database resource for understanding the advanced functions and effects of biological systems, such as cells, organisms and ecosystems, hence providing molecular-level information, especially genome sequencing and other high-throughput, large-scale molecular data (http://www.genome.jp/kegg/). We used clusterProfiler R package to test the statistical enrichment of differential expression genes in KEGG pathways. KEGG enrichment is significantly enriched with a padj of less than 0.05.

Construction of the ceRNA network

The putative lncRNA-miRNA-mRNA interactions were evaluated using the algorithms of Targetscan version 7.2 (http://www.targetscan.org/) and ENCORI (http://starbase.sysu.edu.cn/). The ceRNA network was constructed and visually displayed using the Cystoscope software V3.7.1 (https://cytoscape.org/).

Construction of the PPI network

To study the PPI of malignant transformation of human embryonic esophageal epithelial cells by HPV18, Omicsbean (http://www.omicsbean.cn/) online tools were used to construct a PPI network for significantly differential mRNAs.

Survival analysis

To analyse the prognostic significance of DEGs, survival analysis was conducted using clinical information from TCGA. Online tools GEPIA (http://gepia.cancer-pku.cn/) and OncoLnc (http://www.oncolnc.org/) were used to survival analysis of DEGs.

Real-time qPCR

We selected seven DEGs involved in apoptosis, migration and inflammation for quantitative reverse transcriptase PCR (qRT-PCR) validation of the RNA-Seq results. The PCR primers were designed using Oligo7 software and are listed in Table 1. Total RNA was extracted from cells using trizol reagent (Invitrogen,USA). The 1 µg of total extracted RNA was reverse transcribed in a final volume of 20 µL to synthesise the first-strand cDNA using a PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa Biotechnology, Dalian, China) per the instructions of the manufacturer. Gene expression was detected by qRT-PCR using SYBR® Premix DimerEraser TM (TaKaRa Biotechnology, Dalian, China). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an endogenous control to normalise expression, and SHEE P26 was used as the experimental control. In total, 10 µL of PCR reactions in quadruplicate were carried out by an initial denaturation at 95 °C for 30 s followed by 40 cycles, each consisting of 30 s at 95 °C, 30 s at 55 °C, 30 s at 72 °C and then one cycle for melting curve that consisted of 0 s at 95 °C, 15 s at 65 °C and 0 s at 95 °C. The 2−ΔΔCt method was used to calculate the relative gene expression levels.

Table 1

Primer sequence for real-time qPCR

Symbol Oligonucleotide (5' to 3')
MAGEA3 Sense: TCACTCTTTGAAGCGAGCAG; antisense: TGCTGACGCTCATTCAACCA
IGFBP5 Sense: GCTCAACTTCGAAAATGGCAAC; antisense: TCTGTTCATCCAACTTGCCT
SPP1 Sense: TATGCTGGTTGTAGACCCCAA; antisense: ATTGACCTCAGAAGATGCACT
Linc00958 Sense: CACGTTTTATTTCTGACTGCT; antisense: AGTGGACTCATCTTTGCCT
Spock1 Sense: CCTTTGGCTTCGTAAATAACCC; antisense: CCCAAAACAGAGTTTGCCTA
IL6 Sense: ACTCACCTCTTCAGAACGAATTG; antisense: CCATCTTTGGAAGGTTCAGGTTG
S100A7A Sense: TCCATAATAGGCATGATCGAC; antisense: TCAATCTTGCCATCACGTC
GAPDH Sense: TGTTGCCATCAATGACCCCTTC; antisense: AGCATCGCCCCACTTGATTTTG

Statistical analysis

For the obtained experimental data, the mean (± standard deviation) of at least three independent experiments was calculated. Using an analysis of variance (ANOVA), the statistical differences were considered to be significant when the P value was less than 0.05.


Results

Analysis of DEGs between P26 and P79 cells

Using DESeq2 software to perform a differential expression analysis, the padj was set to less than 0.05 as the difference significance criterion, and a total of 9,990 differential genes were obtained, including 4,882 upregulated genes and 5,108 downregulated genes (Figure 1). We focused on genes with an absolute value >5 of expression fold change (log2Fold Change, log2FC), that is, genes with expression MAGEA3, CAPN6, S100A7A, SPP1, SLCO2B1, IGFBP5, MT1E, S100A7, RIBC2, IQGAP2 were significantly up-regulated. It is worth noting that MAGEA3 is most significant up-regulated gene. Six DEGs were down-regulated, they are LINC00958, IL6, KRT6A, ELFN2, C6, SPOCK1. LINC00958 is most significant down-regulated gene in all differential genes (Table 2, Figure 2).

Figure 1 Differential gene volcano map. The volcano map can visualize the overall distribution of genes with the significant difference, and the horizontal axis represents the fold change in the expression of the gene in different samples (Log2Fold Change), and the vertical axis indicates the significant difference in expression (−Log10Padj). The upward gene represented by the red dot and the downward gene indicated by a green dot.

Table 2

Differential genetic classification of prominent type between SHEE P26 and P79, as measured by RNA-Seq. Log2FC >5.0, Padj ≤0.05

Gene name Gene type Log2FC Gene description
Up-regulated
   MAGEA3 Protein coding 8.027788 MAGE family member A3
   CAPN6 Protein coding 6.882745 Calpain 6
   S100A7A Protein coding 6.309254 S100 calcium binding protein A7A
   SPP1 Protein coding 6.285962 Secreted phosphoprotein 1
   SLCO2B1 Protein coding 6.177653 Solute carrier organic anion transporter family member 2B1
   IGFBP5 Protein coding 6.128139 Insulin like growth factor binding protein 5
   MT1E Protein coding 5.542693 Metallothionein 1E
   S100A7 Protein coding 5.487781 S100 calcium binding protein A7
   RIBC2 Protein coding 5.077229 RIB43A domain with coiled-coils 2
   IQGAP2 Protein coding 5.076827 IQ motif containing GTPase activating protein 2
Down-regulated
   LINC00958 lncRNA −6.44381 Long intergenic non-protein coding RNA 958
   IL6 Protein coding −5.9908 Interleukin 6
   KRT6A Protein coding −5.32287 Keratin 6A type II
   ELFN2 Protein coding −5.24081 Extracellular leucine-rich repeat and fibronectin type III domain containing 2
   C6 Protein coding −5.17227 Complement component 6
   SPOCK1 Protein coding −5.07011 Sparc/osteonectin cwcv and kazal-like domains proteoglycan (testican) 1
Figure 2 Prominent differential gene heat map. Red indicates high gene expression, while Blue indicates low gene expression.

Differentially expressed lncRNAs between P26 and P79 cells

There are 82 lncRNAs were identified in differential expression analysis (Figure 3). Thirteen lncRNAs were found to be aberrantly expression, an absolute value >5 of log2FC and the padj was set to less than 0.05. Six lncRNAs were downregulated and 7 lncRNAs were upregulated (Table 3, Figure 4).

Figure 3 A heatmap of all lncRNAs in DEGs. Red indicates high gene expression, while Blue indicates low gene expression.

Table 3

Differential lncRNAs between SHEE P26 and P79, as measured by RNA-Seq. Log2FC >2.0, Padj ≤0.05

Gene name GeneType Log2FC Padj
LINC00958 LncRNA −6.44381 2.90E-60
LINC00941 LncRNA −3.81448 6.93E-15
LINC01127 LncRNA −2.9287 6.05E-08
LINC01293 LncRNA −2.83693 3.91E-10
LINC01085 LncRNA −2.27758 5.04E-05
LINC00452 LncRNA −2.0481 3.01E-04
CASC15 LncRNA 2.13179 3.86E-26
C7orf65 LncRNA 2.282585 5.18E-13
LINC01460 LncRNA 2.285909 3.32E-06
LINC00877 LncRNA 2.451322 1.46E-30
LINC00885 LncRNA 2.647319 1.01E-65
LINC01234 LncRNA 3.437637 3.48E-60
ZNF503-AS1 LncRNA 4.048393 1.79E-77
Figure 4 A heatmap of differentially expressed lncRNAs with an absolute value >5 of log2FC and the padj <0.05.

Enrichment analysis

GO category and KEGG pathway were analyzed for relationships among all DEGs in human embryonic esophageal epithelial cells were analysed using GO and KEGG. A GO analysis is a comprehensive database describing gene function and can be divided into three parts: biological process (BP), molecular function (MF), and cellular component (CC). GO analysis results showed that changes in BP were enriched in ribosome biogenesis, ribonucleoprotein complex biogenesis, DNA replication, in utero embryonic development, ncRNA processing, response to endoplasmic reticulum stress, cell cycle phase transition, mitotic cell cycle phase transition, response to type I interferon, DNA-dependent DNA replication. Changes in MF mainly enriched in transcription factor binding, damaged DNA binding, structure-specific DNA binding, double-stranded DNA binding, ubiquitin protein ligase binding, ubiquitin-like protein ligase binding, protein serine/threonine kinase activity, magnesium ion binding, ubiquitin-like protein transferase activity, helicase activity. Changes in CC mainly enriched in nucleolar part, mitochondrial inner membrane, organelle inner membrane, nuclear membrane, nuclear envelope, preribosome, mitochondrial matrix, centrosome, focal adhesion, transferase complex and transferring phosphorus-containing groups (Figure 5). KEGG is a comprehensive database that integrates genomic, chemical and system function information. As shown in Figure 6, the three most affected pathways are ribosome biogenesis in eukaryotes, protein processing in the endoplasmic reticulum and the cell cycle. Table 4 shows the expressions of the genes involved in each pathway.

Figure 5 Gene ontology (GO) enrichment analysis of DEGs between SHEE P26 and P29.
Figure 6 The first 10 biological functions identified by KEGG analysis. * for significant differences.

Table 4

Gene expression in significant anomalous pathway in KEGG analysis

KEGG pathway Count P value Genes
Ribosome biogenesis in eukaryotes↓ 61 6.01E-08 REXO2↑, XRN2↑, CSNK2A2↑, NOL6↓, NOP56↓, DKC1↓, WDR43↓, GNL3↓, TCOF1↓, RCL1↓, RPP25↓, GTPBP4↓, UTP6↓, SNU13↓, TBL3↓, WDR36↓, UTP15↓, WDR3↓, UTP14A↓, RRP7A↓, IMP4↓, NXF1↓, POP1↓, RIOK1↓, GNL2↓, NAT10↓, RAN↓, NOP58↓, CIRH1A↓, GAR1↓, MDN1↓, UTP18↓, HEATR1↓, NXT1↓, DROSHA↓, BMS1↓, NHP2↓, REXO1↓, NVL↓, EIF6↓, XPO1↓, RBM28↓, POP5↓, RPP40↓, FBL↓, LSG1↓, POP7↓, RPP25L↓, POP4↓, RIOK2↓, IMP3↓, EFTUD1↓, WDR75↓, SPATA5↓, TAF9↓, GNL3L↓, MPHOSPH10↓, RPP30↓, SBDS↓, EMG1↓, NMD3
Protein processing in endoplasmic reticulum↑ 114 6.82E-06 PDIA4↑, MBTPS1↑, CALR↑, HSPA5↑, HSP90B1↑, PDIA3↑, HERPUD1↑, UBE2G1↑, SEC24D↑, LMAN1↑, AMFR↑, MARCH6↑, CAPN1↑, HYOU1↑, XBP1↑, RNF5↑, MAPK10↑, ERO1B↑, EIF2AK1↑, ERP29↑, UGGT1↑, ATXN3↑, ATF4↑, DDOST↑, OS9↑, MAN1A1↑, NFE2L2↑, P4HB↑, PDIA6↑, EDEM2↑, SSR2↑, UBE4B↑, UBE2D1↑, UBXN6↑, FBXO2↑, UBE2E1↑, GANAB↑, PRKCSH↑, MAN1B1↑, SEL1L↑, WFS1↑, DNAJC3↑, DAD1↑, YOD1↑, ATF6B↑, BCAP31↑, TRAM1↑, CAPN2↑, SEC23A↑, STT3B↑, TUSC3↑, SEC31A↑, UBE2J1↑, CUL1↑, RPN2↑, SEC24A↑, EIF2AK3↑, UFD1L↑, SKP1↑, ERLEC1↑, SEC61G↑, BAG1↓, DNAJA1↓, HSPA1A↓, HSPH1↓, HSP90AA1↓, UBE2G2↓, SEC23B↓, HSPA8↓, HSP90AB1↓, PLAA↓, BAG2↓, VCP↓, HSPBP1↓, DNAJC1↓, HSPA1B↓, CKAP4↓, EIF2S1↓, MAP3K5↓, UBE2J2↓, DNAJB2↓, RAD23A↓, BAX↓, UGGT2↓, PPP1R15A↓, DDIT3↓, EIF2AK2↓, UBE2D4↓, SSR3↓, DNAJB1↓, SSR4↓, STUB1↓, UBE2D2↓, HSPA6↓, ERN1↓, EDEM1↓, TRAF2↓, DNAJB12↓, UBE2E2↓, MBTPS2↓, SIL1↓, RPN1↓, NGLY1↓, RAD23B↓, HSPA4L↓, NSFL1C↓, SEC31B↓, VIMP↓, BCL2↓, PREB↓, SSR1↓, NPLOC4↓, RBX1↓, SEC61A2
Cell cycle↑ 88 1.06E-05 RBL2↑, TFDP2↑, TGFB2↑, CDKN1B↑, ABL1↑, GADD45B↑, CCNB2↑, SMAD4↑, CDC25C↑, SKP2↑, CDKN1C↑, HDAC2↑, BUB1↑, CDK2↑, CDC14A↑, ATM↑, ORC4↑, CDKN2C↑, PTTG1↑, WEE2↑, CDC14B↑, STAG2↑, CDC26↑, GADD45G↑, YWHAG↑, CUL1↑, ORC2↑, ANAPC2↑, SKP1↑, MYC↓, CDC16↓, CDC25B↓, TFDP1↓, ANAPC7↓, CDC6↓, CCND1↓, MCM4↓, CDKN1A↓, GADD45A↓, MCM3↓, STAG1↓, ANAPC5↓, PKMYT1↓, ORC1↓, CDC25A↓, MCM2↓, PRKDC↓, SFN↓, ATR↓, SMC1A↓, E2F4↓, ORC5↓, PCNA↓, CDK6↓, E2F1↓, CDC45↓, MCM5↓, CDC7↓, YWHAH↓, ORC6↓, RBL1↓, ORC3↓, GSK3B↓, CDK4↓, CCNA2↓, MAD1L1↓, CDK1↓, CCNE1↓, MAD2L2↓, E2F2↓, ANAPC11↓, CCNE2↓, CCNH↓, E2F5↓, RBX1↓, CDC23↓, MAD2L1↓, CREBBP↓, YWHAZ↓, BUB3↓, ANAPC4↓, CHEK2↓, CDC20↓, MCM7↓, E2F3↓, PLK1↓, TTK↓, TP53

The arrows ↑ and ↓ indicate up-regulation or down-regulation.

Construction of the ceRNA network

CeRNA network was built on the basis of the differential expression lncRNA, miRNA, and mRNA. As shown in Figure 7, 5 lncRNA nodes, 12 miRNA nodes, 383 mRNA nodes and 1596 edges were in the network.

Figure 7 The lncRNA-miRNA-mRNA Competing endogenous RNA network. light blue nodes represent mRNAs, dark blue nodes represent lncRNAs and red nodes represent miRNAs.

PPI & functional analyses of mRNAs

The top 464 differentially expressed mRNAs were input into PPI network by Omicsbean (Figure 8A). Figure 8B shows that IL6 plays an important role in malignant transformation of human embryonic esophageal epithelial cells by HPV18 E6E7, It is the hub and connects multiple classic cancer pathways. Survival analysis was performed IL6 maybe serve as prognostic biomarkers (Figure 8C).

Figure 8 PPI network and Kaplan–Meier curve analysis. (A) PPI Network of differentially expressed mRNA; (B) the PPI Network of IL6 as a hub nodes; (C) Kaplan-Meier survival curves of IL6 in esophageal cancer.

Real-time QPCR verification

To verify the results of the RNA-Seq, we selected seven genes to corroborate the RNA-seq data using quantitative qRT-PCR. The 26th generation of SHEE cells (P26) was used as reference sample. The result shows qRT-PCR data have a strong correlation with the RNA-Seq data (Figure 9) and confirms that RNA-Seq results are reliable in this study.

Figure 9 Correlation between the results of qRT-PCR and RNA-Seq. Black bars represent RNA-Seq measurements, and red bars represent qRT-PCR experimental results, calculated using the 2−ΔΔCt method and normalized to GAPDH.

Discussion

HPV is a small double-stranded DNA virus that can cause squamous epithelium proliferation in human skin mucosa and is closely related to the occurrence of numerous tumours. HPV is reported to be a necessary causative factor for cervical and anal-genital tumours (such as anal cancer, penile cancer and vulvar cancer) (16). However, the mechanism of HPV in cell transformation is still unclear. The SHEE cell line, which was established by this paper’s authors, consists of immortalised epithelial cells from the embryonic oesophagus after being induced by the E6E7 genes of HPV 18. We previously demonstrated that carcinogenesis of esophageal epithelial cells induced by HPV was a multistage process that went through the initial, immortal, premalignant and malignant transformation stages (14,17-19). Therefore, this cell line may be useful in elucidating the carcinogenic mechanism of HPV. Considering the significant difference in differentiation between the 26th and 79th generations of the cells, we conducted a comparative transcriptomic study on these two generations of cells to analyse the global gene expression and the signal transduction pathways, hence exploring the molecular mechanism of malignant transformation of human embryonic esophageal epithelial cells by HPV18 E6E7. According to the results of the GO analysis (Figure 5), the most significant difference was ribosome biogenesis and ribonucleoprotein complex biogenesis. Ribosomes are complex organelles composed of different proteins and nucleic acids and are responsible for protein synthesis in living cells. Therefore, ribosome biogenesis is a very important biological process in cells. Dysregulated ribosome biogenesis plays an important role in the development and progression of most cancers. Upregulated ribosome biogenesis lead to downregulate the expression and activity of p53, thus promote tumor development. Coincidentally, p53 is also the target of E6 protein. At the same time, some studies have shown that the inhibition of ribosome biogenesis may provide a viable treatment for cancer treatment (20-23). The ribonucleoprotein complex is formed by the interaction of RNA and RNA-binding proteins in the cells and participates in almost all cellular life functions. The ribonucleoprotein complex, such as ribonucleic acid-binding motif protein 8A (RBM8A), has been shown to be involved in transcriptional translation, gene expression, cell cycle regulation and the apoptosis regulation of RNA, and its increased expression may play a carcinogenic role (24,25). Therefore, the imbalance in the biogenesis of the ribonucleoprotein complex plays an important role in cell malignant transformation.

Our analysis of DEGs in the current study demonstrated that the most significant genes involved MAGEA3 which is the most significant up-regulated gene and LINC00958 which is the most significant down-regulated gene. MAGEA3, as a member of the melanoma antigen gene family, has been reported as having a significant inhibitory effect on apoptosis. MAGEA3 will enhance the activity of specific E3 ubiquitin-synthase and AMP-activated protein kinase α1 (AMPKα1) to reduce AMP-activated protein kinase (AMPK) and plays a role in promoting the growth of cancer cells and inhibiting apoptosis (26). Abikhair et al. reported that MAGEA3 is only expressed in cancer and testis, so MAGEA3 has the potential to be a diagnostic marker (27). LINC00958 is a significant downregulated long-chain noncoding RNA in our analysis result. Few studies have been done on LINC00958, but it has been reported that cell viability is observed by silencing LINC00958 in bladder cancer. The migration ability was reduced, and LINC00958 was shown as being combined with the initiation of translation and RNA post-transcriptionally modified proteins and their involvement in regulation. The same results were obtained in a study of LINC00958 in glioma (28-30).

In this study, we used bioinformatics to analyze the expression and function of aberrantly expressed lncRNA, miRNA, and mRNA and establish ceRNA network and PPI network to predict the regulation mechanism of DEGs. Through ceRNA network analysis, we predicted an interesting regulatory mechanism, TP53TG1- miR-6835- IGF2. TP53TG1 is a recently identified lncRNA, study have shown that TP53TG1 may play an important role of tumor suppressor gene or oncogene in different tumors (31). As we know, E6 protein can inactivate tumor suppressor p53 by inducing its degradation (32). The association between TP53TG1 and p53 has been confirmed by CHIP (33,34). Insulin-like growth factor 2 (IGF2) is a mitogenic peptide hormone expressed by liver and many other tissues. IGF2 is highly expressed in many cancers and promotes the development of cancer. IGF is associated with multiple signaling pathways, such as MAPK signaling pathway, Ras signaling pathway and PI3K-AKT signaling pathway (35-37). These signaling pathways are closely related to the development of cancer. IGF2 has been proved to promote proliferation and migration of esophageal cancer. Meanwhile, IGF2 is also an excellent esophageal cancer therapeutic targets and biomarkers (38,39). Therefore, we boldly speculate TP53TG1 may contribute to the malignant transformation of SHEEs by acting as a competing endogenous RNA to competitively bind to miR-6835 and regulate IGF2 expression.

IL6 is the most critical node in the PPI network. IL-6 is one of the most characteristic proinflammatory cytokines that activate JAK-STAT signaling pathway, MAPK signaling pathway, PI3K-Akt signaling pathway and other pathways. In the PI3K-Akt signaling pathway, IL-6 can promote the proliferation of cancer cells and inhibit apoptosis in vitro and in vivo. In the JAK-STAT pathway, IL6 functions by regulating gp130, STAT3 is considered to be the main signal transduction factor downstream of the gp130 signal transduction, acting as an oncogene and key factor in the link between inflammation and cancer (40-42). Simultaneously, Kaplan- Meier survival curves has shown IL6 maybe serve as prognostic biomarkers and therapy target in HPV-induced esophageal cancer.


Conclusions

Using RNA-Seq, we identified DEGs between the 26th and 79th generations of SHEEs and found that the significantly DEGs were involved in the development of cancer. Our results indicate that the activations of multiple biological processes, cellular components and molecular function are involved in the dynamic progressive process of carcinogenesis of esophageal epithelial cells induced by HPV. Furthermore, we predicted TP53TG1-miR-6835-IGF2 and IL6 play an important role in HPV-induced malignant transformation of human embryonic esophageal epithelial cells.


Acknowledgments

Funding: This work is supported by the National Natural Science Foundation of China (No.21677006).


Footnote

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at http://dx.doi.org/10.21037/tcr.2020.02.23). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). Institutional ethical approval and informed consent were waived.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.


References

  1. Torre LA, Bray F, Siegel RL, et al. Global cancer statistics, 2012. CA Cancer J Clin 2015;65:87-108. [Crossref] [PubMed]
  2. Stoner GD, Gupta A. Etiology and chemoprevention of esophageal squamous cell carcinoma. Carcinogenesis 2001;22:1737-46. [Crossref] [PubMed]
  3. Poljak M, Sterbenc A, Lunar MM. Prevention of human papillomavirus (HPV)-related tumors in people living with human immunodeficiency virus (HIV). Expert Rev Anti Infect Ther 2017;15:987-99. [Crossref] [PubMed]
  4. Zhou ZX, Li D, Guan SS, et al. Immunotherapeutic Effects of Dendritic Cells Pulsed with a Coden-optimized HPV 16 E6 and E7 Fusion Gene in Vivo and in Vitro. Asian Pac J Cancer Prev 2015;16:3843-7. [Crossref] [PubMed]
  5. Zhou ZX, Zhao C, Li QQ, et al. A novel mutant of human papillomavirus type 18 E6E7 fusion gene and its transforming activity. Asian Pac J Cancer Prev 2014;15:7395-9. [Crossref] [PubMed]
  6. Zhang K, Li JT, Li SY, et al. Integration of human papillomavirus 18 DNA in esophageal carcinoma 109 cells. World J Gastroenterol 2011;17:4242-6. [Crossref] [PubMed]
  7. Syrjanen K, Parkkinen S, Mantyjarvi R, et al. Human papillomavirus (HPV) type as an important determinant of the natural history of HPV infections in uterine cervix. Eur J Epidemiol 1985;1:180-7. [Crossref] [PubMed]
  8. Cao B, Tian X, Li Y, et al. LMP7/TAP2 gene polymorphisms and HPV infection in esophageal carcinoma patients from a high incidence area in China. Carcinogenesis 2005;26:1280-4. [Crossref] [PubMed]
  9. Han C, Qiao G, Hubbert NL, et al. Serologic association between human papillomavirus type 16 infection and esophageal cancer in Shaanxi Province, China. J Natl Cancer Inst 1996;88:1467-71. [Crossref] [PubMed]
  10. Farhadi M, Tahmasebi Z, Merat S, et al. Human papillomavirus in squamous cell carcinoma of esophagus in a high-risk population. World J Gastroenterol 2005;11:1200-3. [Crossref] [PubMed]
  11. Weston AC, Prolla JC. Association between esophageal squamous cell carcinoma and human papillomavirus detected by Hybrid Capture II assay. Dis Esophagus 2003;16:224-8. [Crossref] [PubMed]
  12. Acevedo-Nuño E, Gonzalez-Ojeda A, Vazquez-Camacho G, et al. Human papillomavirus DNA and protein in tissue samples of oesophageal cancer, Barrett's oesophagus and oesophagitis. Anticancer Res 2004;24:1319-23. [PubMed]
  13. Monsonego J. Contracept Fertil Sex 1995;23:731-40. [Cellular and molecular pathogenesis of cancer of the cervix]. [PubMed]
  14. Shen ZY, Xu LY, Chen MH, et al. Progressive transformation of immortalized esophageal epithelial cells. World J Gastroenterol 2002;8:976-81. [Crossref] [PubMed]
  15. Shen Z, Cen S, Shen J, et al. Study of immortalization and malignant transformation of human embryonic esophageal epithelial cells induced by HPV18 E6E7. J Cancer Res Clin Oncol 2000;126:589-94. [Crossref] [PubMed]
  16. Herrero R, Castellsague X, Pawlita M, et al. Human papillomavirus and oral cancer: the International Agency for Research on Cancer multicenter study. J Natl Cancer Inst 2003;95:1772-83. [Crossref] [PubMed]
  17. Shen ZY, Xu LY, Li C, et al. A comparative study of telomerase activity and malignant phenotype in multistage carcinogenesis of esophageal epithelial cells induced by human papillomavirus. Int J Mol Med 2001;8:633-9. [PubMed]
  18. Shen ZY, Xu LY, Li EM, et al. The multistage process of carcinogenesis in human esophageal epithelial cells induced by human papillomavirus. Oncol Rep 2004;11:647-54. [PubMed]
  19. Shen Z, Shen J, Zeng Y. Biological characteristics of human fetal esophageal epithelial cell line immortalized by the E6 and E7 gene of HPV type 18. Zhonghua Shi Yan He Lin Chuang Bing Du Xue Za Zhi 1999;13:209-12. [PubMed]
  20. Derenzini M, Montanaro L, Trere D. Ribosome biogenesis and cancer. Acta Histochem 2017;119:190-7. [Crossref] [PubMed]
  21. Pelletier J, Thomas G, Volarevic S. Ribosome biogenesis in cancer: new players and therapeutic avenues. Nat Rev Cancer 2018;18:51-63. [Crossref] [PubMed]
  22. Penzo M, Montanaro L, Trere D, et al. The Ribosome Biogenesis-Cancer Connection. Cells 2019;8: [Crossref] [PubMed]
  23. Bianco C, Mohr I. Ribosome biogenesis restricts innate immune responses to virus infection and DNA. Elife 2019;8: [Crossref] [PubMed]
  24. Ishigaki Y, Nakamura Y, Tatsuno T, et al. Depletion of RNA-binding protein RBM8A (Y14) causes cell cycle deficiency and apoptosis in human cells. Exp Biol Med (Maywood) 2013;238:889-97. [Crossref] [PubMed]
  25. Liang R, Lin Y, Ye JZ, et al. High expression of RBM8A predicts poor patient prognosis and promotes tumor progression in hepatocellular carcinoma. Oncol Rep 2017;37:2167-76. [Crossref] [PubMed]
  26. Wu F, Liu F, Dong L, et al. miR-1273g silences MAGEA3/6 to inhibit human colorectal cancer cell growth via activation of AMPK signaling. Cancer Lett 2018;435:1-9. [Crossref] [PubMed]
  27. Abikhair M, Roudiani N, Mitsui H, et al. MAGEA3 Expression in Cutaneous Squamous Cell Carcinoma Is Associated with Advanced Tumor Stage and Poor Prognosis. J Invest Dermatol 2017;137:775-8. [Crossref] [PubMed]
  28. Guo E, Liang C, He X, et al. Long Noncoding RNA LINC00958 Accelerates Gliomagenesis Through Regulating miR-203/CDK2. DNA Cell Biol 2018;37:465-72. [Crossref] [PubMed]
  29. Seitz AK, Christensen LL, Christensen E, et al. Profiling of long non-coding RNAs identifies LINC00958 and LINC01296 as candidate oncogenes in bladder cancer. Sci Rep 2017;7:395. [Crossref] [PubMed]
  30. He W, Zhong G, Jiang N, et al. Long noncoding RNA BLACAT2 promotes bladder cancer-associated lymphangiogenesis and lymphatic metastasis. J Clin Invest 2018;128:861-75. [Crossref] [PubMed]
  31. Zhang Y, Yang H, Du Y, et al. Long noncoding RNA TP53TG1 promotes pancreatic ductal adenocarcinoma development by acting as a molecular sponge of microRNA-96. Cancer Sci 2019;110:2760-72. [Crossref] [PubMed]
  32. Li X, Coffino P. High-risk human papillomavirus E6 protein has two distinct binding sites within p53, of which only one determines degradation. J Virol 1996;70:4509-16. [Crossref] [PubMed]
  33. Diaz-Lagares A, Crujeiras AB, Lopez-Serra P, et al. Epigenetic inactivation of the p53-induced long noncoding RNA TP53 target 1 in human cancer. Proc Natl Acad Sci U S A 2016;113:E7535-44. [Crossref] [PubMed]
  34. Nikulenkov F, Spinnler C, Li H, et al. Insights into p53 transcriptional function via genome-wide chromatin occupancy and gene expression analysis. Cell Death Differ 2012;19:1992-2002. [Crossref] [PubMed]
  35. Pereira SS, Monteiro MP, Costa MM, et al. IGF2 role in adrenocortical carcinoma biology. Endocrine 2019;66:326-37. [Crossref] [PubMed]
  36. Cacheux W, Lievre A, Richon S, et al. Interaction between IGF2-PI3K axis and cancer-associated-fibroblasts promotes anal squamous carcinogenesis. Int J Cancer 2019;145:1852-9. [Crossref] [PubMed]
  37. Ji Y, Wang Z, Li Z, et al. Silencing IGF-II impairs C-myc and N-ras expressions of SMMC-7721 cells via suppressing FAK/PI3K/Akt signaling pathway. Cytokine 2017;90:44-53. [Crossref] [PubMed]
  38. Gao T, He B, Pan Y, et al. Long non-coding RNA 91H contributes to the occurrence and progression of esophageal squamous cell carcinoma by inhibiting IGF2 expression. Mol Carcinog 2015;54:359-67. [Crossref] [PubMed]
  39. Zheng DN, Zhang CJ, Sun GP. Long non-coding RNA MNX1-AS1 promotes migration and invasion of esophageal squamous cell carcinoma by upregulating IGF2. Eur Rev Med Pharmacol Sci 2019;23:6179-85. [PubMed]
  40. Landskron G, De la Fuente M, Thuwajit P, et al. Chronic inflammation and cytokines in the tumor microenvironment. J Immunol Res 2014;2014:149185.
  41. Taniguchi K, Karin M. IL-6 and related cytokines as the critical lynchpins between inflammation and cancer. Semin Immunol 2014;26:54-74. [Crossref] [PubMed]
  42. Nguyen DP, Li J, Tewari AK. Inflammation and prostate cancer: the role of interleukin 6 (IL-6). BJU Int 2014;113:986-92. [Crossref] [PubMed]
Cite this article as: Tang D, Wang B, Khodahemmati S, Li J, Zhou Z, Gao J, Sheng W, Zeng Y. A transcriptomic analysis of malignant transformation of human embryonic esophageal epithelial cells by HPV18 E6E7. Transl Cancer Res 2020;9(3):1818-1832. doi: 10.21037/tcr.2020.02.23

Download Citation