Abstract

   Hepatocellular carcinoma (HCC) is one of the deadliest malignant tumors
   that are harmful to human health. Increasing evidence has underscored
   the critical role of the competitive endogenous RNA (ceRNA) regulatory
   networks among various human cancers. However, the complexity and
   behavior characteristics of the ceRNA network in HCC were still
   unclear. In this study, we aimed to clarify a phosphatase and tensin
   homolog (PTEN)-related ceRNA regulatory network and identify potential
   prognostic markers associated with HCC. The expression profiles of
   three RNAs (long non-coding RNAs [lncRNAs], microRNAs [miRNAs], and
   mRNAs) were extracted from The Cancer Genome Atlas (TCGA) database. The
   DLEU2L-hsa-miR-100-5p/ hsa-miR-99a-5p-TAOK1 ceRNA network related to
   the prognosis of HCC was obtained by performing bioinformatics
   analysis. Importantly, we identified the DLEU2L/TAOK1 axis in the ceRNA
   by using correlation analysis, and it appeared to become a clinical
   prognostic model by Cox regression analysis. Furthermore, methylation
   analyses suggested that the abnormal upregulation of the DLEU2L/TAOK1
   axis likely resulted from hypomethylation, and immune infiltration
   analysis showed that the DLEU2L/TAOK1 axis may have an impact on the
   changes in the tumor immune microenvironment and the development of
   HCC. In summary, the current study constructing a ceRNA-based
   DLEU2L/TAOK1 axis might be a novel important prognostic factor
   associated with the diagnosis and prognosis of HCC.

   Keywords: HCC, ceRNA network, DLEU2L/TAOK1 axis, prognosis

Graphical Abstract

   graphic file with name fx1.jpg
   [75]Open in a new tab
     __________________________________________________________________

   Shi et al. provide evidence that a ceRNA-based DLEU2L/TAOK1 axis has an
   impact on changes in the tumor immune microenvironment and might be a
   novel important prognostic factor associated with the diagnosis and
   prognosis of HCC.

Introduction

   Liver cancer, mainly hepatocellular carcinoma (HCC), is the sixth most
   common cancer type and the fourth deadly cause of cancer in the
   world.[76]^1 It has been estimated that, in 2018, there were 782,000
   deaths and 841,000 new cases of liver cancer in the world, with a
   higher incidence in men than in women in many parts of the
   world.[77]^2^,[78]^3 Statistics indicate that the 5-year survival for
   HCC is 18%, which shows that only 30%–40% of patients were diagnosed in
   the early stage.[79]^4 Curative-intent treatments for early stage HCC
   include surgical resection, ablation, liver transplantation (LT), and
   transarterial chemoembolization (TACE).[80]^5 However, the main
   characteristics of HCC include a high rate of recurrence and
   metastasis, a low detection rate of curable stages, and ineffective
   treatments. Moreover, the incidence and mortality of HCC still continue
   to increase in many countries, and it has been generally recognized as
   a cancer with poor prognosis and reduced survival. Hence, it is
   essential to explore useful prognostic biomarkers and/or therapeutic
   targets for HCC to improve our deficiencies in the diagnosis,
   prevention, and treatment of the disease.

   It is well known that publicly available large-scale cancer omics data,
   such as The Cancer Genome Atlas (TCGA), provide us with clinical
   information and molecular data on various cancer patients, which may be
   very helpful for us to find candidate tumor biomarkers.[81]^6 The
   widespread applications of high-throughput RNA sequencing (RNA-seq)
   resulted in the development of better noninvasive, repeatable,
   sensitive, and accurate tools used in early screening, diagnosis,
   evaluation, and monitoring of patients.

   Long non-coding RNAs (lncRNAs) are broadly defined as a subtype of
   ncRNAs with a length greater than 200 nt that have no or limited
   protein-coding ability.[82]^7 Several studies have shown that lncRNAs
   have direct and indirect regulatory effects in cancer biology.[83]8,
   [84]9, [85]10 It has been reported that these long molecules are
   frequently dysregulated in a variety of cancers, and in different types
   of cancers (including HCC), specific lncRNAs are associated with cancer
   recurrence, metastasis, and poor prognosis.[86]11, [87]12, [88]13,
   [89]14 Therefore, it is thought that the specific lncRNA biomarkers
   related to the prognosis and diagnosis of HCC are of great clinical
   significance.

   MicroRNAs (miRNAs), a kind of small single-stranded ncRNA, include
   19–25 nt and can bind to the 3′ untranslated region (UTR) of its target
   mRNA to inhibit gene degradation or translation.[90]15, [91]16, [92]17
   We have learned that miRNAs regulate approximately 30% of the genes in
   the human genome.[93]^18 Moreover, miRNA-mediated post-transcriptional
   regulation requires multiple RNA-binding proteins, which are conducive
   to the function of miRNA in tumorigenesis.[94]19, [95]20, [96]21 In the
   cytoplasm, Salmena et al.[97]^22 proposed the competitive endogenous
   RNA (ceRNA) hypothesis that lncRNA is described as mainly regulating
   mRNA through a ceRNA regulation mechanism. As far as HCC is concerned,
   for the selected example, Zhan et al.[98]^23 demonstrated that
   HOXA11-AS, as an oncogene, can accelerate HCC growth through the
   miR-506-3p/Slug axis. Zhang et al.[99]^24 identified that high levels
   of NEAT1 promoted the metastasis of HCC cells by acting as a miR-320a
   molecular sponge and targeting LAGE3. However, a study by Xu
   et al.[100]^25 illustrated that overexpression of WWOX inhibited cell
   proliferation and metastasis by sponge miR-20b-5p.

   Phosphatase and tensin homolog (PTEN) is a very commonly mutated tumor
   suppressor gene in human cancers, including HCC.[101]26, [102]27,
   [103]28, [104]29 Accumulating evidence has shed light on PTEN activity,
   which can be regulated by mutations, epigenetic silencing, abnormal
   protein, transcriptional inhibition, localization, and
   posttranslational modification.[105]30, [106]31, [107]32 At the same
   time, some data have indicated that decreased PTEN expression is a
   common event in HCC and is associated with disease stage, tumor grade
   and size, and increased expression of tumor markers.[108]33, [109]34,
   [110]35, [111]36

   In this study, we constructed a ceRNA network related to PTEN in HCC
   ([112]Figure 1). First, through differential expression analysis in two
   groups of PTEN^high and PTEN^low expression in 371 HCC samples, the
   HCC-related lncRNA-miRNA-mRNA triple regulatory networks constructed
   from three differentially expressed RNAs were obtained. Functional
   enrichment analysis was conducted to assess the functional role and
   potential mechanism of the network in HCC. Then, a key ceRNA network
   was identified by expression analysis, survival analysis, and
   nuclear-cytoplasmic localization analysis of RNAs from hub triple
   regulatory networks. Moreover, we performed a correlation analysis
   between these four predictive genes and PTEN that showed that the
   DLEU2L/TAOK1 axis played a vital role in HCC. Finally, Cox regression
   analysis was carried out to obtain the diagnostic and prognostic values
   of TAOK1 in HCC, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and
   Genomes (KEGG) analyses were utilized to get the possible function of
   TAOK1 in HCC, and methylation analysis and immune infiltration analysis
   were further performed to study the potential biological function of
   TAOK1 in HCC.

