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

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. 2020 Translational Cancer Research. All rights reserved.

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 . 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 (). 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 (, ).

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 (). 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 (, ).

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 (). KEGG is a comprehensive database that integrates genomic, chemical and system function information. As shown in , the three most affected pathways are ribosome biogenesis in eukaryotes, protein processing in the endoplasmic reticulum and the cell cycle. 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.

CeRNA network was built on the basis of the differential expression lncRNA, miRNA, and mRNA. As shown in , 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 (). 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 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 () 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 (), 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.