Abstract

   Hepatoblastoma (HB) is one of the most common hepatic malignancies in
   the pediatric population. HB are composed of a variety of tumors, which
   derived from different origins and had varying clinical outcomes.
   However, the unclear underlying mechanisms of HB limited exploring
   novel biomarkers and effective therapeutic targets. We searched
   microarray datasets on Gene Expression Omnibus (GEO) database and
   selected [36]GSE75271 and [37]GSE75283 datasets for comprehensive
   analysis. Weighted gene correlation network analysis (WGCNA) was
   employed to identify genes which were associated with tumor malignant
   phenotypes, including HB subtypes, Cairo classification and tumor
   stage. Coexpression analysis of identified genes was also performed and
   lncRNA-miRNA-mRNA network was finally conducted. Our results showed
   that a total of 22 lncRNAs, 13 miRNAs and 66 mRNAs were identified to
   be associated with tumor malignant phenotypes. Mechanistically, these
   molecules might promote the malignant phenotypes via regulating
   metabolic pathways. Among of them, 6 miRNAs (hsa-miR-106b,
   hsa-miR-130b, hsa-miR-19a, hsa-miR-19b, hsa-miR-20a and hsa-miR-301a),
   8 lncRNAs ([38]NR_102317, [39]XR_245338, [40]XR_428373, [41]XR_924945,
   [42]XR_929728, [43]XR_931611, [44]XR_935074 and [45]XR_946696), and 6
   mRNAs (EGFR, GAREM, INSIG1, KRT81, SAR1B and SDC1) were selected to
   conduct a lncRNA-miRNA-mRNA network. Taken together, our findings
   provide evidence for exploring molecular mechanisms of HB. Those
   identified malignant phenotype-associated molecules might be potential
   biomarkers and anti-cancer therapeutic targets in future.

   Keywords: hepatoblastoma, lncRNA, miRNA, mRNA, malignant phenotype

INTRODUCTION

   Hepatoblastoma (HB) is the most common malignant liver tumor in
   children, accounting for approximately 50% of pediatric hepatic-related
   cancers [[46]1]. The incidence of HB in children with age <15 is about
   one per million (1/1,000,000), and nearly 20% of those patients already
   have a synchronous metastasis at the first diagnosis. Despite recent
   advances in treatment, such as surgical resection, adjuvant
   chemotherapy, and liver transplantation, the prognosis in advanced HB
   stages still remains poor [[47]2, [48]3]. It has been accepted that HB
   are composed of a variety of tumors deriving from different immature
   liver precursors, including hepatocytes, biliary, and other epithelial
   or mesenchymal cells, which caused significant tumor heterogeneity
   [[49]4, [50]5]. For decades, researchers have observed the varying
   clinical outcomes in patients with different histological subgroups
   [[51]5–[52]9]. In addition, other parameters, such as tumor stage,
   distant metastasis, multifocality, patient age and birth weight, have
   also been reported to be associated with prognosis [[53]10, [54]11].
   However, the underlying mechanisms remains unknown. A thorough
   understanding of molecular mechanisms regarding tumor progression is
   essential for exploring effective therapeutic targets against HB.

   Studies have revealed that only a small proportion (1%–2%) of the
   genome encodes proteins, and the majority of the mammalian genome
   encodes plenty of non-coding RNAs [[55]12–[56]14]. As an important
   member of the non-coding RNAs, microRNAs (miRNAs), referred to 18-25
   nucleotides, have been clearly demonstrated to regulate a variety of
   cellular processes via recognizing the 3’-untranslated regions of
   specific mRNAs and suppressing the expression of target genes [[57]15].
   Long non-coding RNAs (lncRNAs), another sort of non-coding RNAs with
   more than 200 nucleotides, are found to play pivotal roles in chromatin
   remodeling, transcription regulation and post-transcriptional mRNA
   processing [[58]16–[59]19]. Recently, increasing studies reported that
   lncRNAs could also function as competing endogenous RNAs (ceRNAs) by
   competitively binding to miRNAs through their miRNA response elements
   (MRE) [[60]20]. These ceRNAs usually share MRE with other coding
   transcripts and therefore act as sponges for that cluster of miRNAs,
   protecting the targeted mRNA transcripts from degradation [[61]21].
   Although the functional mechanism that lncRNAs act as ceRNAs leading to
   miRNA deregulation has been observed in multiple malignancies, the
   lncRNA-miRNA-mRNA network in HB is far from being fully investigated.

   In the present study, we performed a comprehensive analysis based on
   mRNA, lncRNA and miRNA expression profiling data derived from 50 HB
   patients and 5 controls, and identified differentially expressed mRNA,
   lncRNA and miRNA in HB. Next, we employed weighted gene correlation
   network analysis (WGCNA) and coexpression analysis to select the
   modules which are associated with tumor malignant phenotypes. Based on
   above results, we finally conducted several potential lncRNA-miRNA-mRNA
   networks in HB.