Figure 1.

   [113]Figure 1
   [114]Open in a new tab

   Flowchart of construction and analysis of ceRNA

Results

The tumor suppressor role and prognostic value of PTEN overexpression in HCC

   To investigate the possible role of PTEN in HCC, based on the Human
   Protein Atlas (HPA) database, we found that PTEN was overexpressed in
   normal hepatic tissue, but downexpressed in HCC tissues ([115]Figures
   2A and [116]S1). Similar deregulation of PTEN expression was also
   demonstrated by immunohistochemistry (IHC) staining obtained from the
   HPA ([117]Figure 2B), and the patient data from the IHC are listed in
   [118]Table S1.

Figure 2.

   [119]Figure 2
   [120]Open in a new tab

   The tumor suppressor role of PTEN in human hepatocellular carcinoma
   (HCC)

   (A) Expression distribution of PTEN in pan-cancer tissues. (B)
   Validation of the expression of PTEN on the translational level by the
   Human Protein Atlas database (immunohistochemistry). (C) The low (n =
   185) and high expression (n = 184) of PTEN were compared using a
   Kaplan-Meier survival curve. (D) The distribution of PTEN genomic
   alterations in TCGA HCC is shown on a cBioPortal OncoPrint plot. (E and
   F) The association between PTEN copy number and mRNA expression are
   shown in the dot plot (E) and correlation plot (F) by cBioPortal.

   Since PTEN is abnormally downexpressed in tumor specimens, we then
   studied the clinical significance of PTEN expression in HCC patients.
   According to the Kaplan-Meier survival curves, as shown in
   [121]Figure 2C, our data indicated that lower expression of PTEN
   markedly correlated with poor overall survival (OS) in patients of the
   HCC cohort.

   Furthermore, to understand the potential mechanism of the abnormally
   low expression of PTEN in HCC, we analyzed the genomic and copy numbers
   of PTEN. The results of our analysis on cBioPortal are as follows, with
   an OncoPrint plot showing the deletion of PTEN gene in TCGA HCC dataset
   ([122]Figure 2D).

   As shown in [123]Figure 1E, more than one-third of all HCC samples
   harbored PTEN deletion, and, consistently, HCC samples harboring PTEN
   deletion exhibited lower mRNA expression than did those that exhibit
   diploid PTEN. Additionally, a positive correlation between PTEN copy
   number value and mRNA expression was found in HCC samples
   ([124]Figure 2F).

   Taken together, these data demonstrate that PTEN expression is
   downregulated in HCC and that the deletion of PTEN copy numbers is
   likely to be one of the main mechanisms that make a contribution to the
   downregulation of PTEN among HCC patients.

Identification of differentially expressed genes (DEGs) (DEmRNAs, DElncRNAs,
and DEmiRNAs)

   According to the above analysis, the ceRNA network related to PTEN can
   be used as a potential prognostic model for HCC patients. Also, we must
   know that the meaning of the expression levels in HCC samples with
   PTEN^high and PTEN^low expression groups are opposite to those in
   cancer and paracancerous groups.

   To verify this conjecture, we first identified the DElncRNAs, DEmiRNAs,
   and DEmRNAs in HCC samples with PTEN^high and PTEN^low expression
   groups as well as in HCC and adjacent normal tissues using TCGA
   database, with p <0.05 and |log fold change [FC]| >0.5 as the lncRNA
   threshold, p <0.05 and |logFC| >0.3 as the miRNA threshold, and p <0.05
   and |logFC| >0.7/0.5 as the mRNA threshold.

   In total, 860 DElncRNAs (137 upregulated and 723 downregulated), 54
   DEmiRNAs (21 upregulated and 33 downregulated), and 1,871 DEmRNAs (195
   upregulated and 1676 downregulated) were sorted out from HCC samples
   with PTEN^high and PTEN^low expression groups.

   Furthermore, a total of 3,371 DElncRNAs (3,041 upregulated and 330
   downregulated), 420 DEmiRNAs (102 upregulated and 318 downregulated),
   and 8,294 DEmRNAs (7,059 upregulated and 1,235 downregulated) were
   identified between HCC samples and normal liver tissue samples. Volcano
   plots visually displaying the distribution of DElncRNAs, DEmiRNAs, and
   DEmRNAs were generated ([125]Figures 3A–3C and [126]S2A–S2C), and the
   heatmaps depict the expression of 15 significant variable genes in HCC
   samples with PTEN^high and PTEN^low expression, as well as in HCC and
   normal samples ([127]Figures 3D–3F and [128]S2D–S2F).

