DUXAP8 a Pan-Cancer Prognostic Marker Involved in the Molecular Regulatory Mechanism in Hepatocellular Carcinoma: A Comprehensive Study Based on Data Mining, Bioinformatics, and in vitro Validation

Introduction

Despite the rapid advances in medical science in recent decades, cancer still represents the most significant threat to human health.^1–3 Malignant tumors remain one of the leading causes of death worldwide, mainly due to the late diagnosis of most patients who are frequently not diagnosed until advanced stages as well as to the lack of effective therapy.⁴ Accordingly, the identification of novel cancer biomarkers and the understanding of their involvement in the underlying molecular mechanisms of tumorigenesis and cancer progression are urgent and essential for the development of more effective diagnostic methods and therapeutic strategies.

Pseudogenes have long been considered to be non-functional transcriptional relics of the genome of various organisms.⁵ Approximately 15,000 pseudogenes, including unprocessed and processed pseudogenes, have been identified in the human genome according to the latest annotation of the ENCODE project.⁶ The functional role of pseudogenes remains poorly understood and several studies have suggested that pseudogenes could play key roles in cancer progression by acting as competing endogenous RNAs (ceRNAs) in recent years.^7,8

Double homeobox A pseudogene 8 (DUXAP8) is a member of the DUXA homeobox gene family. In 2017, Sun et al first determined by bioinformatics analysis that DUXAP8 was overexpressed in lung cancer tissues and then demonstrated that DUXAP8 promoted cell invasion and proliferation through epigenetic silencing of RHOB and EGR1 in non-small cell lung cancer.⁹ Recently, pseudogene DUXAP8 has been shown to be upregulated in multiple cancers and plays critical regulatory roles in cancer progression. For instance, DUXAP8 was shown to regulate cell growth by epigenetic silencing of PLEKHO1 expression in gastric cancer.¹⁰ In addition, knockdown of DUXAP8 inhibited cell proliferation and invasion in esophageal squamous cell cancer.¹¹ Moreover, expression of DUXAP8 promoted pancreatic carcinoma cell growth by epigenetic silencing of KLF2 and CDKN1A.¹² Also, high levels of DUXAP8 was associated with poor prognosis in patients with lung,⁹ gastric,¹⁰ esophageal,¹¹ pancreatic,¹² glioma,¹³ and bladder¹⁴ cancers. However, the degree and reliability of the prognostic value of DUXAP8 in human cancers has not been thoroughly evaluated due to the limited studies and small sample size. Furthermore, the expression pattern and biological function of DUXAP8 in HCC remain unclear.

In this study, due to the limited published studies on DUXAP8, RNA sequencing (RNA-Seq) count matrix and matched clinical data extracted from the Cancer Genome Atlas (TCGA) database were used to evaluate the expression of DUXAP8 and its clinical significance in human cancers. Additionally, we verified the DUXAP8 expression levels in 66 HCC samples and corresponding adjacent normal liver tissue samples by polymerase chain reaction (PCR) analysis and explored its biological function in HCC through in vitro cellular function assays.

Materials and Methods

Data Acquisition