RESULTS

Identification of the malignant phenotype-associated mRNAs and lncRNAs in HB

   [62]GSE75271 was selected for identify differentially expressed mRNA
   and lncRNAs in HB patients [[63]22]. As shown in Figure [64]1A, Then,
   we performed a probe level analysis of [65]GSE75271 by linear models
   for microarray data (LIMMA). Those probes with P-value less than 0.05
   and foldchange more than 2.0 were identified as differentially
   expressed probes. One hundred ninety-four probes were primary selected,
   which were all down-regulated in HB patients (Figure [66]1B). Next, 194
   probes were annotated by Affymetrix microarray annotation files
   according to previous reported method [[67]23]. Finally, we identified
   a total of 61 lncRNAs and 133 mRNAs from 194 differentially expressed
   probes (Figure [68]1C).

Figure 1. Identification of differentially expressed mRNA and lncRNAs from
the HB dataset.

   [69]Figure 1
   [70]Open in a new tab

   (A) Normalization of dataset [71]GSE75271. (B) A total of 194
   downregulated probes were picked up from HB patients. (C) Within the
   194 downregulated probes, 133 mRNAs and 61 lncRNAs were identified.

   To further explore the malignant phenotype-associated mRNAs and
   lncRNAs, we performed a weighted gene correlation network analysis
   (WGCNA) (Figure [72]2A) and divided 194 differentially expressed probes
   into eight module eigengenes (ME), including MEgrey, MEred, MEblue,
   MEgreen, MEbrown, MEmagenta and MEblack (Figure [73]2B). Next, we
   employed module-trait relationship analysis to conduct the association
   between eight MEs and 13 clinical traits (including race, sex, age, HB
   subtypes, Cairo classification, tumor stage, CTNNB1, NFE2L2 and TERT
   mutation status, histological types, prognosis, tumor recurrent status
   and survival time). Interestingly, we found that one module (MEblue)
   was negatively associated with HB subtypes, Cairo classification and
   tumor stage, with P-values of 1x10^-7, 2x10^-8 and 1x10^-4,
   respectively (Figure [74]2C). The blue module contained a total of 66
   mRNAs and 22 lncRNAs, which were selected as malignant
   phenotype-associated mRNAs and lncRNAs for the further analysis (Table
   [75]1).

Figure 2. The WGCNA analysis of the malignant phenotype-associated mRNAs and
lncRNAs.

   [76]Figure 2
   [77]Open in a new tab

   (A) The cluster dendrogram of differentially expressed mRNAs and
   lncRNAs derived from the dataset [78]GSE75271. (B) The cluster
   dendrogram of module eigengenes. (C) The module-trait relationship
   analysis between the 8 odules and clinical characteristics. ME: module
   eigengenes.