Figure 3.

   [129]Figure 3
   [130]Open in a new tab

   Volcano plots and heatmap plots of DElncRNAs, DEmiRNAs, and DEmRNAs
   between the expression of PTEN^high and PTEN^low in HCC samples

   Red represents upregulated genes and blue indicates downregulated
   genes. (A–C) The volcano plots describe (A) 860DElncRNAs (|log[2]fold
   change| > 0.5 and adjusted p value < 0.05), (B) 54DEmiRNAs (|log[2]fold
   change| > 0.3 and adjusted p value < 0.05), and (C) 1871DEmRNAs
   (|log[2]fold change| > 0.5 and adjusted p value < 0.05). (D–F) The
   horizontal axis of the heatmap indicates the samples, and the vertical
   axis of the heatmap indicates 15 significant DEGs.

Construction of the lncRNA-miRNA-mRNA triple regulatory network

   To establish the lncRNA-miRNA-mRNA triple regulatory network in HCC, we
   conducted a joint analysis in the PTEN^high and PTEN^low expression
   groups as well as in HCC and normal hepatic tissue groups. First, we
   put the remaining DElncRNAs into the TarBase database to identify
   potential miRNAs targeting lncRNAs. Four out of the predicted miRNAs
   were selected after taking the intersection with 54 DEmiRNAs. We then
   used the databases of miRDB and TargetScan to identify the downstream
   target mRNAs with reference to the four DEmiRNAs. In addition, we
   sought candidate mRNAs that were only shared by the two databases to
   enhance the veracity of the prediction. The results revealed that 373
   out of 8,294 DEmRNAs were identified. Finally, a total of 5 lncRNAs (4
   upregulated and 1 downregulated), 4 miRNAs (1 upregulated and 3
   downregulated), and 372 mRNAs (326 upregulated and 46 downregulated)
   were included to form the HCC-associated lncRNA-miRNAmRNA triple
   regulatory network by using Cytoscape software ([131]Figure 4A).

Figure 4.

   [132]Figure 4
   [133]Open in a new tab

   Construction and functional enrichment analysis of the lnRNA-miRNA-mRNA
   triple regulatory network

   The ellipses denote miRNAs, diamonds denote mRNAs, and round rectangles
   denote lncRNAs. (A) The triple regulatory network in HCC. Red indicates
   upregulated, and blue represents downregulated. (B) Thirteen hub genes
   in this network with a score of >2. (C) Functional enrichment analysis
   (GO and KEGG) of the DEmRNAs in the network

   The Cytoscape plug-in cytoHubba was utilized to determine the hub
   triple regulatory network. The results showed that three lncRNAs
   (DLEU2L, FAM99A, and ARRDC1-AS1), four miRNAs (miR-99a-5p, miR-100-5p,
   miR-9-5p, and miR-125b-5p), and six mRNA (TAOK1, HS3ST3B1, RHOQ, AGO2,
   BAZ2A, and NR6A1) were identified (see [134]Figure 4B).

   To further explore the potential functions associated with the triple
   regulatory network, functional enrichment analysis (including GO and
   KEGG) was utilized by Metascape. The results showed that the DEmRNAs
   participating in the network were particularly enriched in the “protein
   kinase activity,” “cell morphogenesis involved in differentiation,” and
   “GTPase binding” ([135]Figure 4C).

Construction and validation of the ceRNA network and selection of a model
with HCC-specific prognostic value

   To establish a crucial ceRNA of great prognostic value in HCC, we first
   analyzed the expression levels of RNAs from the hub triple regulatory
   network in HCC samples with PTEN^high and PTEN^low expression groups
   and in HCC and adjacent normal liver tissues.

   As mentioned above, the expression levels in HCC samples with PTEN^high
   and PTEN^low expression groups are opposite to those in cancer and
   paracancerous groups. Our results showed one downregulated (DLEU2L) and
   two undifferentiated (ARRDC1-AS1 and FAM99A) lncRNAs, one upregulated
   (miR-9-5p) and three downregulated (miR-99a-5p, miR-100-5p, and
   miR-125b-5p) miRNAs, three downregulated (TAOK1, HS3ST3B1, and AGO2)
   and three undifferentiated (RHOQ, BAZ2A, and NR6A1) mRNAs in HCC
   samples with PTEN^high and PTEN^low expression groups ([136]Figure 5A),
   while we found two upregulated (DLEU2L and ARRDC1-AS1) and one
   downregulated (FAM99A) lncRNAs, one upregulated (miR-9-5p) and three
   downregulated (miR-99a-5p, miR-100-5p, and miR-125b-5p) miRNAs, and
   five upregulated (TAOK1, RHOQ, AGO2, BAZ2A, and NR6A1) and one
   downregulated (HS3ST3B1) mRNAs in HCC and adjacent normal liver tissues
   ([137]Figure 5B).

