Research Paper Volume 12, Issue 19 pp 19399—19420
Comprehensive analysis of TGF-β-induced mRNAs and ncRNAs in hepatocellular carcinoma
- 1 Hepatic Surgery Center, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
- 2 Department of Anesthesiology, Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology, Wuhan, China
Received: March 7, 2020 Accepted: July 14, 2020 Published: October 4, 2020https://doi.org/10.18632/aging.103826
How to Cite
Copyright: © 2020 Liang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Transforming growth factor β (TGF-β) is a potent inducer of epithelial-mesenchymal transition (EMT) in hepatocellular carcinoma (HCC), and plays a critical role in its tumorigenesis and progression. Accumulating evidence indicates that protein-coding mRNAs, as well as non-coding RNAs (ncRNAs), may play key roles in the tumorigenesis and progression of HCC. In this study, we first report on the differential expression of lncRNAs, mRNAs, miRNAs, and circRNAs in Huh7 cells treated with TGF-β or DMSO for 7 days. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed for significantly differentially expressed RNAs (DE RNAs). Then the competing endogenous RNA (ceRNA) network based on these DE RNAs was predicted and constructed. Among them, we identified that lncRNA SLC7A11-AS1 and hsa_circ_0006123 are involved in the EMT process induced by TGF-β and may promotes the metastasis of HCC. This knowledge may pave the way to develop novel clinical diagnostics and therapeutic approaches. Our study might open a new avenue for future investigations of the molecular mechanisms driving the EMT process induced by TGF-β in HCC.
Hepatocellular carcinoma (HCC) is one of the most common and aggressive human malignancies and the second leading cause of cancer-related deaths worldwide . Effective treatments are limited due to the poor prognosis and high-recurrence rate, which is largely due to the high rate of tumor metastasis . Therefore, improvements in HCC prevention and treatment strategies require a better understanding of the underlying pathophysiological mechanisms.
The transforming growth factors β(TGF-β) signaling pathway is a key player in tumorigenesis and cancer progression [3, 4]. TGF-β exerts a tumor-suppressive effect by inducing a cellular cytostatic program [5–7], but paradoxically, TGF-β is also known to function as a tumor promoter because of its prominent role in enhancing proliferation, migration and invasion [4, 8, 9]. It is well known that TGF-β promotes tumor invasion and metastasis in large part by inducing the epithelial-mesenchymal transition (EMT) [10, 11], which is a critical step in tumor invasion and dissemination [12, 13]. During the EMT process, cells lose some components of epithelial cell junctions and gain invasive properties . TGF-β is a potent inducer of the EMT in HCC, and there are many underlying mechanisms that should be explored to fully understand the relationship between the TGF-β signaling pathway and EMT [15, 16].
Competing endogenous RNAs (ceRNAs) are RNAs that share miRNA recognition elements (MREs), and consequently compete for miRNA binding, contributing to mutual regulation lncRNAs, circRNAs, mRNAs, and pseudogenes [17, 18]. The lncRNAs are functionally defined as transcripts of more than 200 nucleotides in length with no protein coding potential . The circRNAs are evolutionarily conserved transcripts featuring covalently linked 5‘ and 3’ ends that are derived from pre-mRNA back-splicing . Numerous recent studies have demonstrated that ceRNAs play an important role in tumorigenesis and cancer progression by affecting the expression of key oncogenes or tumor suppressor genes. For example, a recent quantitative study demonstrated that ciRs7 acts as a sink that absorbs miR-7, influencing its target genes in many human cancers [21–23]. Jiang et al. reported that long noncoding RNA 00976 promotes pancreatic cancer progression through OTUD7B by sequencing miR-137 involved in the EGFR/MAPK pathway .
In this study, we performed high-throughput RNA-sequencing (RNA-Seq) in four paired Huh7 hepatoma cell cultures continuously treated with 10 ng/ml TGF-β or DMSO for 7 days to comprehensively identify differentially expressed lncRNAs, circRNAs, miRNAs, and mRNAs. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)pathway analyses were performed for the differentially expressed mRNAs (DE mRNAs) and target mRNAs of differentially expressed lncRNAs, miRNAs, and circRNAs. The ceRNA network in the EMT process of HCC induced by TGF-β will provide new insights into the potential therapeutics or novel diagnostics approaches for HCC. Then, we confirmed that the lncRNA SLC7A11-AS1 and hsa_circ_0006123 are induced by TGF-β and may be promotes the metastasis of HCC. Individual knockdown of lncRNA SLC7A11-AS1 or hsa_circ_0006123, significantly reduced the migration and invasion ability of HCC cells. Our study might offer a new avenue for future investigations to better understand the molecular mechanisms of the EMT process induced by TGF-β in HCC.
Identification of an EMT cell model
To identify ncRNAs that are regulated by TGF-β and involved in the EMT process induced by TGF-β, we treated huh7 and HLF cells with 10ng/ml TGF-β or DMSO for 7 days. TGF-β treated huh7 cells show a spindle-shipped appearance (Figure 1A). Furthermore, the in vivo migration and invasion ability of TGF-β treated cells was significantly increased (Figure 1B). Moreover, TGF-β treatments reduced the mRNA and protein expression levels of known epithelial-associated gene E-cadherin, while the expression of mesenchymal associated genes such as ZEB1, SNAIL1, Twist and N-cadherin was upregulated (Figure 1D, 1E). As shown in the (Supplementary Figure 1), Although TGF-β treated HLF cells did not show significant morphological differences, their in vivo migration and invasion ability were significantly increased. In the case of TGF-β stimulation, E-cadherin was downregulated and ZEB1, SNAIL1 and N-cadherin was upregulated.
Figure 1. Identification of the EMT cell model based on TGF-β treatment. (A) Phase-contrast micrographs of Huh7 cells treated with 10 ng/ml TGF-β for 7 days or left untreated. Scale bar = 100 μm; (B) In vitro migration and invasion of Huh7 cells treated with 10 ng/ml TGF-β or left untreated. The number of migrated cells per field was counted after 48 h. (C) Relative mRNA levels of E-cadherin, ZEB1, SNAIL1, Twist and N-cadherin in Huh7 cells treated with 10 ng/ml TGF-β for 7 days or left untreated. (D) Relative protein levels of E-cadherin, ZEB1, SNAIL1, Twist and N-cadherin in Huh7 cells treated with 10 ng/ml TGF-β for 7 days or left untreated.
Identification of differentially expressed lncRNAs, mRNAs miRNAs and circRNAs in TGF-β-treated huh7 cells
Second-generation sequencing technology was used to detected the differentially expressed lncRNAs, mRNAs, miRNAs and circRNAs between TGF-β-treated and DMSO treated Huh7 cells. Volcano plots were used to reveal differentially expressed lncRNAs, mRNAs miRNAs and circRNAs with various p-values and fold changes (Figure 2A). Among them, we identified 158 up- and 124 downregulated lncRNAs (Supplementary Table 1), 2037 up- and 1749 downregulated mRNAs (Supplementary Table 2), 69 up- and 34 downregulation miRNAs (Supplementary Table 3), as well as 31 up- and 26 downregulation circRNAs (Supplementary Table 4). Hierarchical clustering of these DE lncRNAs, DE mRNAs, DE miRNAs and DE circRNAs illustrated an obvious difference between TGF-β treated and DMSO treated huh7 cells (Figure 2B). Among the DE coding and non-coding RNAs, we uncovered 202 novel lncRNAs and 33 novel circRNAs that have not been reported before (Supplementary Tables 5, 6). These novel DEncRNAs suggest a promising new pattern in the EMT process induced by TGF-β in HCC.
Figure 2. Expression profiles of differentially expressed RNAs. (A) Volcano plot showing differentially expressed lncRNAs, mRNAs, miRNAs and circRNAs with various p-values and fold changes. x axis: log2 ratio of RNA expression levels between treated and untreated Huh7 cells. y axis: false discovery rates (-log10 transformed) of different RNAs. Red points, upregulated RNAs; blue points, downregulated RNAs. (B) Hierarchical clustering analysis based on the significantly differentially expressed lncRNAs, mRNAs, miRNAs, and circRNAs, respectively (Fold Change > 2 and P-Value < 0.05). Red, high relative expression; blue, low relative expression; white, no change in gene expression. Color brightness reflects the degree of expression increase or decrease.
GO and KEGG enrichment analysis of differentially expressed lncRNAs, mRNAs, miRNAs and circRNAs
LncRNAs usually act on neighboring target genes, which is known as the cis role of lncRNAs [25, 26]. We searched for DE mRNAs derived from 100kb upstream and downstream of DE lncRNAs to predict putative cis target genes of DE lncRNAs, followed by analyzing the functions of these DE mRNAs to annotate the lncRNAs. The prediction results indicate that 282 lncRNAs may act on 450 target genes, forming a total of 754 lincRNA–gene connections (Supplementary Table 7). To explore the potential biological functions of these identified lncRNAs, we performed statistical enrichment of GO and KEGG terms among the target genes of the DE lncRNAs. Based on the results of GO analysis, several important GO terms were enriched, including “trigeminal nerve development”, “TAP complex”, “bone morphogenesis”, “retina layer formation”, and “regulation of oxidoreductase activity” (Figure 3A). KEGG pathway analysis showed that DE lncRNAs target genes are involved in several pathways such as, “neomycin”, “kanamycin” and “gentamicin biosynthesis”, “complement” and “coagulation cascades”, “riboflavin metabolism” and “staphylococcus aureus infection” (Supplementary Figure 2A). To predict the potential functional implications of the DE mRNAs, we performed GO and KEGG functional enrichment analyses. For DE mRNAs, the most relevant GO terms included “extracellular region”, “plasma membrane”, “type I interferon signaling pathway”, and” response to virus” (Figure 3B). The most relevant KEGG pathways were “PI3K-Akt signaling pathway”, “ECM-receptor interaction”, “small cell lung cancer” and “ovarian steroidogenesis” (Supplementary Figure 2B). To predict the potential functional implications of the DE miRNAs, we predicted target mRNAs of DE miRNAs using Miranda and TargetScan, after which GO and KEGG pathway enrichment analyses were conducted based on the target mRNAs. The result of GO enrichment analysis indicated that the DE miRNAs are closely associated with the terms “protein binding”, “cytoplasm”, “cytosol” and “transferase activity” (Figure 3C). Similarly, the data of KEGG pathway enrichment analysis indicated that the DE miRNAs were closely associated with the “p53 signaling pathway”, “pathways in cancer”, “pyrimidine metabolism” and “autophagy – animal” (Supplementary Figure 2C). Using the same approach, we also performed GO and KEGG analysis or the DE circRNAs, which suggested that the GO terms “voltage-gated potassium channel activity”, “negative regulation of viral entry into host cell”, and “positive regulation of ubiquitin-protein transferase activity” are enriched (Figure 3D), while “lysine degradation” was the only significantly enriched KEGG pathway (Supplementary Figure 2D). Many of the GO terms and KEGG pathway above have been reported correlated with EMT process, such as regulation of oxidoreductase activity, protein binding, transferase activity, positive regulation of ubiquitin-protein transferase activity, PI3K-Akt signaling pathway, p53 signaling pathway, pathways in cancer and so on [5, 27, 28].
Validation of sequencing data using real-time PCR
We randomly selected differentially expressed transcripts for validation via qRT-PCR, including 5 DE lncRNAs, 5 DE mRNAs, 5 DE miRNAs and 5 DE circRNAs. The qRT-PCR was performed to determine the expression of these lncRNAs, mRNAs, miRNAs and circRNAs in paired TGF-β-treated and DMSO-treated Huh7 cell cultures. As shown in (Figure 4A–4D), the qRT-PCR results were consistent with the RNA-seq data.
Identification of lncRNAs SLC7A11-AS1 and hsa_circ_0006123 as potentially related to the metastasis of HCC
To screen lncRNAs and circRNAs that are involved in the TGF-β-induced EMT of HCC, we specifically selected the top 10 DE lncRNAs and known DE circRNAs for qRT-PCR validation in paired TGF-β-treated and DMSO-treated Huh7 and HLF cell cultures. The results indicated that 7 of 10 DE lncRNAs, as well as 6 of 10 known DE circRNAs were differentially expressed in the paired TGF-β-treated and DMSO-treated Huh7 cell cultures (Figure 5A and 5B). Similarly, 6 of 10 DE lncRNAs and 5 of 10 known DE circRNAs were differentially expressed in the paired TGF-β-treated and DMSO-treated HLF cell cultures (Figure 5C and 5D). Among these, we selected 3 DE lncRNAs and 3 DE circRNAs with the most significant expression differentials for further validation in 81 pairs of HCC samples and adjacent non-cancerous tissues using qRT-PCR. Among these six no-coding RNAs, lncRNA SLC7A11-AS1, LINC01224, hsa_circ_0006123 and hsa_circ_0005480 were found to be upregulated in HCC, while the expression of the others was similar to that in non-cancerous tissues (Figure 6A and 6B). Next, we divided the 81 patients into two groups based on the presence or absence of vascular invasion of the tumor. Reprocessing of the qRT-PCR data revealed that lncRNA SLC7A11-AS1 and hsa_circ_0006123 were upregulated in the patients with vascular invasion compared with the patients without vascular invasion. The other 4 no-coding RNAs did not show significant differences between the two groups of patients (Figure 6C and 6D). As shown in Table 1, overexpression of lncRNA SLC7A11-AS1 was related to vascular invasion and TNM stage. Moreover, overexpression of hsa_circ_0006123 was related to vascular invasion, TNM stage and BCLC stage. Baseline information on cancer history of all the cases are shown in Supplementary Table 11. Among them, the expression of the upregulated lncRNA SLC7A11-AS1 was confirmed by comparing with the TCGA database (Figure 6E). At the same time, Kaplan-Meier analysis using the TCGA database revealed that high lncRNA SLC7A11-AS1 levels in HCC tissues are significantly correlated with shorter overall survival (OS) (P = 0.0031; log-rank test; Figure 6F) and recurrence-free survival (RFS) (P = 0.0014; log-rank test; Figure 6G). Therefore, our data indicate that lncRNA SLC7A11-AS1 and hsa_circ_0005480 may play an important role in TGF-β-mediated HCC invasion and metastasis.
Figure 5. Expression levels of the top ten dysregulated lncRNAs and known circRNAs assessed by quantitative real time-PCR (qRT-PCR). Relative expression levels of ten dysregulated lncRNAs (A–C) and known circRNAs (B–D) in Huh7 and HLF cells treated with 10 ng/ml TGF-β for 7 days or left untreated.
Figure 6. The lncRNAs SLC7A11-AS1 and hsa_circ_0006123 may be related to the metastasis of HCC. Expression levels of three selected lncRNA s (A) and three selected circRNAs (B) in 81 pairs of HCC samples and adjacent non-cancerous tissues using qRT-PCR. Expression levels of three selected lncRNAs (C) and three selected circRNAs (D) in patients with vascular invasion (Tumor-M) and patients without vascular invasion (Tumor-N). (E) Analysis of the lncRNA SLC7A11-AS1 expression level in TCGA liver cancer samples. (F) Kaplan–Meier analysis of overall survival in liver cancer patients stratified according to their SLC7A11-AS1 expression levels in the TCGA liver cancer cohort. (G) Kaplan–Meier analysis of diseases-free survival of liver cancer patients stratified according to their lncRNA SLC7A11-AS1 expression levels in the TCGA liver cancer cohort.
Table 1. Clinicopathologic characteristics of patients with hepatocellular carcinoma.
|Clinicopathological variables||SLC7A11-AS1 expression||P value||hsa_circ_0006123 expression||P value|
|Low (n=45)||High (n=36)||Low (n=38)||High (n=43)|
|Tumor size (cm)|
|II -III -IV||25||28||20||33|
Construction and analysis of ceRNA networks of lncRNA LNCSLC7A11-AS1 and hsa_circ_0006123
According to the ceRNA hypothesis, lncRNAs and circRNAs can sequester relevant miRNAs through MREs to post-transcriptionally regulate gene expression . We used our RNA-seq data to construct ceRNA networks according to the interaction mechanism of lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA complexes (Supplementary Table 8). GO analysis was performed, and several GO terms were found to be significantly enriched, including “type I interferon signaling pathway”, “cell cortex region” and “galactoside binding” (Figure 7A). KEGG pathway analysis showed that the terms “adherens junction”, “arrhythmogenic right ventricular cardiomyopathy” (ARVC) and “focal adhesion” were significantly enriched (Figure 7B). To further uncover the functions of lncRNA SLC7A11-AS1 and hsa_circ_0006123, two theoretical ceRNA networks derived from the global ceRNA network were predicted. The first network contained lncRNA SLC7A11-AS1, 4 downregulated miRNAs and 338 upregulated mRNAs (Figure 7C; Supplementary Table 9). The second network contained hsa_circ_0006123, 1 downregulated miRNA and 108 upregulated mRNAs (Figure 7D; Supplementary Table 10). GO (Supplementary Figure 3A and 3B), as well as KEGG pathway analyses (Supplementary Figure 3C and 3D) were performed to identify the potential functions of the two partial ceRNA networks. These results suggested that lncRNA SLC7A11-AS1 and hsa_circ_0006123 may function as ceRNAs that sequester miRNAs, thus implicating them in the pathogenesis of HCC.
Figure 7. LncRNA SLC7A11-AS1 and hsa_circ_0006123 ceRNA network. GO (A) and KEGG (B) pathway analysis of the global ceRNA network. The ceRNA network of lncRNA SLC7A11-AS1 (C) and hsa_circ_0006123 (D) with largest degrees. The yellow nodes represent lncRNAs, the green nodes represent circRNAs, the red nodes represent miRNAs, and the white nodes represent mRNAs.
Knockdown of lncRNA SLC7A11-AS1 and hsa_circ_0006123 suppresses the migration and invasion of HCC cells
Based on the connection between higher levels of lncRNA SLC7A11-AS1 or hsa_circ_0006123 and vascular invasion in HCC patients, we further assessed the possible roles of lncRNA SLC7A11-AS1 and hsa_circ_0006123 in modulating HCC cell migration and invasion. Individual knockdown of lncRNA SLC7A11-AS1 or hsa_circ_0006123, significantly reduced the migration and invasion ability of Huh7 and HLF cells (Figure 8A, 8B and Supplementary Figure 4A, 4B). These results clearly demonstrate that the lncRNA SLC7A11-AS1 and hsa_circ_0006123 promote the migration and invasion ability of HCC cell lines. Moreover, knockdown of lncRNA SLC7A11-AS1 in Huh7 and HLF cells reduced protein expression levels of ZEB1 and N-cadherin, while the expression of E-cadherin was upregulated (Figure 8C and Supplementary Figure 4C). Knockdown of hsa_circ_0006123 increased E-cadherin while it decreased ZEB1 and snail1 in Huh7 and HLF cells. All these data indicated that lncRNA SLC7A11-AS1 and hsa_circ_0006123 promote EMT processes of HCC cells (Figure 8D and Supplementary Figure 4D). qRT-PCR analysis showed that mRNA expression of lncRNA SLC7A11-AS1 and hsa_circ_0006123 were positively correlated with the ZEB1 and snail1 mRNA in 81 HCC patients (Figure 8E, 8F).
Figure 8. Knockdown of lncRNA SLC7A11-AS1 and hsa_circ_0006123 suppresses the migration and invasion of Huh7 cell. Migration and invasion of Huh7 following treatment with si-lncRNA SLC7A11-AS1 (A) or si-hsa_circ_0006123 (B). Relative protein levels of EMT markers in Huh7 cells following treatment with si-lncRNA LC7A11-AS1 (C) or si-hsa_circ_0006123 (D). qRT-PCR analysis showed that mRNA expression of lncRNA SLC7A11-AS1 and hsa_circ_0006123 were positively correlated with the ZEB1 (E) and snail1 (F) mRNA in 81 HCC patients
Although HCC is the third leading cause of cancer-related deaths worldwide , the underlying molecular mechanisms of its progression remain poorly characterized. The TGF-β signaling pathway has been increasingly recognized to be involved in tumorigenesis and progression of HCC . The EMT is an essential event in cancer progression as well as a trigger of migration, invasion, and metastasis of HCC cells, whereby TGF-β is the strongest and most well-characterized stimulator of EMT and metastasis in HCC [4, 30]. TGF-β expression is usually elevated in many advanced or metastatic tumors, thus promoting the development of advanced tumors . To explore the underlying molecular mechanisms of TGF-β-induced EMT in HCC, we constructed an EMT cell model by treating Huh7 cells with 10 ng/ml TGF-β or DMSO for 7 days. We artificially increase the concentration of TGF-β to simulate the environment in which TGF-β is highly expressed in advanced or metastatic tumors. Thus, the effect of TGFB on tumor cells can be further amplified to explore its mechanism of action. In addition to mRNAs, ncRNAs such as miRNAs, lncRNAs, and circRNAs were also demonstrated to play an essential role in the tumorigenesis and progression of HCC [32–34]. In fact, some pivotal TGF-β-regulated ncRNAs have been reported to be involved in tumorigenesis and cancer progression, including lncRNA-ATB, nc886, miR-29, miR-100 and miR-125b [15, 35–37]. However, these studies usually only focused on a single gene or protein in a single molecular pathway rather than on systematic or coordinated expression of other genes during HCC development. A comprehensive analysis of DE RNA profiles of HCC cells treated with TGF-β has, to our best knowledge, not been reported to date. To explore the functions and complex interactions of TGF-β-regulated ncRNAs and further explore and improve our knowledge of the TGF-β signaling pathway, we characterized the extensive transcription landscape related to the TGF-β-induced EMT process of HCC.
To screen lncRNAs and circRNAs that contribute to the TGF-β-induced EMT of HCC, we identified three most highly differentially expressed lncRNAs (lncRNA SLC7A11-AS1, LINC01224 and AL590004) and circRNAs (has_circ_0006123, has_circ_0005480 and has_circ_0004420). We found that lncRNA SLC7A11-AS1, LINC01224, hsa_circ_0006123 and hsa_circ_0005480 were upregulated in 81 pairs of HCC samples compared with adjacent non-cancerous tissues. The EMT induced by TGF-β is closely related to tumor metastasis, and further results showed that lncRNA SLC7A11-AS1 and hsa_circ_0006123 were upregulated in patients with vascular invasion compared with the patients without vascular invasion. LncRNA SLC7A11-AS1 has been reported to promote chemoresistance by blocking SCFβ-TRCP-mediated degradation of NRF2 in pancreatic cancer . These results suggested that lncRNA SLC7A11-AS1 and hsa_circ_0006123 may be involved in the EMT process induced by TGF-β, and therefore closely associated with the metastasis of HCC.
The “ceRNA hypothesis”, which describes how mRNAs, transcribed pseudogenes, and lncRNAs “talk” to each other using MREs as letters of a new language, was first proposed by Salmena et al. . Many lncRNAs and circRNAs have been reported to function as ceRNAs to sequester miRNAs and thereby further regulate the expression of mRNAs and their related proteins during carcinogenesis and tumor progression, including metastasis [15, 33]. Here, we showed for the first time that the ceRNA activity of these DE transcripts is related to the EMT induced by TGF-β in HCC cells (Supplementary Table 8). GO and KEGG pathway analyses were performed to identify the potential functions of the global ceRNA network. At the same time, a portion of the theoretical ceRNA network associated with lncRNA SLC7A11-AS1 and hsa_circ_0006123 was predicted. These results suggested that lncRNA SLC7A11-AS1 and hsa_circ_0006123 may act as ceRNAs.
Only a few ceRNAs were reported to be induced by TGF-β in HCC in earlier studies, and details of the ceRNA network induced by TGF-β in HCC remained unclear [15, 39, 40]. Here, we systematically investigated the ceRNA network based on differentially expressed lncRNAs, mRNAs, circRNAs and miRNAs for the first time. A part of these miRNAs had previously been linked to HCC, and they all target a number of lncRNAs, circRNAs and mRNAs, including the majority of the novel identified non-coding RNAs [41–44]. These findings illustrate the functional complexity of mRNAs and non-coding RNAs, and also indicate that these novel lncRNAs and circRNAs may be involved in the EMT process induced by TGF-β in HCC, meriting further investigation. Nevertheless, there are some limitations in our study. More complex animal experiment should be carried out to confirm that lncRNA SLC7A11-AS1 and hsa_circ_0006123 promote the migration and invasion ability of HCC cells in vivo. Furthermore, overexpression of the SLC7A11-AS1 and hsa_circ_0006123 lncRNAs should be performed to elucidate the underlying molecular mechanisms. It is possible that some of the differences in RNA expression found in this study are due to innate differences of cell lines. Our future work will focus on further investigating the molecular mechanisms of these ceRNAs, and exploring the role of the corresponding network in the progression and metastasis of HCC.
In conclusion, this is the first report on the differential expression of lncRNAs, mRNAs, miRNAs, and circRNAs in Huh7 cells treated with TGF-β or DMSO for 7 days. GO and KEGG pathway enrichment analyses, as well as constructing a ceRNA network allowed us to assess the roles and potential mechanisms of the differentially expressed transcripts. Among them, we identified that lncRNA SLC7A11-AS1 and hsa_circ_0006123 are likely involved in the EMT process induced by TGF-β, and may promotes the metastasis and of HCC. These findings will hopefully pave the way to develop novel clinical diagnostic and therapeutic approaches. Our study might therefore open a new avenue for future investigations of the molecular mechanisms driving the TGF-β-induced EMT of HCC.
Materials and Methods
Patients and specimens
This study included patients with human primary liver cancer who underwent liver resection at the Hepatic Surgery Center of Tongji Hospital between Sept 2012 and Oct 2017 and who provided their informed consent. The study was approved by the Medical Ethics Commission of Medical Ethics Committee of Tongji Hospital.
Cell lines and culture
Huh7 and HLF cell lines were purchased from the China Center for Type Culture Collection (Wuhan, China). They are cultured in Dulbecco's modified Eagle's medium (Invitrogen Corporation, Carlsbad, CA, USA) supplemented with 10% FBS (Life Technologies Inc., Gibco/Brl Division, Grand Island, NY, USA) under the culture conditions of 37°C and 5% CO2.
Construction of EMT cell model
According to the previous experimental results and recent detection results in our laboratory, we found that the TGF-β expression levels of Huh7 and HLF were relatively low. In addition, Huh7 show a well-differentiated epithelial morphology . Therefore, we believe that they have greater potential for TGF-β induced EMT. To establish the model of cancer cells undergoing EMT, cells were treated with 10 ng/ml TGF-β for 7 days with TGF-β replenishment every day . Compared with HLF, Huh7 showed more obvious morphological and migration changes after TGF-β treatment. Therefore, we chose Huh7 for the sequencing.
RNA library construction and sequencing
Total RNA was isolated and purified using Trizol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's procedure. The RNA amount and purity of each sample was quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA). The RNA integrity was assessed by Agilent 2100 with RIN number >7.0. Approximately 5 ug of total RNA was used to deplete ribosomal RNA according to the manuscript of the Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, USA). After removing ribosomal RNAs, the left RNAs were fragmented into small pieces using divalent cations under high temperature. Then the cleaved RNA fragments were reverse-transcribed to create the cDNA, which were next used to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase I, RNase H and dUTP.. An A-base is then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single-or dual-index adapters are ligated to the fragments, and size selection was performed with AMPureXP beads. After the heat-labile UDG enzyme treatment of the U-labeled second-stranded DNAs, The ligated products are amplified with PCR by the following conditions: initial denaturation at 95°C for 3 min; 8cycles of denaturation at 98°C for 15 sec, annealing at 60°C for 15 sec, and extension at 72°C for 30 sec; and then final extension at 72°C for 5 min. The average insert size for the final cDNA library was 300 bp (±50 bp). At last, we performed the paired-end sequencing on an IlluminaX Ten (LC Bio, China) following the vendor's recommended protocol.
Different expression analysis of lncRNAs, mRNAs, miRNAs and cicrRNAs
StringTie was used to perform expression level for lncRNAs, mRNA and miRNA by calculating FPKM. The differentially expressed mRNAs and lncRNAs were selected with log2 (fold change) >1 or log2 (fold change) <-1 and with statistical significance (p value < 0.05) by R package Ballgown .
Quantitative real-time (qRT)-PCR
Trizol regent (Invitrogen) was used to extracted total RNA, and FastQuant cDNA Synthesis Kit (TIANGEN, Beijing, China) was used to performed reverse transcription from 2 μg total RNAs. Quantitative real-time PCR was carried out with SuperReal PreMix Plus (SYBR Green, TIANGEN, Beijing, China) according to the manufacture’s instruction on a CFX Connect™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). For all qRT-PCR experiments, data analysis was performed byΔΔCt method and GAPDH was used as controls.
Migration and invasion assays
For transwell assays, cells were trypsinized and resuspended in culture medium without FBS after transfected with 100 nM gene-speicific siRNA (RiboBio, China) or non-targeting control employing Lipo3000 (Invitrogen) following the manufacturer’s instruction. For migration assay, 5–10 × 104 normal liver cells or HCC cells were added to the top 8-μm chamber without Matrigel. For invasion assays, 10–20 × 104 cells were added to the Matrigel-coated upper chamber. The lower chambers were filled with medium containing 10–20% serum. After 24–48 h incubation, the invaded cells at the bottom surface of the upper chamber membrane were fixed in 4% paraformaldehyde, stained with 1% crystal violet and counted under a microscope in a blinded manner.
Gene ontology and kyoto encyclopedia of genes and genomes analysis
Gene Ontology (GO) analysis annotates differentially expressed genes in terms of cell composition, molecular function, and biological processes. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis is an effective method for predicting the potential biological function of differentially expressed genes. Gene Ontology (GO) classification and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment were conducted using the clusterProfiler of R/Bioconductor1. In the GO and KEGG analyses, P value < 0.05 was used as the screening standard.
ceRNA network analysis
We searched the sequences of the lncRNAs, circRNAs and mRNAs to identify potential MREs. We used miRanda (http://www.microrna.org/microrna/home.do) to predict miRNA-binding seed-sequence sites, and an overlap of the same miRNA-binding sites on both circRNAs and mRNAs indicated potential circRNAmiRNA-mRNA interaction. ceRNA networks were constructed and visually displayed using Cytoscape software v.3.5.0 (San Diego, CA, USA).
All statistical analyses were performed using GraphPad Prism 7.0 and Statistical Program for Social Science version 22 (SPSS 22.0). All the statistical data on normal distribution are displayed as mean value ± standard deviation and and were compared with t-test. Data on abnormal distribution were compared with a nonparametric test. Clinical correlations were analyzed using the v2 test. All statistical tests were 2-sided, and p < 0.05 was regarded as statistical significance.
All authors searched the literature, designed the study, interpreted the findings and revised the manuscript. Junnan Liang, Tongtong Liu, Yu Wang, and Jingyuan Wen carried out data management and statistical analysis and drafted the manuscript. Jingyu Liao, Ning Cai performed the experiments. Zhao Huang, Weiqi Xu, Ganxun Li helped with cohort identification and data management. Zeyang Ding and Bixiang Zhang performed project administration.
We gratefully acknowledge all the participants of the study. We would like to thank the Ms. Lianjun Mu for the helpful support in writing this paper.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
This work was supported by the State Key Project on Infectious Diseases of China (2018ZX10723204-003), the National Nature Science Foundation of China (Nos. 81874065, 81500565, 81874149, 81572427, and 81401997), the Hepato-Biliary-Pancreatic Malignant Tumor Investigation Fund of Chen Xiao-ping Foundation for the Development of Science and Technology of Hubei Province (CXPJJH11800001-2018356).
- 1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018; 68:394–424. https://doi.org/10.3322/caac.21492 [PubMed]
- 2. El-Serag HB. Hepatocellular carcinoma. N Engl J Med. 2011; 365:1118–27. https://doi.org/10.1056/NEJMra1001683 [PubMed]
- 3. Majumdar A, Curley SA, Wu X, Brown P, Hwang JP, Shetty K, Yao ZX, He AR, Li S, Katz L, Farci P, Mishra L. Hepatic stem cells and transforming growth factor β in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. 2012; 9:530–38. https://doi.org/10.1038/nrgastro.2012.114 [PubMed]
- 4. Massagué J. TGFbeta in cancer. Cell. 2008; 134:215–30. https://doi.org/10.1016/j.cell.2008.07.001 [PubMed]
- 5. Derynck R, Akhurst RJ, Balmain A. TGF-beta signaling in tumor suppression and cancer progression. Nat Genet. 2001; 29:117–29. https://doi.org/10.1038/ng1001-117 [PubMed]
Jr, Bascom CC, Sipes NJ, Graves-Deal R, Weissman BE, Moses HL. Selective inhibition of growth-related gene expression in murine keratinocytes by transforming growth factor beta. Mol Cell Biol. 1988; 8:3088–93. https://doi.org/10.1128/mcb.8.8.3088 [PubMed]
- 7. Laiho M, DeCaprio JA, Ludlow JW, Livingston DM, Massagué J. Growth inhibition by TGF-beta linked to suppression of retinoblastoma protein phosphorylation. Cell. 1990; 62:175–85. https://doi.org/10.1016/0092-8674(90)90251-9 [PubMed]
- 8. Fabregat I, Fernando J, Mainez J, Sancho P. TGF-beta signaling in cancer treatment. Curr Pharm Des. 2014; 20:2934–47. https://doi.org/10.2174/13816128113199990591 [PubMed]
- 9. Batlle E, Massagué J. Transforming growth factor-β signaling in immunity and cancer. Immunity. 2019; 50:924–40. https://doi.org/10.1016/j.immuni.2019.03.024 [PubMed]
- 10. Heldin CH, Vanlandewijck M, Moustakas A. Regulation of EMT by TGFβ in cancer. FEBS Lett. 2012; 586:1959–70. https://doi.org/10.1016/j.febslet.2012.02.037 [PubMed]
- 11. Oft M, Peli J, Rudaz C, Schwarz H, Beug H, Reichmann E. TGF-beta1 and ha-ras collaborate in modulating the phenotypic plasticity and invasiveness of epithelial tumor cells. Genes Dev. 1996; 10:2462–77. https://doi.org/10.1101/gad.10.19.2462 [PubMed]
- 12. Thiery JP, Acloque H, Huang RY, Nieto MA. Epithelial-mesenchymal transitions in development and disease. Cell. 2009; 139:871–90. https://doi.org/10.1016/j.cell.2009.11.007 [PubMed]
- 13. Chaffer CL, Weinberg RA. A perspective on cancer cell metastasis. Science. 2011; 331:1559–64. https://doi.org/10.1126/science.1203543 [PubMed]
- 14. Ocaña OH, Córcoles R, Fabra A, Moreno-Bueno G, Acloque H, Vega S, Barrallo-Gimeno A, Cano A, Nieto MA. Metastatic colonization requires the repression of the epithelial-mesenchymal transition inducer Prrx1. Cancer Cell. 2012; 22:709–24. https://doi.org/10.1016/j.ccr.2012.10.012 [PubMed]
- 15. Yuan JH, Yang F, Wang F, Ma JZ, Guo YJ, Tao QF, Liu F, Pan W, Wang TT, Zhou CC, Wang SB, Wang YZ, Yang Y, et al. A long noncoding RNA activated by TGF-β promotes the invasion-metastasis cascade in hepatocellular carcinoma. Cancer Cell. 2014; 25:666–81. https://doi.org/10.1016/j.ccr.2014.03.010 [PubMed]
- 16. Giannelli G, Bergamini C, Fransvea E, Sgarra C, Antonaci S. Laminin-5 with transforming growth factor-beta1 induces epithelial to mesenchymal transition in hepatocellular carcinoma. Gastroenterology. 2005; 129:1375–83. https://doi.org/10.1053/j.gastro.2005.09.055 [PubMed]
- 17. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013; 495:384–88. https://doi.org/10.1038/nature11993 [PubMed]
- 18. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell. 2011; 146:353–58. https://doi.org/10.1016/j.cell.2011.07.014 [PubMed]
- 19. Iyer MK, Niknafs YS, Malik R, Singhal U, Sahu A, Hosono Y, Barrette TR, Prensner JR, Evans JR, Zhao S, Poliakov A, Cao X, Dhanasekaran SM, et al. The landscape of long noncoding RNAs in the human transcriptome. Nat Genet. 2015; 47:199–208. https://doi.org/10.1038/ng.3192 [PubMed]
- 20. Chen LL. The biogenesis and emerging roles of circular RNAs. Nat Rev Mol Cell Biol. 2016; 17:205–11. https://doi.org/10.1038/nrm.2015.32 [PubMed]
- 21. Weng W, Wei Q, Toden S, Yoshida K, Nagasaka T, Fujiwara T, Cai S, Qin H, Ma Y, Goel A. Circular RNA ciRS-7-a promising prognostic biomarker and a potential therapeutic target in colorectal cancer. Clin Cancer Res. 2017; 23:3918–28. https://doi.org/10.1158/1078-0432.CCR-16-2541 [PubMed]
- 22. Hansen TB, Kjems J, Damgaard CK. Circular RNA and miR-7 in cancer. Cancer Res. 2013; 73:5609–12. https://doi.org/10.1158/0008-5472.CAN-13-1568 [PubMed]
- 23. Li RC, Ke S, Meng FK, Lu J, Zou XJ, He ZG, Wang WF, Fang MH. CiRS-7 promotes growth and metastasis of esophageal squamous cell carcinoma via regulation of miR-7/HOXB13. Cell Death Dis. 2018; 9:838. https://doi.org/10.1038/s41419-018-0852-y [PubMed]
- 24. Lei S, He Z, Chen T, Guo X, Zeng Z, Shen Y, Jiang J. Long noncoding RNA 00976 promotes pancreatic cancer progression through OTUD7B by sponging miR-137 involving EGFR/MAPK pathway. J Exp Clin Cancer Res. 2019; 38:470. https://doi.org/10.1186/s13046-019-1388-4 [PubMed]
- 25. Zhao F, Qu Y, Liu J, Liu H, Zhang L, Feng Y, Wang H, Gan J, Lu R, Mu D. Microarray profiling and co-expression network analysis of LncRNAs and mRNAs in neonatal rats following hypoxic-ischemic brain damage. Sci Rep. 2015; 5:13850. https://doi.org/10.1038/srep13850 [PubMed]
- 26. Dong R, Du J, Wang L, Wang J, Ding G, Wang S, Fan Z. Comparison of long noncoding RNA and mRNA expression profiles in mesenchymal stem cells derived from human periodontal ligament and bone marrow. Biomed Res Int. 2014; 2014:317853. https://doi.org/10.1155/2014/317853 [PubMed]
- 27. Shaul YD, Freinkman E, Comb WC, Cantor JR, Tam WL, Thiru P, Kim D, Kanarek N, Pacold ME, Chen WW, Bierie B, Possemato R, Reinhardt F, et al. Dihydropyrimidine accumulation is required for the epithelial-mesenchymal transition. Cell. 2014; 158:1094–109. https://doi.org/10.1016/j.cell.2014.07.032 [PubMed]
- 28. Su J, Morgani SM, David CJ, Wang Q, Er EE, Huang YH, Basnet H, Zou Y, Shu W, Soni RK, Hendrickson RC, Hadjantonakis AK, Massagué J. TGF-β orchestrates fibrogenic and developmental EMTs via the RAS effector RREB1. Nature. 2020; 577:566–71. https://doi.org/10.1038/s41586-019-1897-5 [PubMed]
- 29. Llovet JM, Zucman-Rossi J, Pikarsky E, Sangro B, Schwartz M, Sherman M, Gores G. Hepatocellular carcinoma. Nat Rev Dis Primers. 2016; 2:16018. https://doi.org/10.1038/nrdp.2016.18 [PubMed]
- 30. Giannelli G, Koudelkova P, Dituri F, Mikulits W. Role of epithelial to mesenchymal transition in hepatocellular carcinoma. J Hepatol. 2016; 65:798–808. https://doi.org/10.1016/j.jhep.2016.05.007 [PubMed]
- 31. Lippitz BE. Cytokine patterns in patients with cancer: a systematic review. Lancet Oncol. 2013; 14:e218–28. https://doi.org/10.1016/S1470-2045(12)70582-X [PubMed]
- 32. Yang F, Zhang L, Huo XS, Yuan JH, Xu D, Yuan SX, Zhu N, Zhou WP, Yang GS, Wang YZ, Shang JL, Gao CF, Zhang FR, et al. Long noncoding RNA high expression in hepatocellular carcinoma facilitates tumor growth through enhancer of zeste homolog 2 in humans. Hepatology. 2011; 54:1679–89. https://doi.org/10.1002/hep.24563 [PubMed]
- 33. Wei Y, Chen X, Liang C, Ling Y, Yang X, Ye X, Zhang H, Yang P, Cui X, Ren Y, Xin X, Li H, Wang R, et al. A noncoding regulatory RNAs network driven by circ-CDYL acts specifically in the early stages hepatocellular carcinoma. Hepatology. 2020; 71:130–47. https://doi.org/10.1002/hep.30795 [PubMed]
- 34. Roy S, Hooiveld GJ, Seehawer M, Caruso S, Heinzmann F, Schneider AT, Frank AK, Cardenas DV, Sonntag R, Luedde M, Trautwein C, Stein I, Pikarsky E, et al. microRNA 193a-5p regulates levels of nucleolar- and spindle-associated protein 1 to suppress hepatocarcinogenesis. Gastroenterology. 2018; 155:1951–66.e26. https://doi.org/10.1053/j.gastro.2018.08.032 [PubMed]
- 35. Ahn JH, Lee HS, Lee JS, Lee YS, Park JL, Kim SY, Hwang JA, Kunkeaw N, Jung SY, Kim TJ, Lee KS, Jeon SH, Lee I, et al. Nc886 is induced by TGF-β and suppresses the microRNA pathway in ovarian cancer. Nat Commun. 2018; 9:1166. https://doi.org/10.1038/s41467-018-03556-7 [PubMed]
- 36. Lyu G, Guan Y, Zhang C, Zong L, Sun L, Huang X, Huang L, Zhang L, Tian XL, Zhou Z, Tao W. TGF-β signaling alters H4K20me3 status via miR-29 and contributes to cellular senescence and cardiac aging. Nat Commun. 2018; 9:2560. https://doi.org/10.1038/s41467-018-04994-z [PubMed]
- 37. Ottaviani S, Stebbing J, Frampton AE, Zagorac S, Krell J, de Giorgio A, Trabulo SM, Nguyen VT, Magnani L, Feng H, Giovannetti E, Funel N, Gress TM, et al. TGF-β induces miR-100 and miR-125b but blocks let-7a through LIN28B controlling PDAC progression. Nat Commun. 2018; 9:1845. https://doi.org/10.1038/s41467-018-03962-x [PubMed]
- 38. Yang Q, Li K, Huang X, Zhao C, Mei Y, Li X, Jiao L, Yang H. lncRNA SLC7A11-AS1 promotes chemoresistance by blocking SCFβ-TRCP-mediated degradation of NRF2 in pancreatic cancer. Mol Ther Nucleic Acids. 2020; 19:974–85. https://doi.org/10.1016/j.omtn.2019.11.035 [PubMed]
- 39. Wang B, Hsu SH, Majumder S, Kutay H, Huang W, Jacob ST, Ghoshal K. TGFbeta-mediated upregulation of hepatic miR-181b promotes hepatocarcinogenesis by targeting TIMP3. Oncogene. 2010; 29:1787–97. https://doi.org/10.1038/onc.2009.468 [PubMed]
- 40. Yang P, Li QJ, Feng Y, Zhang Y, Markowitz GJ, Ning S, Deng Y, Zhao J, Jiang S, Yuan Y, Wang HY, Cheng SQ, Xie D, Wang XF. TGF-β-miR-34a-CCL22 signaling-induced treg cell recruitment promotes venous metastases of HBV-positive hepatocellular carcinoma. Cancer Cell. 2012; 22:291–303. https://doi.org/10.1016/j.ccr.2012.07.023 [PubMed]
- 41. Huang W, Huang F, Lei Z, Luo H. LncRNA SNHG11 promotes proliferation, migration, apoptosis, and autophagy by regulating hsa-miR-184/AGO2 in HCC. Onco Targets Ther. 2020; 13:413–21. https://doi.org/10.2147/OTT.S237161 [PubMed]
- 42. Hu X, Feng Y, Sun L, Qu L, Sun C. Roles of microRNA-330 and its target gene ING4 in the development of aggressive phenotype in hepatocellular carcinoma cells. Dig Dis Sci. 2017; 62:715–22. https://doi.org/10.1007/s10620-016-4429-2 [PubMed]
- 43. Tu J, Zhao Z, Xu M, Lu X, Chang L, Ji J. NEAT1 upregulates TGF-β1 to induce hepatocellular carcinoma progression by sponging hsa-mir-139-5p. J Cell Physiol. 2018; 233:8578–87. https://doi.org/10.1002/jcp.26524 [PubMed]
- 44. Wong CC, Wong CM, Tung EK, Au SL, Lee JM, Poon RT, Man K, Ng IO. The microRNA miR-139 suppresses metastasis and progression of hepatocellular carcinoma by down-regulating rho-kinase 2. Gastroenterology. 2011; 140:322–31. https://doi.org/10.1053/j.gastro.2010.10.006 [PubMed]
- 45. Ding ZY, Jin GN, Wang W, Chen WX, Wu YH, Ai X, Chen L, Zhang WG, Liang HF, Laurence A, Zhang MZ, Datta PK, Zhang B, Chen XP. Reduced expression of transcriptional intermediary factor 1 gamma promotes metastasis and indicates poor prognosis of hepatocellular carcinoma. Hepatology. 2014; 60:1620–36. https://doi.org/10.1002/hep.27273 [PubMed]
- 46. Lin X, Chai G, Wu Y, Li J, Chen F, Liu J, Luo G, Tauler J, Du J, Lin S, He C, Wang H. RNA m6A methylation regulates the epithelial mesenchymal transition of cancer cells and translation of snail. Nat Commun. 2019; 10:2065. https://doi.org/10.1038/s41467-019-09865-9 [PubMed]
- 47. Frazee AC, Pertea G, Jaffe AE, Langmead B, Salzberg SL, Leek JT. Ballgown bridges the gap between transcriptome assembly and expression analysis. Nat Biotechnol. 2015; 33:243–46. https://doi.org/10.1038/nbt.3172 [PubMed]