Table 1. Identification of HB-associated mRNAs, miRNAs and lncRNAs.

   HB-associated mRNAs HB-associated
   miRNAs HB-associated
   lncRNAs
   AC005523.3 GPRC5C REPS2 [79]NR_045387 hsa-miR-106b
   ADHFE1 GSTA1 RET [80]NR_024548 hsa-miR-130a
   AKR7A3 GSTZ1 RIPK4 [81]NR_027005 hsa-miR-130b
   ALDH6A1 HAO2 RNF144B [82]NR_040515 hsa-miR-17
   AMBP INSIG1 RPS6KA4 [83]NR_102317 hsa-miR-18a
   APBB1IP KMO SAR1B [84]NR_102357 hsa-miR-19a
   APOA1 KNG1 SDC1 [85]XR_171896 hsa-miR-19b
   CBS KRT81 SEBOX [86]XR_245338 hsa-miR-20a
   CD53 LCAT SLC22A7 [87]XR_427361 hsa-miR-25
   CDC37L1 LGR4 SLC7A2 [88]XR_428373 hsa-miR-301a
   CRAT LPA SMARCA4 [89]XR_922928 hsa-miR-301b
   CTAG2 MGAT5 STS [90]XR_923052 hsa-miR-449a
   CTAGE15 MIR6778 SYDE1 [91]XR_924107 hsa-miR-451
   CYP1A2 MUT TTC39C [92]XR_924945
   CYP2C18 NEUROG3 UGT1A1 [93]XR_924990
   CYP4F2 NFKBIZ VNN3 [94]XR_929728
   DHTKD1 OAF WFDC3 [95]XR_931611
   EGFR PCCB ZFYVE28 [96]XR_931899
   ETS2 PIGC ZG16 [97]XR_933428
   GAREM PTCRA ZNF584 [98]XR_935074
   GLYAT RARG ZSCAN5A [99]XR_946696
   GPLD1 RDH16 ZYG11A [100]XR_949848
   [101]Open in a new tab

Identification of the malignant phenotype-associated miRNAs in HB

   We also identified tumor-associated miRNAs using WGCNA method based on
   [102]GSE75283, which was miRNA profiling dataset derived from the same
   HB patients as [103]GSE75271 (Figure [104]3A). Similarly, we divided
   887 miRNAs into ten MEs, including MEgreen, MEpink, MEblack, MEtan,
   MEgreenyellow, MEmidnoghtblue, MEpurple, MEsalmon, MEblue and MEmagenta
   (Figure [105]3B). Interestingly, the module-trait relationship analysis
   revealed that MEmidnoghtblue was positively associated with HB
   subtypes, Cairo classification and tumor stage, with P-values of
   2x10^-4, 3x10^-4 and 2x10^-4, respectively (Figure [106]3C). The
   midnoghtblue module contained 13 miRNAs, which were identified as
   malignant phenotype-associated miRNAs (Table [107]1).

Figure 3. The WGCNA analysis of the malignant phenotype-associated miRNAs.

   [108]Figure 3
   [109]Open in a new tab

   (A) The cluster dendrogram of differentially expressed miRNAs derived
   from the dataset [110]GSE75283. (B) The cluster dendrogram of module
   eigengenes. (C) The module-trait relationship analysis between the 10
   modules and clinical characteristics. ME: module eigengenes.