Figure 5.

   [138]Figure 5
   [139]Open in a new tab

   The distribution of 13 hub-RNA expression patterns from the triple
   regulatory network in TCGA HCC dataset

   (A and B) The expression patterns of three hub-DElncRNAs, four
   hub-DEmiRNAs, and six hub-DEmRNAs in HCC samples with PTEN^high and
   PTEN^low expression groups (A), and in HCC and adjacent normal liver
   tissues (B).

   Additionally, the expression levels of these RNAs were verified in 50
   paired HCC samples from TCGA cohort. An almost consistent finding with
   the above analysis was elucidated ([140]Figure S3).

   Then, to determine whether these RNAs were associated with HCC
   prognosis, we used Kaplan-Meier analysis and a log-rank test to perform
   OS analysis of HCC patients. In total, one DElncRNA (DLEU2L), two
   DEmiRNAs (miR-99a-5p and miR-100-5p), and three DEmRNAs (TAOK1,
   HS3ST3B1, and RHOQ) were found to be related to prognosis based on p
   <0.05 ([141]Figure 6).

Figure 6.

   [142]Figure 6
   [143]Open in a new tab

   Overall survival analysis for the RNAs in the hub triple regulatory
   network

   (A–C) The high- and low-expression values of three hub-lncRNAs (A),
   four hub-miRNAs (B), and six hub-mRNAs (C) were compared by a
   Kaplan-Meier survival curve for TCGA HCC patient cohort. The horizontal
   axis indicates the overall survival time in days, and the vertical axis
   represents the survival rate.

   Moreover, considering that cellular localization of lncRNAs determined
   the underlying mechanisms, we analyzed the subcellular localization of
   the three DElncRNAs by performing the lncLocator. As shown in
   [144]Figure 7A, DLEU2L was mainly located in the cytoplasm, but the
   other two lncRNAs (FAM99A and ARRDC1-AS1) were mainly distributed in
   the cytosol. These data indicate that DLEU2L may act as a ceRNA to
   improve the expression of TAOK1 through sponging miR-99a-5p/miR-100-5p.
   Thus, a DLEU2L-miR-100-5p/miR-99a-5p-TAOK1 ceRNA network was
   constructed ([145]Figure 7B). The target sites in the DLEU2L and TAOK1
   3′ UTRs were predicted to pair with miR-99a-5p and miR-100-5p by
   TarBase and TargetScan, respectively ([146]Figure 7C).

Figure 7.

   [147]Figure 7
   [148]Open in a new tab

   Construction and correlation analysis of the ceRNA network

   (A) The cellular localization for three hub-lncRNAs (DLEU2L, FAM99A,
   and ARRDC1-AS1) was predicted using lncLocator. (B) Schematic model of
   ceRNA. Blue indicates downregulated; red indicates upregulated. (C)
   Base pairing between miR-99a-5p and miR-100-5p and the target site in
   the DLEU2L and TAOK1 3′ UTR predicted by TarBase and TargetScan,
   respectively. (D) Correlation analysis between these four predictive
   RNAs and PTEN in HCC.

   Furthermore, expression correlation analysis indicated a positive
   relationship between DLEU2L expression and TAOK1 expression
   ([149]Figure 7D). Therefore, the DLEU2L/TAOK1 axis in the ceRNA network
   was selected as a potential prognostic model for the following step
   analysis.

Clinical relevance of the DLEU2L/TAOK1 axis in HCC Patients

   To determine whether the expression levels of DLEU2L and TAOK1 were
   affected by clinical characteristics, we explored the correlation of
   DLEU2L expression and TAOK1 expression with clinical factors. Our study
   demonstrated that both DLEU2L and TAOK1 expression levels were
   positively correlated with the diameter of the tumor (DLEU2L, p =
   0.031; TAOK1, p = 0.029) and tumor-node-metastasis (TNM) stage (DLEU2L,
   p = 0.048; TAOK1, p = 0.032), as shown in [150]Table S2. However, we
   did not observe any association between the expression levels of DLEU2L
   and TAOK1 and age, sex, lymph node metastasis, distant metastasis, and
   BMI (p > 0.05; [151]Table S2)

   To further analyze the prognostic significance of clinical features, we
   conducted univariate and multivariate Cox regression analyses to
   determine OS-related features. In univariate Cox regression analysis
   models of DLEU2L and TAOK1, some prognostic factors (TNM stage, tumor
   diameter, lymph node metastasis, and distant metastasis) were closely
   relevant to OS (p < 0.05; [152]Tables S3 and [153]S4) in HCC patients
   in TCGA cohorts.

   Importantly, both DLEU2L (hazard ratio [HR] = 1.572, p = 0.011) and
   TAOK1 (HR = 1.537, p = 0.016) overexpression levels were significantly
   associated with a worse prognosis ([154]Table S4). In multivariate Cox
   regression analysis of TAOK1, tumor diameter (HR = 2.245, p = 0.000),
   distant metastasis (HR = 1.325, p = 0.003), and TAOK1^high expression
   (HR = 1.622, p = 0.007) were still relevant to OS ([155]Table S4) in
   HCC patients. However, multivariate Cox regression analysis of DLEU2L
   revealed that elevated DLEU2L expression was not correlated with a poor
   prognosis (HR = 1.369, p = 0.078; [156]Table S3). Thus, TAOK1 may
   become an independent prognostic factor for HCC patients.