Data from TCGA database (https://cancergenome.nih.gov/), RNA-Seq count matrices of human pan-cancers were extracted using the R software (https://www.r-project.org/). Datasets that lacked sufficient samples or survival information were excluded. Gene symbols corresponding to each Ensembl ID were acquired based on the gene transfer format file (Homo_sapiens.GRCh38.96.chr) using the R software. In addition, studies concerning DUXAP8 in the PubMed, EMBASE and Web of Science databases were comprehensively evaluated in this study.

DUXAP8 Expression and Its Clinical Value in Human Cancers

Expression levels of DUXAP8 in various human cancers were acquired from TCGA database to assess the expression patterns of DUXAP8 in pan-cancers. The area under the curve (AUC), specificity, and sensitivity of each dataset were acquired by receiver operating characteristic (ROC) curve analysis. Then, the pooled AUC was calculated by summary ROC (SROC) curve analysis to evaluate the diagnostic value of DUXAP8 in pan-cancers. The Deeks’ funnel plot asymmetrical test was performed to assess the potential publication bias in the SROC analysis.

The hazard ratio (HR) with 95% confidence interval (CI) was calculated to evaluate the association of DUXAP8 expression with overall survival (OS) in each study, and the pooled HR was calculated by meta-analysis to assess the prognostic importance of DUXAP8 in pan-cancers. Potential publication bias was assessed using the funnel plots with Begg’s test. High DUXAP8 expression was considered to be related to poor prognosis if the pooled HR>1 and the 95% CI did not overlap 1.

Tissue Specimens and HCC Cells

A total of 66 pairs of HCC-tissue and adjacent liver-tissue samples were collected from the Department of General Surgery at Beijing Tongren Hospital between 2011 and 2017. No patients received other treatments before sample collection. Written informed consent, which was in accordance with the Declaration of Helsinki, was obtained from each patient. The protocol of this study was approved by the institutional review board of Beijing Tongren Hospital, Capital Medical University.

HCC cell lines (BEL-7404, HEP-G2, BEL-7402, SMMC-7721, HUH-7, and HEP-3B) and the normal human liver cell line (LO-2) were obtained from the Cell Bank of Shanghai Institute of Biochemistry & Cell Biology (Shanghai, China). All the cells were cultured in Dulbecco’s modified Eagle’s medium containing 10% fetal bovine serum.

Real-Time Quantitative PCR (RT-qPCR) Analysis

Total RNA was extracted from the tissue samples using Trizol (Invitrogen/Thermo Fisher Scientific, Waltham, MA, USA). A reverse-transcription kit was used to synthesize cDNA. The real-time PCR analysis was performed using a real-time PCR system (Applied Biosystems, Foster City, CA, USA) according to the manufacturer’s instructions. Gene expression level was determined by comparing Ct (2^−ΔΔCT) values.¹⁵

Cell Transfection

Three siRNA sequences (si-DUXAP8-1#: 5ʹ- AAGATAAAGGTGGTTTCCACAAGAA-3ʹ; si-DUXAP8-2#: 5ʹ-GGAACTTCCCAAACCTCCATGATTT-3ʹ; si-DUXAP8-3#: 5ʹ-CAGCATACTTCAAATTCACAGCAAA-3ʹ) were designed according to the cDNA sequence of DUXAP8 obtained from GeneBank. HCC cells were transfected with siRNA using lipofectamine (Invitrogen) according to the manufacturer’s instructions. In addition, a negative control sequence (si-NC) was synthesized for use as the negative control group.

Cell Viability and Colony Formation Assays

Cells from si-NC and si-DUXAP8 groups were seeded into 96-well plates. Cells viability was accessed at 24, 48, 72, and 96 h using the Cell Counting Kit 8 (CCK-8) assays according to the manufacturer’s instructions (Dojindo Laboratories, Kumamoto, Japan). For the colony formation assay, the colonies containing more than 50 cells were counted two weeks after seeding the cells into 6-well plates.

Transwell Assay and Wound Healing assay

Transwell chambers with Matrigel were used to evaluate cell invasion ability. Briefly, cell suspensions of the si-DUXAP8 and si-NC groups were added to the upper Transwell chambers containing serum-free medium. The corresponding lower chambers containing 10% FBS medium served as a chemo-attractant. Cells that migrated onto the lower surface of the membrane after 24 h incubation were fixed with formaldehyde and stained with hematoxylin for 30 min. Cells which migrated to the lower compartment of the chamber were counted using a microscope in 5 randomly selected fields. For the wound healing assay, cells from the 2 groups were seeded into 6-well plates. A scratch wound was made in each well by scratching the cell layer once confluence was reached. Cells were cultured in medium after washing away the detached cells with PBS. Photos were captured at 0, 8 and 24 h to determine the number of cells migrating into the scratch wound.

Statistical Analysis

Kaplan-Meier survival analysis and ROC curve analysis were performed using the GraphPad Prism software (GraphPad Software Inc., La Jolla, CA, USA). Survival data was extracted from the Kaplan–Meier curve using the Engauge Digitizer software if the published article did not provide them, and then the HR was calculated using Tierney’s method.¹⁶ Meta-analysis and SROC curve analysis were conducted using the Stata software (College Station, TX, USA). All experiments in the current study were conducted in triplicate, and the results are represented as the mean ± SD. Comparison among different groups was performed by one-way analysis of variance (ANOVA) test or Student’s t-test, and P<0.05 indicates statistical significance.

Results

DUXAP8 Was Overexpressed and Had Excellent Clinical Value in Human Cancers

A total of 28 studies including 8301 patients were enrolled in this study. Among them, 17 TCGA datasets contained the expression data of DUXAP8 in sufficient tumor and normal samples. HRs with 95% CI for OS of cancer patients were calculated based on DUXAP8 expression and survival information extracted from 22 TCGA datasets. Six published articles detected OS comparisons between patients with high and low expression of DUXAP8. Detailed information of the collected studies is shown in Tables S1 and S2.

The data of DUXAP8 expression in human cancer -tissue samples and normal tissue samples from patients with various types of cancers were extracted from each RNA-Seq count matrix and then log₂-transformed. DUXAP8 was significantly overexpressed in 15 types of human cancers, as shown in Figure 1. The AUC was obtained by ROC curve analysis for each TCGA dataset (Table S1) to determine the diagnostic value of DUXAP8 in each type of cancer. Then, the pooled specificity, sensitivity, and AUC were obtained by SROC curve analysis. As shown in Figure 2A and B, the pooled AUC was 0.82 [0.78–0.85], pooled sensitivity was 0.70 [0.63–0.75], and the pooled specificity was 0.84 [0.76–0.90]. The above results indicated that DUXAP8 had an excellent diagnostic value in pan-cancers. The results of Deeks’ funnel plot asymmetry test indicated that no significant publication bias existed for the studies in this SROC analysis (P =0.25) (Figure 2C).

Figure 1 DUXAP8 expression levels in pan-cancers according to TCGA database (*P<0.05). Abbreviation: TCGA, The Cancer Genome Atlas.

Figure 2 SROC curve analysis according to TCGA database. Notes: (A) Pooled sensitivity and pooled specificity were calculated to assess the diagnostic value of DUXAP8 in pan-cancers. (B) The results of pooled AUC (0.82 [0.78–0.85]) indicated that DUXAP8 has a good diagnostic value in pan-cancers. (C) No publication bias existed in the SROC analysis based on Deeks’ funnel plot asymmetrical test (P=0.25). Abbreviations: SROC, summary receiver operating characteristic curve; TCGA, the cancer genome atlas; AUC, area under the curve.

Next, we calculated the HR of each TCGA dataset based on the DUXAP8 expression data and matched OS information of the patients. For the published articles that were included in this study, the HR was obtained using the survival data excerpted from the Kaplan–Meier curves. The pooled HR obtained by meta-analysis suggested that DUXAP8 had an excellent prognostic value in human cancers (pooled HR=1.20[1.06–1.34], P<0.01, random), as shown in Figure 3. No significant publication bias existed for the studies in this meta-analysis according to the Begg’s funnel plot test (P =0.18). The results of this meta-analysis were reliable according to the sensitivity analysis, which indicated that all the datasets were within the 95% CI of the combined analysis.

Figure 3 Pooled HR was calculated to assess the prognostic value of DUXAP8 in pan-cancers. Note: Meta-analysis was performed to obtain the pooled HR. Abbreviation: HR, hazard ratio.

The above results suggested that DUXAP8 was upregulated in human tumor tissues and could be used as a prognostic and diagnostic marker in human cancers.

Correlation Between DUXAP8 Expression and Clinicopathological Characteristics in HCC Patients

Given that the expression pattern and clinical significance of DUXAP8 in HCC were not yet characterized, we investigated the correlation between the DUXAP8 expression level of HCC tissue and the clinicopathologic characteristics of HCC patients based on the information extracted from TCGA database. The results revealed that DUXAP8 was significantly overexpressed in HCC tissue (Figure 4A and B), and DUXAP8 had excellent diagnostic value in HCC patients, according to the ROC curve analysis (AUC=0.84[0.80–0.88]) (Figure 4C). In addition, a high DUXAP8 expression level indicated poor OS in HCC patients (Figure 4D). Moreover, DUXAP8 was significantly overexpressed in patients over 60 years of age (Age≥ 60, 5.61 ± 0.15 vs. Age< 60, 4.73 ± 0.18; P< 0.01), with higher pathological stage (Stage III+IV, 5.67 ± 0.23 vs. StageI+II, 4.97 ± 0.14; P= 0.04), and with vascular invasion (Positive, 5.47 ± 0.23 vs. Negative 4.78 ± 0.15; P= 0.03), as shown in Table 1.

Table 1 Correlation Between DUXAP8 Expression and Clinicopathological Characteristics in TCGA

Figure 4 DUXAP8 expression patterns in HCC according to TCGA database. Notes: (A) DUXAP8 expression levels in HCC tissues compared with normal liver tissues, according to unpaired student’s t-test. (B) DUXAP8 expression level in HCC tissue compared with adjacent liver tissue from the same patient according to paired student’s t-test. (C) Significance of DUXAP8 in HCC according to the ROC curve analysis. (D) Overall survival (OS) plots of DUXAP8. Abbreviations: HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; ROC, receiver operating characteristic curve.

Identification of DUXAP8 Related Protein-Coding Genes in HCC

To investigate the potential molecular regulatory mechanisms of DUXAP8 in HCC tumorigenesis, the protein-coding genes whose expression level was closely correlated with that of DUXAP8 were identified. Two-sided Pearson correlation coefficient analysis and the z-test were performed using the R software based on the gene expression data extracted from TCGA. |Pearson correlations| >0.50 and z-test p<0.001 were used as cutoff criteria. Ultimately, a total of 177 protein-coding genes were identified as DUXAP8-related genes (Table 2).

Table 2 DUXAP8-Related Protein-Coding Genes in HCC Based on TCGA Database

Gene-Annotation and Pathway Enrichment Analysis

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and gene ontology (GO) annotation enrichment analyses were performed using the DAVID online tool (https://david.ncifcrf.gov/) to identify the potential signaling pathways and biological processes involved in HCC tumorigenesis and progression in which DUXAP8 participates.

As shown in Figure 5A, regarding the GO enrichment analysis, there were 25 significant biological process annotations (P<0.01), which included G1/S and transition of mitotic cell cycle, cell division, G2/M transition of mitotic cell cycle, cell proliferation, etc. The GO enrichment analysis also revealed that there were 22 significant cellular component terms (P<0.01), which were correlated with condensed chromosome kinetochore, kinetochore, nucleoplasm, centrosome, etc. The GO analysis also revealed 8 significant molecular function terms (P<0.01), associated with protein binding, microtubule motor activity, ATP binding, DNA-dependent ATPase activity, etc. The detailed data are listed in Table 3.

Table 3 The GO Analysis of DUXAP8-Related Protein-Coding Genes in HCC Based on TCGA Database

Figure 5 Gene-annotation and pathway enrichment analysis. Notes: (A) GO enrichment analysis; (B) KEGG pathway enrichment analysis.Abbreviations: GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Regarding the KEGG pathway enrichment analysis, there were 5 signaling pathways that were statistically significant (P<0.01), including cell cycle, Fanconi anemia pathway, homologous recombination, oocyte meiosis and progesterone-mediated oocyte maturation (Figure 5B).

DUXAP8 Was Overexpressed in HCC Tissue and Cell Lines

To verify the expression pattern of DUXAP8 in HCC, qPCR analysis was conducted to determine the expression levels of the 66 pairs of HCC tissue samples. As shown in Figure 6A, DUXAP8 was significantly overexpressed in HCC tissue, according to the paired Student’s t-test. Furthermore, high DUXAP8 expression was significantly associated with poor OS in HCC patients, according to the Kaplan–Meier curve analysis (Figure 6B), which was consistent with the result of the comprehensive analysis based on TCGA database.

Figure 6 DUXAP8 expression pattern in HCC was verified by qPCR analysis. Notes: (A) DUXAP8 expression levels in 66 pairs of HCC tissue and adjacent normal liver tissue samples. (B) Overall survival (OS) plots of DUXAP8 based on the 66 HCC patients’ clinical information. (C) DUXAP8 expression levels in normal liver cells and 6 HCC cell lines. Abbreviation: HCC, hepatocellular carcinoma.

In addition, the expression levels of DUXAP8 in the 6 HCC cell lines were higher than that in the normal liver LO-2 cells (Figure 6C). The transfection of the Hep-G2 cells with three siRNAs using a lentiviral vector system revealed that the siRNA-DUXAP8-3# (5ʹ-CAGCATACTTCAAATTCACAGCAAA-3ʹ) had the highest interference efficiency according to qPCR analysis. (Figure 7A), thus it was used in subsequent experiments.

Figure 7 Effects of DUXAP8 on HCC cell proliferation and invasion. Notes: (A) Interference efficiency of the 3 si-RNAs in Hep-G2 cells. (B, C) Downregulation of DUXAP8 inhibited HCC cell proliferation and viability according to the CCK-8 assay and colony formation assay, respectively. (D, E) Invasion and migration ability of HCC cells was decreased after DUXAP8 knockdown according to the Transwell assay and wound healing assay (*P<0.05). Abbreviations: HCC, hepatocellular carcinoma; CCK-8, Cell Counting Kit-8.

Effects of DUXAP8 Expression on HCC Cells Proliferation and Invasion

Cell viability in the DUXAP8 knockdown group was significantly inhibited compared with that in the si-NC group according to the growth curves obtained with the CCK-8 assay (Figure 7B). Additionally, according to the colony formation assay, HCC cells in the si-DUXAP8 group produced fewer colonies than those in the negative control group (Figure 7C).

The results of the Transwell assay suggested that downregulation of DUXAP8 inhibited the invasion ability of HCC cells (Figure 7D). The wound healing assay also revealed that the migration ability of HCC cells in the si-DUXAP8 group was reduced compared with that in the si-NC group (Figure 7E).

The results of the in vitro cellular function assays described above suggested that overexpression of DUXAP8 in HCC promoted cell proliferation and invasion.

Discussion

As a result of the development of next-generation sequencing technology, the integrated analysis through data mining, bioinformatics analysis, and meta-analysis has become one of the most powerful methods in cancer research.¹⁷ For example, Zhang et al found that the lncRNA MIR22HG could act as a cancer suppressor in HCC based on data extracted from TCGA portal.¹⁸ In addition, Ge et al showed that LINC00996 is significantly decreased in colorectal cancer tissues and the depletion of LINC00996 is associated with a poor outcome in colorectal cancer patients according to a meta-analysis using data extracted from the Gene Expression Omnibus (GEO) and TCGA databases.¹⁹

In our previous research, we used an integrated analysis approach based on the microarray data from the GEO and TCGA databases and found that DUXAP10, which is another member of the DUXA homeobox gene family, could act as a potential marker in pan-cancers.²⁰ Additionally, we also showed that DUXAP10 could regulate cell growth by activating p-AKT in HCC.²¹ After comprehensively searching the studies on DUXAP8, we only found a few published articles in electronic databases. Therefore, we decided to mine the expression data of DUXAP8 from public databases, including TCGA and GEO.

Although, DUXAP8 has been identified in previous studies as an aberrantly expressed gene in lung cancer and gastric cancer by bioinformatics analysis based on GEO datasets,^9,10 we could not find any data on DUXAP8 expression in human cancers from the GEO database, which might be due to the update of probe annotation files. Nevertheless, we successfully extracted the expression information of DUXAP8 from TCGA database. Comprehensive analysis including data mining, OS analysis, ROC curve analysis, SROC curve analysis, and meta-analysis were performed in this study, and the results revealed that DUXAP8 was overexpressed in human cancers tissue samples and had both excellent diagnostic and prognostic values in human cancers.

We also found that the DUXAP8 expression level in HCC tissue was significantly upregulated compared with that in adjacent normal liver tissue (log₂FC>1 and P<0.01) based on the data extracted from TCGA. Also, high DUXAP8 expression was highly related to poor OS in HCC (P<0.01), and had excellent diagnostic value (AUC=0.84, 95% CI 0.80–0.88) in HCC according to the ROC curve analysis. Moreover, high DUXAP8 levels were detected in HCC patients with vascular invasion or high pathologic stage according to our analysis. However, to date, no study has been reported on the roles of DUXAP8 in HCC tumorigenesis and progression. Accordingly, we identified the DUXAP8 related protein-coding genes to investigate the potential biological process and signaling pathway in which DUXAP8 is involved. The results of the GO and KEGG pathway enrichment analysis suggested that DUXAP8-related genes were most significantly involved in the biological processes of cell division, cell cycle, etc. Also, according to the above results, we decided to experimentally examine the roles of DUXAP8 in HCC. The results of the qPCR analysis and in vitro cell function assays showed that DUXAP8 was upregulated in HCC tissue, and its downregulation inhibited HCC cell proliferation and invasion.

This study had a number of limitations that should be discussed. First, no clinicopathological characteristics were closely related with the DUXAP8 expression levels in the 66 HCC tissue samples that were collected in this study, which might be explained by the limited sample size. Second, the signaling pathways through which DUXAP8 might exert its regulatory effect on HCC progression were not explored in the current study. Thus, undoubtedly, this issue needs to be addressed through further investigation and verification.

Conclusions

To summarize, this study revealed that DUXAP8 was overexpressed in pan-cancers and could serve as a potential prognostic and diagnostic marker for cancer patients. In addition, we verified the DUXAP8 expression pattern in HCC and showed that DUXAP8 regulated HCC cell growth. To the best of our knowledge, this study is the first investigation of the regulatory roles of DUXAP8 in HCC and may contribute to the development of novel diagnosis markers and therapeutic approach for HCC.