Coexpression analysis of malignant phenotype-associated molecules

   To explore the possible relationship between malignant
   phenotype-associated molecules, including mRNAs, lncRNAs and miRNAs, we
   performed coexpression analysis. As shown in Figure [111]4, we found
   that the lncRNAs were positively associated with nearly almost mRNAs.
   Interestingly, we also found that identified miRNAs were negatively
   associated with both mRNAs and lncRNAs. These results suggested that
   these mRNAs, lncRNAs and miRNAs might form lncRNA-miRNA-mRNA network,
   which affected the progression of HB patients. In order to predict
   biological function of above identified lncRNA-miRNA-mRNA network, we
   performed Gene ontology (GO) and Kyoto Encyclopedia of Genes and
   Genomes (KEGG) pathway enrichment analysis. Biological process analysis
   suggested these molecules were mainly involved in metabolic biological
   processes, such as lipoprotein metabolic process and
   oxidation-reduction process (Figure [112]5A). Cellular components
   analysis showed that they were located at endoplasmic reticulum (ER)
   membrane, the latter of which played a critical role in liver
   metabolism (Figure [113]5B). Besides, molecular function analysis also
   showed these molecules mostly belonged to heme-binding and
   oxygen-binding protein families (Figure [114]5C). KEGG pathway analysis
   further revealed that these malignant phenotype-associated molecules
   might be involved in regulating the metabolic related pathways to
   influence carcinogenesis and tumor progression (Figure [115]5D).

Figure 4. Coexpression analysis of malignant phenotype-associated molecules.

   [116]Figure 4
   [117]Open in a new tab

   The heatmap was presented by correlation coefficient between each pair
   of malignant phenotype-associated molecules. Range of colors (red to
   blue) shows the correlation coefficient (high to low).

Figure 5. Function annotation of malignant phenotype-associated molecules.

   [118]Figure 5
   [119]Open in a new tab

   Gene ontology analysis of biological processes (A), cellular components
   (B) and molecular function (C) of those molecules. (D) The KEGG pathway
   enrichment analysis.