Validation of TAOK1 abnormally high expression

   For the sake of better understanding the role of the DLEU2L/TAOK1 axis
   in HCC, we next analyzed TAOK1 in detail. First, based on the online
   tool Cancer Cell Line Encyclopedia (CCLE), we found that TAOK1 was high
   expression in HCC cell lines ([157]Figure S4A). Additionally, we
   performed the difference analysis of TAOK1 in HCC patients with 20
   paired HCC samples from the GEO: [158]GSE41804 HCC cohorts in the GEO
   database to verify the expression of TAOK1 further ([159]Figure S4B).
   The results revealed that the expression of TAOK1 in the HCC tissues
   was markedly higher than that of the paired nontumorous liver tissues,
   which coincided with our results from TCGA data ([160]Figures 5A and
   [161]6).

   In addition, considering that the genetic mutation might act as a
   potential mechanism to promote the abnormal overexpression of TAOK1, we
   performed the following analysis. In [162]Figure S5A, an OncoPrint plot
   shows the amplification of TAOK1 gene in TCGA HCC dataset. However, no
   significant associations were observed for TAOK1 expression and copy
   number values among HCC samples ([163]Figures S5B and S5C). Hence,
   these data demonstrated no correlation between TAOK1 abnormal
   expression and copy number alterations in HCC.

Relationship between methylation and expression of TAOK1

   To further elucidate the abnormal upregulated mechanisms of TAOK1 in
   HCC tissues, we also explored the correlation of TAOK1 expression
   levels and their methylation status by various methods. First,
   differential expression analysis of three DNA methyltransferases
   (DNMT1, DNMT3A, and DNMT3B) between TAOK1^high and TAOK1^low groups of
   HCC showed that expression of the DNMT1, DNMT3A, and DNMT3B in
   TAOK1^high group was markedly increased compared with that in TAOK1^low
   group ([164]Figure 8A). Second, the analysis of UALCAN showed that
   DNMT1 had a trend of higher methylation level in normal liver tissues
   than in HCC tissues (p = 0.101, [165]Figure 8B). Similarly, DiseaseMeth
   version 2.0 analysis demonstrated that the methylation of TAOK1 was
   significantly higher in paracancerous normal tissues compared with HCC
   tissues (p < 0.0001; [166]Figure 8C). Additionally, we found two
   methylation sites (cg16008158 and cg2570978) in the DNA sequences of
   TAOK1 that were negatively associated with their expression levels
   ([167]Figure 8D). Third, the differential methylation regions related
   to TAOK1 were presented as heatmaps ([168]Figure S6). Interestingly, a
   NEURL1B-associated methylation site cg2570978 was found in the open sea
   region and 3′ UTR region ([169]Figure S6).

Figure 8.

   [170]Figure 8
   [171]Open in a new tab

   Methylation analysis of TAOK1

   (A) Differential expression of three DNA methyltransferases (DNMT1,
   DNMT3A, and DNMT3B). (B) Methylation was evaluated using UALCAN. (C)
   Methylation was assessed using DiseaseMeth version 2.0. (D) The
   methylation site of TAOK1 DNA sequence association with gene expression
   was visualized using MEXPRESS. The expression of TAOK1 is illustrated
   by the blue line in the center of the plot. Pearson’s correlation
   coefficients and p values for methylation sites and query gene
   expression are shown on the right side.

Correlation between immune infiltration and expression of TAOK1 in HCC

   It has been reported that tumor-infiltrating lymphocytes (TILs) are
   independent predictors of sentinel lymph node (SLN) status and survival
   in certain cancers.[172]^37 To evaluate the potential relationship
   between TAOK1 expression and immune infiltration levels in HCC, we
   conducted the following analysis by using TIMER. First, the “SCNA”
   module analysis showed that several immune cell infiltration levels
   seemed to associate with altered TAOK1 gene copy numbers, including B
   cells, CD4^+ T cells, macrophages, and dendritic cells (DCs) in HCC
   ([173]Figure 9A). Then, the “Gene” module analysis showed that TAOK1
   expression has no notable correlation with tumor purity, but it is
   markedly positively correlated with infiltrating levels of CD4^+
   T cells, macrophages, neutrophils, and dendritic cells in HCC
   ([174]Figure 9B). Finally, we further evaluated the impact of immune
   infiltration on the clinical prognosis of patients with HCC.

Figure 9.

   [175]Figure 9
   [176]Open in a new tab

   Correlation analysis of TAOK1 expression and immune infiltration in HCC

   (A) Association between TAOK1 gene copy number and immune cell
   infiltration levels in HCC cohorts. (B) Correlation of TAOK1 expression
   with immune infiltration level in HCC. (C) Kaplan-Meier plots were used
   to analyze the immune infiltration and overall survival rate of HCC.
   ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001.

   The results showed that high levels of CD4^+ cells, macrophages, and
   neutrophils were associated with poor prognosis of HCC patients with a
   survival period of less than 24 months (p = 0.036, 0.007, and 0.01,
   respectively). These results indicate that the DLEU2L/TAOK1 axis may
   affect HCC and clinical prognosis by regulating the level of
   tumor-infiltrating immune cells.

   To further validate the relationship between TAOK1 and diverse
   immune-infiltrating cells, we also explored the associations between
   TAOK1 and immune marker sets of 16 immune cells.[177]^38 Interestingly,
   we found that the expression levels of five markers (BDCA-4 [NRP1],
   STAT1, STAT6, STAT3, and STAT5B) of dendritic cells, T helper (Th)1,
   Th2, Th17, and regulatory T cells (Tregs) have markedly positive
   correlations with LAYN expression in HCC ([178]Tables S5 and [179]S6).

   We also investigated the relationship between TAOK1 expression and the
   above markers of dendritic cells, Th cells (Th1, Th2 and Th17), and
   Tregs in the GEPIA database in HCC. Correlation analysis results showed
   a similar trend to TIMER ([180]Tables S5 and [181]S6). These results
   suggest that tumor-infiltrating immune cells may have an influence on
   the clinical outcome of the DLEU2L/TAOK1 axis in HCC.

TAOK1-related functional enrichment analysis in HCC

   To further investigate the possible function of TAOK1 in HCC, we
   utilized GO and KEGG pathway enrichment analysis of the top 200
   correlated genes of TAOK1 in HCC ([182]Figure S7). The enrichment term
   related to TAOK1 is the mitogen-activated protein kinase (MAPK)
   signaling pathway. Moreover, GO enrichment analyses, including
   biological process (BP), cellular component (CC), and molecular
   function (MF), related to TAOK1 were mainly enriched in “mitotic cell
   cycle phase transition,” “histone modification,” and “kinase activity.”

Discussion

   HCC is one of the leading causes of cancer mortality worldwide due to
   its poor prognosis and high aggressiveness.[183]39, [184]40, [185]41,
   [186]42 Although various therapeutic strategies have been adopted for
   HCC patients, the prognosis is still poor.[187]43, [188]44, [189]45,
   [190]46 Clarifying the molecular mechanisms and processes of the
   pathogenesis of HCC and identifying promising biomarkers are essential
   for identifying new therapeutic targets and improving patients’ outcome
   with this disease.[191]47, [192]48, [193]49 According to reports, the
   ceRNA regulatory network participates in the occurrence and development
   of many human cancers, including lung cancer,[194]^50 gastric
   cancer,[195]^51 and pancreatic cancer.[196]^52 However, to our
   knowledge, few studies have focused on a comprehensive ceRNA regulatory
   network to predict the prognosis of HCC. Accordingly, in this study, we
   attempted to establish a PTEN-related ceRNA triple network in HCC, as
   well as associate it with prognosis. Increasing evidence has
   demonstrated that PTEN acts as a tumor suppressor functionally involved
   in many types of cancers, including HCC.[197]53, [198]54, [199]55
   Interestingly, we verified similar results in this study through
   survival analysis, IHC, and copy number variation analysis.

   In this study, first, the lncRNA-miRNA-mRNA triple regulatory network
   containing 5 lncRNAs, 4 miRNAs, and 372 mRNAs was obtained through in
   silico analysis. Enrichment analysis showed that the network was mainly
   concentrated in “protein kinase activity,” “cell morphogenesis involved
   in differentiation” and “GTPase binding.” Next, a key triple regulatory
   network was obtained by hub analysis based on a score >2, including
   three lncRNAs, four miRNAs, and six mRNAs. Finally, the expression
   analysis and survival analysis were performed on this hub regulatory
   network. Moreover, because the interactions in the ceRNA network exist
   only in the cytoplasm, we also performed subcellular localization
   analysis of the three lncRNAs in the network. In summary, the
   DLEU2L-hsa-miR-100-5p/ hsa-miR-99a-5p-TAOK1 overexpressed ceRNA network
   related to the prognosis of HCC was obtained.

   By searching for these genes in PubMed, we observed that miR-100-5p,
   miR-99a-5p, and TAOK1 have been studied for their roles in cancer or in
   the relationship with cancer. Stroese et al.[200]^56 provided results
   suggesting that the miRNA-99 family is an unfavorable prognostic factor
   for patients with pancreatic ductal adenocarcinoma (PDAC), and Dhayat
   et al.[201]^57 performed bioinformatics analysis and concluded that
   miR-99a upregulation is involved in the control of PDAC cell death and
   cycle pathways. However, we found that the low expression of miR-99a-5p
   is related to the poor prognosis of breast cancer, lung adenocarcinoma,
   cervical cancer,[202]^58 and similar phenomena in HCC.[203]^59 Song
   et al.[204]^60 study showed that low expression of miR-100-5p was
   relevant to gross vascular invasion (VI) and worse outcome of HCC
   patients after hepatic resection. Additionally, another study suggested
   that miR-100 inhibits the migration and invasion of nasopharyngeal
   carcinoma cells by targeting insulin-like growth factor 1
   receptor.[205]^61 Agosto-Arroyo et al.[206]^62 found that the
   overexpression of TAOK1 in the ERBB2 subtype could be associated with
   the occurrence of invasive breast cancer. Thean et al.[207]^63 have
   concluded that TAOK1 as a target oncogene of miR-24 was coordinately
   upregulated with PIM2 and MIA-RAB4B to contribute to colorectal cancer.
   However, according to searches in PubMed, few studies have investigated
   the role of DLEU2L in cancer.

   In our study, DLEU2L and TAOK1 expression were significantly higher in
   HCC than in normal tissues, and high expression of DLEU2L and TAOK1 was
   associated with poor outcome by survival analysis. Opposite to DLEU2L
   and TAOK1, miR-100-5p and miR-99a-5p expression were downregulated, and
   the low expression level of IGF2BP1 indicated worse OS prognosis in
   HCC.

   In the current study, we found no correlation between TAOK1 abnormal
   overexpression and copy number alterations in HCC. According to
   research, DNA methylation plays an important role in regulating gene
   expression.[208]64, [209]65, [210]66, [211]67 Therefore, we referred to
   some websites to explore DNA methylation patterns that could account
   for the abnormal expression of TAOK1 in HCC. We found that TAOK1 was
   hypomethylated in HCC samples compared with adjacent normal ones, which
   is consistent with the observed upregulation of TAOK1 in HCC combined
   with DNA methyltransferase analysis (DNMT1, DNMT3A and DMNT3B). We also
   found that the TAOK1^high group clearly co-occurred with higher
   expression of DNMT1, DNMT3A, and DNMT3B, which promoted us to
   understand the upregulated expression of NEURL1B in HCC. Additionally,
   we found a negative correlation between some methylated sites and the
   prognosis of HCC patients. Finally, associations between TAOK1
   expression and genome-wide methylation were exploited. The results
   showed that more hypermethylation sites were nearer to the open sea
   regions ([212]Figure 5C). Moreover, all hypermethylated sites fell in
   UTR regions (3′ UTR and 5′ UTR), while more hypomethylation sites were
   in the TSS1500 and TSS200 regions. Importantly, the NEURL1B-associated
   methylation site cg2570978 was in the open sea region and 3′ UTR
   region. These findings suggest that the abnormal methylation may be
   related to HCC poor prognosis.

   Previous studies have reported that immune infiltration can affect the
   prognosis of the patient.[213]68, [214]69, [215]70 In this study, we
   found that TAOK1 gene copy numbers were inversely associated with the
   infiltration levels of B cells, CD8^+ T cells, CD4^+ T cells,
   macrophages, neutrophils, and dendritic cells in lung cancer by using
   TIMER. The immune cell infiltration levels changed along with the
   ACK1/TNK2 gene copy numbers. We found that several immune cell
   infiltration levels (B cells, CD8^+ T cells, CD4^+ T cells,
   macrophages, and dendritic cells) were inversely correlated with TAOK1
   gene copy numbers in HCC. We also found that the expression of TAOK1
   was highly correlated with immune infiltration of HCC. Additionally,
   several types of tumor-infiltrating immune cells were notably related
   to the prognosis of HCC patients.[216]71, [217]72, [218]73, [219]74,
   [220]75 Furthermore, significant positive correlations were revealed
   between TAOK1 expression and some immune marker sets in dendritic
   cells, Th cells (Th1, Th2, and Th17) and Tregs in HCC. These findings
   together suggest that these differences induced by the DLEU2L/TAOK1
   axis may have an impact on the changes in the tumor immune
   microenvironment and the development of HCC.

   To further explore the biological functions of TAOK1, we conducted GO
   and KEGG enrichment analyses for TAOK1. The enrichment term related to
   TAOK1 is MAPK signaling pathway. GO enrichment analysis (including BP,
   CC, and MF) related to TAOK1 mainly enriched in mitotic cell cycle
   phase transition, histone modification, and kinase activity.

   Although we have constructed the ceRNA-based DLEU2L/TAOK1 axis, which
   appears to be a potential prognostic biomarker in clinical
   applications, several limitations must also be noted. First, the
   binding affinity of lncRNAs, miRNAs, and mRNAs obtained from the
   database should be further experimentally investigated. Second, we need
   to further study the function and mechanism of the DLEU2L/TAOK1 axis in
   HCC by experiments.

   In conclusion, we established a ceRNA
   (DLEU2L-hsa-miR-100-5p/hsa-miR-99a-5p-TAOK1) overexpressed network
   related to the prognosis of HCC, which is better for knowing the
   correlation among lncRNA-miRNA-mRNA. Additionally, we identified that
   the ceRNA-based DLEU2L/TAOK1 axis can be a novel significant prognostic
   factor participating in HCC, and the prognostic model is useful for
   exploring the pathogenesis of HCC.