Construction of lncRNA-miRNA-mRNA network in HB

   We further conducted lncRNA-miRNA-mRNA network in HB based on the
   coexpression analysis of malignant phenotype-associated molecules. The
   construction of lncRNA-miRNA-mRNA network included three steps: (a)
   candidate malignant phenotype-associated lncRNAs, miRNAs and mRNAs were
   selected according to WGCNA results; (b) the correlation coefficients
   between lncRNAs, miRNAs and mRNAs were calculated based on their
   corresponding profiling data; (c) targets of miRNAs were predicted from
   Targetscan ([120]http://www.targetscan.org/). Finally, a total of 6
   miRNAs (hsa-miR-106b, hsa-miR-130b, hsa-miR-19a, hsa-miR-19b,
   hsa-miR-20a and hsa-miR-301a), 8 lncRNAs ([121]NR_102317,
   [122]XR_245338, [123]XR_428373, [124]XR_924945, [125]XR_929728,
   [126]XR_931611, [127]XR_935074 and [128]XR_946696), and 6 mRNAs (EGFR,
   GAREM, INSIG1, KRT81, SAR1B and SDC1) were selected to conduct
   lncRNA-miRNA-mRNA network. Heatmap was conducted to present the
   correlation coeffecient between molecules (Figure [129]6A), and
   Cytoscape were used for visualization of lncRNA-miRNA-mRNA network
   (Figure [130]6B).

Figure 6. Construction of lncRNA-miRNA-mRNA network in HB.

   Figure 6
   [131]Open in a new tab

   (A) Selected 6 mRNAs, 8 lncRNAs and 6 miRNAs for network. (B) Predicted
   lncRNA-miRNA-mRNA network in HB.

DISCUSSION

   HB is a common malignant hepatic tumor in children and the prognosis
   varies among different categories. To date, there are limited effective
   methods to treat HB. Hence, it’s essential to explore underlying
   molecular mechanism of HB. In the present study, we for the first time
   identified the malignant phenotype-associated lncRNAs, mRNAs and
   miRNAs, which were associated with HB subtypes, Cairo Classification
   and tumor stage. Next, we conducted the lncRNA-miRNA-mRNA network based
   on their expression patterns to provide a molecular explanation for HB
   patients.

   With the advancement of whole-genomic sequencing technologies, lncRNAs
   have attracted increasing attention [[132]24]. Amount studies have
   reported that aberrant lncRNA expressions could serve as novel
   molecular biomarkers in the cancer diagnosis and prognosis prediction
   [[133]25, [134]26]. LncRNAs could function as oncogenes, regulating
   alternation of key signaling pathways, and promoting tumor growth and
   metastasis [[135]27, [136]28]. Recently, the miRNA sponge role of
   lncRNAs and lncRNA-miRNA-mRNA network have been widely accepted
   [[137]20]. For instance, lncRNAs SPRY4-IT1 could sponge miR-101-3p and
   consequently increase EZH2 expression to promote cellular proliferation
   and metastasis of bladder cancer cells [[138]29]. LncRNA-UCC was found
   to promote colorectal cancer progression by sponging miR-143 [[139]30],
   and lncRNA-CCAT1 promotes hepatocellular carcinoma progression by
   sponging let-7 [[140]31]. However, whether the lncRNA-miRNA-mRNA
   network plays an important role in HB patients have not been fully
   explored. A recent study reported that TUG1/miR-34a-5p/VEGFA network
   was involved in regulating hypervascularity and hepatoblastoma
   progression [[141]32]. In the present study, we, for the first time,
   explored malignant phenotype-associated lncRNAs, miRNAs and mRNAs in HB
   by using WGCNA methods, respectively. A total of 22 lncRNAs, 13 miRNAs
   and 66 mRNAs were identified, which were strongly associated with HB
   subtypes, Cairo classification and tumor stage. These findings may
   provide us potential biomarkers and/or anti-cancer targets in future.

   One hallmark of cancer is the metabolic reprogramming in cancer cells,
   including glycometabolism, lipid metabolism and amino acid metabolism
   [[142]33–[143]35]. In HB cell lines HepG2, Hep3B and HuH-6, the
   metabolic reprogramming has been defined as one important hallmark and
   contributes to tumor progression [[144]36]. Several driver genes of HB,
   such as β-catenin, YAP and c-Myc [[145]37, [146]38], could also affect
   intracellular energy metabolism. Moreover, it was also be found that
   that molecular targeting of mitochondrial metabolism holds promise as a
   novel and effective therapeutic approach for HB [[147]39]. In the
   present study, we also performed pathway enrichment analysis of above
   identified molecules. Our results revealed that most molecules
   associated with HB classification and tumor stage belonged to metabolic
   pathway, indicating a close interaction between metabolic pathways and
   HB malignant phenotypes.

   In this study, we also conducted the lncRNA-miRNA-mRNA network
   consisting of 6 miRNAs, 8 lncRNAs and 6 mRNAs, which was associated
   with HB classification and tumor stage. In this conducted network, we
   identified two genes in epidermal growth factor (EGF) signaling, EGFR
   and SDC1 [[148]40], acting as effector genes. Recently, a whole
   transcriptome analysis based on HB patients suggested that aberrant EGF
   signaling was associated with HB classification [[149]41]. The loss of
   EGFR signaling members were shown to be more present in less
   differentiated embryonal and undifferentiated small cells subtypes of
   HB. Therefore, our findings might also provide evidence for explaining
   the underlying molecular variation of different HB subtypes.

   In conclusion, we identified malignant phenotype-associated lncRNAs,
   miRNAs and mRNAs by employing WGCNA method. Mechanistically, these
   molecules might promote the malignant phenotypes of HB via regulating
   metabolic pathways. Moreover, we also conducted a lncRNAs-miRNAs-mRNAs
   network based on above findings. To our best knowledge, this is the
   first comprehensive study investigating lncRNA-miRNA-mRNA network in
   HB. If validated, our findings might provide evidence for exploring
   anti-cancer target in HB patients.

MATERIALS AND METHODS

Microarray data

   Gene and miRNA profiling data of HB patients were searched from the
   Gene Expression Omnibus (GEO, [150]http://www.ncbi.nlm.nih.gov/geo),
   which is a public database containing freely available profiling
   datasets. We finally selected two datasets for bioinformatic analysis,
   including [151]GSE75271 (mRNA profiling data based on Affymetrix Human
   Genome U133 Plus 2.0 Array platform) and [152]GSE75283 (miRNA profiling
   data based on Agilent-029297 Human miRNA Microarray v15 platform),
   which were derived from the same 50 HB patients and 5 controls
   [[153]22]. All raw data were downloaded from GEO database.

Data normalization and probe annotation

   The raw data of [154]GSE75271 was stored as probe-level CEL files, and
   was quantile normalized using Robust Multi-array Average (RMA) method.
   After data normalization, expression level of each probe in
   [155]GSE75271 was obtained. The probe sequences with corresponding
   Affymetrix probe IDs were downloaded from the Affymetrix website
   ([156]http://www.affymetrix.com). The annotation of lncRNA transcripts
   were then performed according to previous reported method [[157]23]. A
   total of 8240 lncRNA transcripts were generated with RefSeq transcript
   IDs. The probe ID-centric gene expression profile was also generated
   according to Affymetrix annotation files.

Weighted gene correlation network analysis

   WGCNA is an algorithm for constructing a co-expression network, defined
   by the similarity of gene co-expression [[158]42]. In data processing,
   the genome-wide gene expression data was initially filtrated with
   measuring the consistency of gene expression profiles by Pearson
   correlation, then we utilized the power adjacent function to Pearson
   correlation matrix to transform data into weighted gene co-expression
   networks. Network module represents a cluster of closely interconnected
   genes. Finally, the adjacency matrix, a measurement of topology
   similarity, is converted into the topological overlap matrix (TOM), and
   modules are detected by cluster analysis [[159]43].

Pathway enrichment analysis

   GO analysis is freely available for users in the annotation and
   biological properties of genes, gene products and sequences [[160]44].
   KEGG is a knowledge base for systematic analysis of gene functions in
   terms of the networks of genes and molecules [[161]45]. We used GO
   analysis and KEGG analysis to identify the function of the aberrantly
   expressed candidate lncRNAs-associated genes. The P-value of each
   enriched pathway was assigned with -log10 transformation. P<0.05
   (-log10(P-value)=1.30) was considered as statistically significant.

Statistical analysis

   All data were analyzed by R software 3.4.1
   ([162]https://www.r-project.org/). For the pair of molecule i and j,
   Pearson correlation coefficient was computed. P < 0.05 was considered
   statistically significant.

Abbreviations

   ceRNAs
          competing endogenous RNAs

   ER
          endoplasmic reticulum

   EGF
          epidermal growth factor

   GEO
          Gene Expression Omnibus

   GO
          Gene ontology

   HB
          Hepatoblastoma

   KEGG
          Kyoto Encyclopedia of Genes and Genomes

   LIMMA
          linear models for microarray data

   lncRNA
          long non-coding RNAs

   miRNA
          microRNAs

   MRE
          miRNA response elements

   ME
          module eigengenes

   RMA
          Robust Multi-array Average

   TOM
          topological overlap matrix

   WGCNA
          weighted gene correlation network analysis

Footnotes

   Author contributions

   Kai Qu and Shunbin Dong: Designed the research and analyzed the data;
   Sida Liu and Fujing Xie: Collected and analyzed data; Sinan Liu:
   Constructed figures; Xiaohong Xiang and Kai Qu: Drafted and revised the
   manuscript.

   CONFLICTS OF INTEREST

   The authors declare no conflicts of interest.

   FUNDING

   This work was supported by National Science Foundation of China (No.
   81201549); Natural Science Basic Research Plan in Shaanxi Province of
   China (No. 2017JM8039); Fundamental Research Fund for the Central
   Universities (No. 2016qngz05); the Clinical Research Award of the First
   Affiliated Hospital of Xi'an Jiaotong University, China (No.
   XJTU1AF-CRF-2015-011).

REFERENCES