Materials and methods

Data preparation and processing

   [221]Figure 1 shows the workflow chart of this study. We downloaded
   primitive data (expressing profiles) and clinic information of human
   HCC from TCGA database ([222]https://portal.gdc.cancer.gov/). Available
   mRNA sequencing (mRNA-seq) data for 371 HCC samples (including 50
   paired HCC samples) and miRNA sequencing (miRNA-seq) data from 367 HCC
   samples (including 50 paired HCC samples) were gained from TCGA
   database. All raw RNA-seq data (lncRNAs, miRNAs, and mRNAs) were
   normalized as fragments per kilobase of exon model per million mapped
   fragment reads. Transformation of miRNA sequences into human mature
   miRNA names was accomplished using the starBase version 2.0 database
   ([223]http://starbase.sysu.edu.cn).

   Meanwhile, we also downloaded the gene expression profiles from GEO:
   [224]GSE41804 (including 20 HCC samples and 20 adjacent nontumor
   samples; platform, GEO: [225]GPL570), which were extracted from the GEO
   database ([226]https://www.ncbi.nlm.nih.gov/geo/), to further validate
   our results. CCLE ([227]https://portals.broadinstitute.org/ccle) and
   HPA ([228]http://www.proteinatlas.org/) searches were performed to
   verify the expression of PTEN and TAOK1 in HCC cell lines and at the
   protein level. We obtained the mutation status of PTEN and TAOK1
   through the cBioPortal for Cancer Genomics
   ([229]http://www.cbioportal.org/).

Screening of DEGs

   When performing the differential expression analysis in PTEN^high and
   PTEN^low HCC samples, we determined the DElncRNAs with thresholds of
   |logFC| >0.5 and p <0.05, the DEmiRNAs with thresholds of |logFC| >0.3
   and p <0.05, and the DEmRNAs with thresholds of |logFC| >0.5 and p
   <0.05. When performing differential expression analysis between HCC
   samples and adjacent nontumorous samples, the critical value of
   DElncRNAs was set to p <0.05 and |logFC| >0.5, the critical value of
   DEmiRNAs was set to p <0.05 and |logFC| >0.3, and the critical value of
   DEmRNAs was set to p <0.05 and |logFC| >0.7. Volcano plots of the
   DERNAs (including DElncRNAs, DEmiRNAs, and DEmRNAs) were visualized
   using GraphPad Prism 8 software (version 8.3.0). The heatmap clustering
   was drawn by TBtools software (version 1.051).

Establishment of the ceRNA network in HCC

   According to the hypothesis that lncRNA could indirectly regulate mRNA
   expression by competing with miRNA as a natural sponge in the
   cytoplasm, the ceRNA network was constructed by the following steps:
   (1) TarBase (TarBase v7.0 database, [230]http://microrna.gr/tarbase/)
   was used to predict the potential miRNAs targeted by DElncRNAs and the
   lncRNA-miRNA interaction pairs; (2) miRDB ([231]http://www.mirdb.org/)
   and TargetScan (version 7.2, [232]http://www.targetscan.org/) were used
   to forecast the target genes of the DEmiRNAs and build the miRNA-mRNA
   interaction pairs; (3) the VennDiagram package in R software was
   utilized to compare the target genes with DEmRNAs, and the target genes
   that overlapped with DEmRNAs in this study were selected for the next
   analysis; and (4) by integrating the lncRNA-miRNA pairs with miRNA-mRNA
   pairs, we built the lnRNA-miRNA-mRNA triple regulatory network.

   LNCipedia ([233]https://lncipedia.org/) was used to obtain DElncRNA
   sequences, and the lncLocator
   ([234]http://www.csbio.sjtu.edu.cn/bioinf/lncLocator/) database was
   used to identify the DElncRNA cellular localization based on its
   sequence. The Cytoscape plug-in cytoHubba was performed to identify the
   hub triple regulatory network. The generated networks were visualized
   by Cytoscape software (version 3.7.0, [235]https://www.cytoscape.org/).

Functional enrichment analysis

   For the sake of understanding the possible biological processes and
   pathways of the network, we conducted the following analysis. To begin
   with, we conducted a functional enrichment analysis of the DEmRNAs in
   the lncRNA-miRNA-mRNA triple regulatory network in Metascape
   ([236]http://metascape.org/gp/index.html). Then, the top 200
   TAOK1-associated genes in LIHC was obtained from GEPIA
   ([237]http://gepia.cancer-pku.cn/), and GO enrichment (including BP,
   CC, and MF) and KEGG pathway analyses of these genes were performed
   using Metascape and visualized by the R package ggplot2. An adjusted p
   <0.05 was determined to be statistically significant.

Survival analysis and construction of a specific prognosis model for HCC

   The survival status and time of HCC patients were gained from TCGA
   clinical dataset. We used GraphPad Prism 8 software to perform
   Kaplan-Meier analysis and a log-rank test on DElncRNAs, DEmiRNAs, and
   DEmRNAs in the ceRNA network to determine the relationship with the OS
   of HCC patients in TCGA database. Univariate Cox regression analysis
   was utilized to analyze the association between candidate genes in the
   ceRNA network and OS to determine the prognostic-related biomarkers.
   Multivariate Cox regression analysis was performed on TAOK1 to
   determine the independent prognostic factors of HCC. A log-rank p <0.05
   was determined to be statistically significant.

Methylation and expression analysis of TAOK1

   Studies showed that DNA methylation is a significant epigenetic
   mechanism, which is able to regulate gene expression by three DNA
   methyltransferases (DNMT1, DNMT3A, and DNMT3B) and influence the
   behavior of cancer cells. First, we investigated the expression level
   of three DNA methyltransferases in TAOK1^high and TAOK1^low groups by
   TCGA database. Then, we utilized the human disease methylation database
   DiseaseMeth version 2.0
   ([238]http://bio-bigdata.hrbmu.edu.cn/diseasemeth/) and UALCAN
   ([239]http://ualcan.path.uab.edu/) to assess methylation levels of
   TAOK1 between the HCC and paracancerous normal tissues. Moreover, we
   investigated the association between hub genes’ expression and their
   DNA methylation status using MEXPRESS ([240]https://mexpress.be).
   Finally, MethSurv ([241]https://biit.cs.ut.ee/methsurv/) was used to
   perform multivariate survival analysis to evaluate the scattering of
   different CpG islands.

Immune infiltrate levels and expression analysis of TAOK1

   To investigate the association of the expression of TAOK1 and
   tumor-infiltrating immune cells, we applied TIMER
   ([242]https://cistrome.shinyapps.io/timer/), which is an online tool
   for the analysis and visualization of the correlation between immune
   infiltrate levels and a number of variables across diverse cancer
   types. We explored the correlation of TAOK1 expression with the
   abundance of tumor-infiltrating immune cells (including B cells, CD4^+
   T cells, CD8^+ T cells, neutrophils, macrophages, and dendritic cells),
   the prognostic value, and TAOK1 copy numbers in HCC.

   Furthermore, we estimated the correlation of TAOK1 with the markers of
   16 tumor-infiltrating immune cells (including markers of CD8^+ T cells,
   T cells, B cells, monocytes, TAMs, M1 macrophages, M2 macrophages,
   neutrophils, natural killer [NK] cells, dendritic cells, Th1 cells, Th2
   cells, follicular Th [Tfh] cells, Th17 cells, Tregs, and exhausted
   T cells).

Statistical analysis

   The obtained data were analyzed using GraphPad Prism (version 8.0) and
   SPSS 23.0 software (SPSS, Chicago, IL, USA). The results are presented
   as the median and 95% confidence interval (CI). The Mann-Whitney test
   and an independent t test were used to calculate the difference between
   two groups of data, while a one-way ANOVA with the Kruskal-Wallis test
   and chi-square test were utilized to evaluate the difference among
   different groups. p <0.05 was considered as statistically difference.

Acknowledgments