In silico analysis of lncRNA-miRNA-mRNA signatures related to Sorafenib effectiveness in liver cancer cells

Abstract

BACKGROUND

Hepatocellular carcinoma (HCC) is the most common subtype of primary liver cancer with varied incidence and epidemiology worldwide. Sorafenib is still a recommended treatment for a large proportion of patients with advanced HCC. Different patterns of treatment responsiveness have been identified in differentiated hepatoblastoma HepG2 cells and metastatic HCC SNU449 cells.

AIM

To define the long non-codingRNA-microRNA-mRNA (lncRNA-miRNA-mRNA) predicted signatures related to selected hallmarks of cancer (apoptosis, autophagy, cell stress, cell dedifferentiation and invasiveness) in RNAseq studies using Sorafenib-treated HepG2 and SNU449 cells. Various available software analyses allowed us to establish the lncRNA-miRNA-mRNA regulatory axes following treatment in HepG2 and SNU449 cells.

METHODS

HepG2 and SNU449 cells were treated with Sorafenib (10 μmol/L) for 24 hours. Total RNA, including small and long RNA, was extracted with a commercial miRNeasy kit. RNAseq was carried out for the identification of changes in lncRNA-miRNA-mRNA regulatory axes.

RESULTS

MALAT, THAP9-AS1 and SNGH17 appeared to coordinately regulate miR-374b-3p and miR-769-5p that led to upregulation of SMAD7, TIRARP, TFAP4 and FAXDC2 in HepG2 cells. SNHG12, EPB41 L4A-AS1, LINC01578, SNHG12 and GAS5 interacted with let-7b-3p, miR-195-5p and VEGFA in SNU449 cells. The axes MALAT1/hsa-mir-374b-3p/SMAD7 and MALAT1/hsa-mir-769-5p/TFAP4 were of high relevance for Sorafenib response in HepG2 cells, whereas PVT1/hsa-miR-195-5p/VEGFA was responsible for the differential response of SNU449 cells to Sorafenib treatment.

CONCLUSION

Critical lncRNAs acting as sponges of miRNA were identified that regulated mRNA expression, whose proteins mainly increased the antitumor effectiveness of the treatment (SMAD7, TIRARP, TFAP4, FAXDC2 and ADRB2). However, the broad regulatory axis leading to increased VEGFA expression may be related to the side effect of Sorafenib in SNU449 cells.

Key Words: Cell culture; Hepatocellular carcinoma; Non-coding RNA; RNAseq; Sorafenib

Core Tip: In the current study, we investigated the differential lncRNA-miRNA-mRNA regulatory axes in HepG2 and SNU449 under Sorafenib treatment. The study identified lncRNA-miRNA regulatory axes leading to increased expression mRNA with positive (SMAD7, TIRARP, TFAP4, FAXDC2 and ADRB2) and negative (VEGFA) therapeutic effects in HepG2 and SNU449 under Sorafenib treatment.

INTRODUCTION

Hepatocellular carcinoma (HCC) is the most common type of primary liver cancer, accounting for 75%-86% of cases[1]. The impact of risk factors on HCC has a variable geographic distribution, including hepatitis C virus and hepatitis B virus infection, alcohol, aflatoxin B1, metabolic-associated fatty liver disease (MAFLD), tobacco and congenital diseases. Furthermore, only one third of HCC patients are diagnosed at the initial stages susceptible to receiving curative treatment, the remaining patients are at intermediate/advanced stage with poor prognosis and a 5-year survival rate below 10%[2]. The Barcelona Clinic Liver Cancer classification considers the number and size of tumor nodules, liver function, general condition of the patient and the presence of vascular invasiveness and metastases[3]. According to this classification, HCC stage C includes patients with preserved liver function but with vascular invasion and/or the presence of extrahepatic nodules. The landscape of treatments for advanced stage HCC has been changing due to the approval of Atezolizumab-Bevacizumab and Durvalumab-Tremelimumab as first-line therapies, extending the median survival up to 19.2 months[3]. However, patients without preserved liver function, a high-risk of bleeding, vascular disorders, and arterial hypertension, as well as severe autoimmune disorders and prior transplantation are not recommended to receive immune checkpoint and antiangiogenic treatments[3]. Furthermore, subgroup analysis of several phase III clinical trials suggests that Sorafenib tends to be more effective than immune checkpoint inhibitor-based therapy in patients with non-viral etiologies (alcohol, MAFLD or unknown)[4].

Sorafenib is an orally administered multityrosine kinase inhibitor approved by the FDA for the treatment of HCC in 2008[5,6]. It targets several proliferation receptors such as vascular endothelial growth factor (VEGFR)-1, VEGFR-2, VEGFR-3, platelet-derived growth factor receptor (PDGFR)-β, c-KIT, FLT-3 and RET. Furthermore, it also inhibits downstream kinases with serine/threonine phosphorylating activity[7-11]. Sorafenib induces early endoplasmic reticulum stress with activation of the JNK/AMPK/autophagy-dependent pathway that leads to a switch towards apoptotic cell death in HepG2 cells[7]. Given that Sorafenib exerts a greater proapoptotic activity in HepG2 cells than in SNU449 cells[12], a link between treatment resistance processes and induction of the epithelial-mesenchymal transition (EMT) process has also been demonstrated. In this regard, SNU449, HLF and HLE liver cancer cell lines expressing mesenchymal markers (CD44, Vimentin and Snail) are refractory to Sorafenib treatment compared to HepG2, Hep3B and PLC/PRF/5, which express epithelial markers (E-cadherin and CK-18)[13]. Sorafenib resistance might also involve overexpression of ATP-binding box (ABC) transporters that export drugs and reduce the efficacy of treatment[14].

The human genome is made up of 70% intergenic regions, while genes account for only 30%. Among non-coding RNAs (ncRNAs), microRNAs (miRNAs), long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) have been related to the initiation, progression and metastasis of HCC[15]. In broad terms, miRNAs are a type of ncRNAs, 22 nucleotides in length, which bind to the complementary sequence of the mRNA and induce its degradation mediated by the RNA-Induced Silencing Complex. Although several miRNA signatures have been associated with the diagnosis of HCC, few models have assessed the risk associated with treatment effectiveness[16,17]. We have recently identified miR-200c-3p, miR-222-5p, and miR-512-3p as prognostic miRNA markers in patients under Sorafenib treatment[18]. On the other hand, lncRNAs and circRNAs are more than 200 nucleotides long; however, lncRNAs are linear, while circRNAs are ring-shaped. Both can be transcribed from exons, introns, intergenic regions or 5′/3′ non-translational regions and folded into complicated second structures, facilitating their interactions with DNA, RNA, and proteins. lncRNAs have been implicated in the regulation of transcriptional activity, and they are specifically expressed in cells in response to numerous stimuli and serve as molecular signals[19]. They can also act as scaffolds, central platforms on which the relevant molecular components will be assembled[20].

Given that Sorafenib exerts a greater proapoptotic activity in HepG2 cells than in SNU449 cells[12], the present study aimed to determine a treatment-associated lncRNA profile directly correlated to differential signatures of lncRNAs-miRNA-mRNAs related to apoptosis, autophagy, endoplasmic reticulum stress, differentiation, and metastasis in RNA-seq data obtained from Sorafenib-treated HepG2 and SNU449 cells.

MATERIALS AND METHODS

Cell culture

The study involved the differentiated hepatoblastoma cell line HepG2 (HB-8065™, ATCC-LGC Standards, S.L.U., Barcelona, Spain) and the metastatic HCC SNU449 cell line (CRL-2234 HB-8065™, ATCC-LGC Standards). Cells were cultured in MEM containing Earle's salts with L-glutamine and 10% FBS (F7524, lot BCBX9154, Merck-Sigma, St. Louis, MO, United States), sodium pyruvate (1 mmol/L) (Ref. 11360070, Thermo Fisher Scientific, Waltham, MA, United States), non-essential amino acids (Ref. 11140035, Thermo Fisher Scientific), penicillin-streptomycin solution (100 U/mL-100 μg/mL) (Ref. 15640055, Thermo Fisher Scientific) at 37 °C, 21% O₂ and 5% CO₂. Cells were plated at 100000 cells/cm², and after 24 hours of stabilization Sorafenib (10 μM) (FS10808, Carbosynth, United Kingdom) was added. Cells were harvested at 6 and 24 hours after Sorafenib treatment. Sorafenib dose, treatment times and response were performed as optimized in the study by Rodríguez-Hernández et al[7].

RNA-seq

Total RNA, including small and long RNA, was extracted with the miRNeasy kit (Ref. 217004, Qiagen, Hilden, Germany). Lysis was performed with Qiazol and the RNA fraction was bound to the RNeasy mini spin column, washed, and eluted in RNase-free water. The DNA was removed by DNase I digestion. RNA was then quantified using the NanoDrop™ One/OneC microvolume UV-Vis UV-Vis spectrophotometer (Thermo Fisher Scientific). RNA was quantified using the RNA Qubit™ HS assay (Q32852, Thermo Fisher Scientific), and the quality was assessed with the Bioanalyzer^® 2100 Eukaryote Total RNA Nano Chip (Agilent Technologies, Santa Clara, CA, United States). All samples had an RNA Integrity Number ≥ 8.8. Total RNA libraries were prepared with the Illumina Stranded Total RNA Prep Ligation with Ribo-Zero™ Plus kit (Ref. 20040529, Illumina Inc., San Diego, CA, United States). Briefly, rRNA was depleted and the total transcriptome was fragmented. cDNA was synthesized, adapters and dual indices were joined and then the library was amplified. The libraries were then quantified with the Qubit™ DNA HS assay and the quality was analyzed with the highly sensitive Bioanalyzer^® 2100 DNA chip (Agilent Technologies). The size of the total RNA libraries was 260-280 bp. They were pooled at 2.4 nM, denatured and diluted to 1.2 nM. PhiX was used as an internal quality control. Paired-end sequencing (75 cycles) was performed in S1 flow cells on a NovaSEq 6000 instrument (Illumina. Inc.).

Analysis of RNA-seq data

Long RNA-seq data analysis: Primary analysis of RNA sequencing data, including quality testing, trimming and demultiplexing, was conducted using the Illumina DRAGEN FASTQ Generation v.3.8.4 app of BaseSpace. Next, alignment to the human genome version GRCh38 was performed with the RNA-Seq Alignment tool v.2.02, using the STAR aligner. Differential expression analysis (DEA) was carried out with the DRAGEN Differential Expression app v.3.9.0 by running the DESeq2 algorithm[21]. BioMart data mining tool (Ensembl)[22] was used to determine gene type. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways annotation and enrichment of protein coding genes was performed with Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8. Gene Set Enrichment Analysis (GSEA) was carried out with the Hallmark and GO Biological Process gene sets[23,24].

Small RNA-seq data analysis: Primary analysis of small RNA sequencing data, including quality testing, trimming and demultiplexing, was conducted using the Illumina DRAGEN FASTQ Generation v.3.8.4 app of BaseSpace. Next, the adaptor sequence was eliminated, and data were filtered by size of 23 bp, taking into account a minimum size insert of 15 bp. Also, readings were cut first and last 4 nucleotides, following the manufacturer’s instructions. Mapping and alignment to the human genome version GRCh38 using bowtie aligner was performed in Galaxy RNA[25] environment with the miRDeep Mapper tool[26]. Quantitation of reads mapping to known miRBase hairpin miRNA precursor and mature miRNA sequences was done with the miRDeep Quantifier tool[26]. Finally, DEA was performed with R Studio using the DSeq2 package.

Heatmaps were constructed using R Studio with ComplexHeatmap, viridis and dplyr packages. GO terms and KEGG pathways enrichment analysis was performed with the DAVID online tool[27] using all genes with counts > 0 as background. Biological connections among selected targets were estimated using Search Tool for the Retrieval of Interacting Genes/Proteins (STRING)[28]. STRING analysis was performed with medium confidence (interaction score 0.4). Validation of proposed interacting networks based on tumor and non-tumor liver tissue expression data was performed using GEPIA using the Correlation tool (Pearson method)[29].

In silico establishment of lncRNA-miRNA-mRNA regulatory axis

lncRNA as a sponge of miRNAs: mRNAs of relevance in Sorafenib response were selected among GSEA GO Biological Process terms related to apoptosis, autophagy, endoplasmic reticulum stress, differentiation, and metastasis (including development and angiogenesis related terms). Moderately and highly expressed mRNAs in relevant significantly enriched GO terms in GSEA with log₂FC > 1 or log₂FC < -1 and P-adj < 0.05 were used as candidates of miRNA regulation. miRNAs that could target these candidates were predicted with the miRDB database, with a cut-off of Target Score ≥ 80. From this miRNA list, we selected those that were down- or up-regulated accordingly in our study. Then, miRNA-lncRNA interactions were predicted with DIANA-LncBase v3 (DIANA tools)[30,31]. lncRNAs significantly up-regulated in our study were used for further analysis.

lncRNA-mRNA interactions to regulate stability: mRNAs previously selected to be regulated by miRNAs were used to test possible interactions with regulatory lncRNAs in our study as described above. LncRRI search web server was used to test interactions with an energy threshold of -12Kcal/mol[32], using specific transcripts.

lncRNA-protein interactions: Finally, to identify lncRNA interactions with proteins whose mRNA is regulated in our study we used the RNAct tool[33], using both predicted and experimentally confirmed targets.

RESULTS

Sorafenib exerts differential anti-tumor properties in different cell lines. Whereas HepG2 cells are considered to be sensitive to Sorafenib, the poorly-differentiated SNU449 cell line showed a lower pattern of treatment response. Therefore, we investigated RNA profiles to unravel potential candidates of Sorafenib response in vitro. First, to identify pathways related to differential treatment response in HepG2 and SNU449 cells, we selected differentially expressed (DE) genes with a base mean > 30, log₂FC > 0.5 (upregulated) or < -0.5 (downregulated) and P-adjusted (P-adj) < 0.05 for Venn analysis. The study of long RNA species identified 8345 DE genes, among which 243 RNAs were common to all conditions studied, indicating a common mechanism of Sorafenib in both cell lines (Figure 1A). Most of these genes were mRNAs coding for proteins, followed by lncRNAs (Figure 1B). Hierarchical clustering allowed for the grouping of expression profiles into 4 clusters (Figure 1C). Cluster 1 corresponded with genes that were up-regulated in response to Sorafenib, and this effect was more prominent in HepG2 cells compared to SNU449 cells. Cluster 2 was composed of genes up-regulated in SNU449 cells, whereas cluster 3 was representative of genes with higher expression in HepG2 cells. Cluster 4 involved down-regulated genes after Sorafenib treatment. As in the case of cluster 1, the effect of Sorafenib was more prominent in HepG2 cells. These results confirmed the differential effectiveness of Sorafenib in these cells lines. GO term and KEGG pathway analysis was performed with genes of each cluster (Figure 1D-G). Cluster 1 was focused on the cell nucleus, controlling transcription and RNA processing (Figure 1D). Genes with increased expression in SNU449 cells (cluster 2) were involved in the regulation of focal adhesions, extracellular matrix and pathways involving MAPK and serine-threonine protein kinase activity (Figure 1E). On the other hand, cluster 3 that contained genes with higher expression in HepG2 cells were related to lipid, carbon and amino acid metabolism, mitochondria, and EMT (Figure 1F). Cluster 4, which was downregulated in response to Sorafenib, was involved in DNA repair and epigenetic chromatin remodeling processes (Figure 1G).

Figure 1 Expression profile of long RNA species and enrichment analysis. A: Venn diagram of differentially expressed (DE) genes in HepG2 and SNU449 cells treated with Sorafenib for 6 and 24 hours; B: Pie chart showing identification of RNA species; C: Heatmap of the 8345 DE genes grouped into 4 row clusters; D: Top 4 Gene Ontology terms and Kyoto Encyclopedia of Genes and Genome pathways in cluster 1; E: Cluster 2; F: Cluster 3; G: Cluster 4 identified in the heatmap.

Next, we sequenced small RNA species to identify a miRNA pattern related to Sorafenib response in each cell line. DE genes with a base mean > 15, log₂FC > 0.5 (upregulated) or < -0.5 (downregulated) and P-value < 0.05 were selected. In this case, we found 110 DE miRNAs in HepG2 and SNU449 cells (Figure 2A). miRNAs could also be clustered according to expression profiles (Figure 2B). Cluster 1 contained miRNAs with higher expression in SNU449 cells. Cluster 2 was up-regulated in response to Sorafenib. Cluster 3 was more heterogeneous and included miRNAs that were slightly up-regulated in HepG2 cells. Cluster 4 included miRNAs with higher expression in HepG2 cells.

Figure 2 Expression profile of small RNAs. A: Venn diagram of differentially expressed (DE) miRNAs in HepG2 and SNU449 cells treated with Sorafenib for 6 and 24 hours; B: Heatmap of the 110 DE miRNAs grouped into 4 row clusters.

Next, we investigated whether the differences in Sorafenib response in HepG2 and SNU449 cells were due to a differential RNA integration network. Therefore, we determined regulatory loops among mRNAs, lncRNAs and miRNAs. In this sense, we explored the role of lncRNAs as sponges of miRNAs that regulate mRNA expression, as well as their role in controlling mRNAs and lncRNA-protein interactions. We performed a targeted search of mRNAs according to previously described Sorafenib anti-tumor properties. We selected mRNAs from significantly enriched GSEA terms related to apoptosis, autophagy, endoplasmic reticulum stress, differentiation, and metastasis, including differentiation and angiogenesis (Figure 3A). In HepG2 cells, Sorafenib induced the deregulation of terms related to cell death induction (endoplasmic reticulum stress, apoptosis, autophagy) more prominently at 24 hours of treatment (34/118, 28.8%), compared to 6 hours (13/79, 16.4%) (Figure 3B, Supplementary Figure 1A, Supplementary Table 1). Surprisingly, Sorafenib had a greater impact on the processes selected at 6 hours than at longer treatment times in SNU449 cells. In particular, the effect on cell death was enriched at 6 hours (18/56, 32.1%), compared to 24 hours (3/28, 10.7%), when terms related to migration and invasion were dominant (Figure 3C, Supplementary Figure 1B, Supplementary Table 1). A total of 441 mRNAs related to Sorafenib properties were deregulated in this study, whose expression is shown in Figure 3D.

Figure 3 Identification of mRNA candidates. A: Diagram representing the study design; B: Bubble plots of enriched GSEA terms in HepG2 cells; C: SNU449 cells treated with Sorafenib 24 hours. Bubble size is representative of gene count and color is representative of the normalized enrichment score; D: Heatmap of selected mRNAs in HepG2 and SNU449 cells.

An in silico prediction of miRNAs targeting selected mRNAs in Figure 3D was carried out to search for sources of post-transcriptional mRNA regulation (Figure 4A and B; Supplementary Figure 2). From these miRNAs, we selected for further analysis those that were DE in our study with a P-value < 0.05 and log₂FC < -1 or > 1 (Figure 4C). A total of 49 putative miRNA-mRNA interactions were identified, in which 23 unique miRNAs participated. Different miRNA profiles were observed in each condition, being hsa-miR-19b-3p and hsa-miR-29b-3p specific to HepG2 cells treated with Sorafenib at 6 hours, compared to 24 hours, a timepoint that involved the specific regulation of hsa-miR-27b-5p, hsa-miR-193b-3p, hsa-miR-194-3p, hsa-miR-374a-3p, hsa-miR-374b-3p, hsa-miR-769-5p, hsa-miR-3187-3p, hsa-miR-4488, hsa-miR-4521 and hsa-miR-7974. HRK, JMY, TNFAIP3, RORA and ZNF385B mRNAs were specifically regulated at 6 hours of Sorafenib treatment, whereas MBNL1, CREB3 L3, ARHGAP24, HOXA1, SH3GLB1, OXR1, GOLGA4, OPTN, SIAH1, LAMC2, PLK2, ATOH8, SMAD7, TIRARP, TFAP4, FAXDC2, and MELTF were specific to 24 hours. ARRDC3 was the only common mRNA in HepG2 cells at 6 and 24 hours of treatment. Although both cell lines shared miR-374 family (hsa-miR-374c-5p in SNU499 cells), hsa-miR-4488 was the only common miRNA shared by treatments at 24 hours, constituting hsa-let-7b-3p, hsa-miR-29b-1-5p, hsa-miR-32-3p, hsa-miR-34a-5p, hsa-miR-195-5p, hsa-miR-12136, hsa-miR-4492 and hsa-miR-5701 a specific signature of SNU449 cells (24 hours). Coherently, no mRNAs were common in both cell lines at 24 hours. HEY1, CEBPB, VEGFA, ADM, YOD1, TNFAIP3, KLF4, ADRB2, HSPA1B, SYNGR3, RFLNB, ACVR2B, NR4A2, CSPG5, AREG, NRARP, MAF, IGFBP5, and PDGFB were specifically related to SNU449 cells (24 hours). As in the case of HepG2 cells, Sorafenib at 6 hours had a low effect on miRNA-mRNA based interactions (hsa-miR-15a-5p, hsa-miR-12135) (CHAC1,FOXJ1) (Table 1).

Figure 4 Identification of miRNA candidates. A: Gene regulatory network showing miRNA-mRNA target predictions in HepG2 cells; B: SNU449 cells treated with Sorafenib for 24 hours; C: Heatmap of selected miRNAs.

Table 1 miRNA-mRNA interactions observed in HepG2 and SNU449 cells treated with Sorafenib for 6 and 24 hours.

Condition	hsa-miR	log2 (Fold-change)	P value	P-adj	mRNA target	Transcript ID	log2 (Fold-change)	P-adj
HepG2 6 hours	hsa-miR-29b-3p	-1.374	3.23E-02	0.9961	HRK	ENSG00000135116.9	1.967	8.926E-08
HepG2 6 hours	hsa-miR-29b-3p	-1.374	3.23E-02	0.9961	JMY	ENSG00000152409.9	1.105	2.253E-13
HepG2 6 hours	hsa-miR-29b-3p	-1.374	3.23E-02	0.9961	TNFAIP3	ENSG00000118503.15	1.186	4.131E-07
HepG2 6 hours	hsa-miR-19b-3p	-1.277	4.29E-02	0.9961	RORA	ENSG00000069667.16	1.048	1.506E-03
HepG2 6 hours	hsa-miR-19b-3p	-1.277	4.29E-02	0.9961	TNFAIP3	ENSG00000118503.15	1.186	4.131E-07
HepG2 6 hours	hsa-miR-19b-3p	-1.277	4.29E-02	0.9961	ZNF385B	ENSG00000144331.20	1.159	1.233E-02
HepG2 6 hours	hsa-miR-19b-3p	-1.277	4.29E-02	0.9961	ARRDC3	ENSG00000113369.9	1.697	4.636E-30
HepG2 24 hours	hsa-miR-4521	-4.496	1.09E-04	0.1484	MBNL1	ENSG00000152601.17	1.242	2.599E-16
HepG2 24 hours	hsa-miR-7974	-3.432	1.97E-42	0.0000	CREB3 L3	ENSG00000060566.14	1.062	4.001E-06
HepG2 24 hours	hsa-miR-194-3p	-1.463	2.73E-08	0.0000	ARHGAP24	ENSG00000138639.18	1.281	1.186E-08
HepG2 24 hours	hsa-miR-194-3p	-1.463	4.47E-05	0.0020	HOXA1	ENSG00000105991.9	1.192	2.484E-06
HepG2 24 hours	hsa-miR-374a-3p	-1.327	4.47E-05	0.0020	SH3GLB1	ENSG00000097033.14	1.296	1.800E-13
HepG2 24 hours	hsa-miR-374a-3p	-1.327	8.95E-03	0.1859	OXR1	ENSG00000164830.18	1.050	2.089E-08
HepG2 24 hours	hsa-miR-374a-3p	-1.327	8.95E-03	0.1859	GOLGA4	ENSG00000144674.16	1.230	8.419E-03
HepG2 24 hours	hsa-miR-374a-3p	-1.327	8.95E-03	0.1859	ARRDC3	ENSG00000113369.9	2.849	2.203E-41
HepG2 24 hours	hsa-miR-27b-5p	-1.236	8.95E-03	0.1859	OPTN	ENSG00000123240.17	1.531	6.346E-20
HepG2 24 hours	hsa-miR-193b-3p	-1.193	3.74E-04	0.0143	SIAH1	ENSG00000196470.12	1.581	1.144E-38
HepG2 24 hours	hsa-miR-193b-3p	-1.193	2.29E-02	0.3864	LAMC2	ENSG00000058085.15	1.327	2.152E-03
HepG2 24 hours	hsa-miR-3187-3p	-1.169	2.29E-02	0.3864	PLK2	ENSG00000145632.15	1.634	3.009E-08
HepG2 24 hours	hsa-miR-3187-3p	-1.169	3.59E-02	0.4538	ATOH8	ENSG00000168874.13	1.691	6.502E-08
HepG2 24 hours	hsa-miR-374b-3p	-1.064	3.59E-02	0.4538	SMAD7	ENSG00000101665.9	1.999	1.669E-42
HepG2 24 hours	hsa-miR-374b-3p	-1.064	7.47E-03	0.1693	TIRARP	ENSG00000163659.13	1.070	1.555E-18
HepG2 24 hours	hsa-miR-769-5p	-1.058	7.47E-03	0.1693	TFAP4	ENSG00000090447.12	1.010	1.161E-11
HepG2 24 hours	hsa-miR-769-5p	-1.058	8.79E-03	0.1859	FAXDC2	ENSG00000170271.11	1.565	1.249E-13
HepG2 24 hours	hsa-miR-4488	9.059	8.79E-03	0.1859	MELTF	ENSG00000163975.12	-1.085	2.722E-05
SNU 6 hours	hsa-miR-15a-5p	-1.526	9.20E-03	0.6536	CHAC1	ENSG00000128965.13	2.262	5.787E-28
SNU 6 hours	hsa-miR-12135	-1.298	2.11E-02	0.6536	FOXJ1	ENSG00000129654.8	1.240	2.099E-06
SNU 24 hours	hsa-miR-29b-1-5p	-2.126	2.50E-05	0.0007	HEY1	ENSG00000164683.17	1.692	2.304E-07
SNU 24 hours	hsa-miR-374c-5p	-1.884	5.25E-03	0.0787	CEBPB	ENSG00000172216.6	1.002	5.306E-17
SNU 24 hours	hsa-miR-374c-5p	-1.884	5.25E-03	0.0787	VEGFA	ENSG00000112715.22	1.448	2.131E-45
SNU 24 hours	hsa-miR-374c-5p	-1.884	5.25E-03	0.0787	ADM	ENSG00000148926.10	2.340	1.782E-143
SNU 24 hours	hsa-miR-32-3p	-1.821	9.41E-04	0.0175	YOD1	ENSG00000180667.10	1.658	1.063E-111
SNU 24 hours	hsa-miR-32-3p	-1.821	9.41E-04	0.0175	TNFAIP3	ENSG00000118503.15	2.316	3.471E-10
SNU 24 hours	hsa-let-7b-3p	-1.541	1.72E-02	0.1785	KLF4	ENSG00000136826.15	1.572	1.166E-31
SNU 24 hours	hsa-let-7b-3p	-1.541	1.72E-02	0.1785	ADRB2	ENSG00000169252.5	1.743	9.094E-48
SNU 24 hours	hsa-miR-195-5p	-1.456	1.51E-02	0.1770	VEGFA	ENSG00000112715.22	1.448	2.131E-45
SNU 24 hours	hsa-miR-195-5p	-1.456	1.51E-02	0.1770	ADRB2	ENSG00000169252.5	1.743	9.094E-48
SNU 24 hours	hsa-miR-34a-5p	-1.316	2.62E-02	0.2364	KLF4	ENSG00000136826.15	1.572	1.166E-31
SNU 24 hours	hsa-miR-34a-5p	-1.316	2.62E-02	0.2364	HSPA1B	ENSG00000204388.7	1.676	2.822E-82
SNU 24 hours	hsa-miR-5701	1.535	3.44E-02	0.2740	SYNGR3	ENSG00000127561.15	-1.505	1.627E-08
SNU 24 hours	hsa-miR-5701	1.535	3.44E-02	0.2740	RFLNB	ENSG00000183688.4	-1.340	1.864E-27
SNU 24 hours	hsa-miR-12136	3.334	4.93E-16	0.0000	ACVR2B	ENSG00000114739.14	-1.075	3.211E-13
SNU 24 hours	hsa-miR-12136	3.334	4.93E-16	0.0000	NR4A2	ENSG00000153234.14	-1.754	9.174E-06
SNU 24 hours	hsa-miR-12136	3.334	4.93E-16	0.0000	CSPG5	ENSG00000114646.10	-1.089	1.709E-08
SNU 24 hours	hsa-miR-12136	3.334	4.93E-16	0.0000	AREG	ENSG00000109321.11	-1.016	2.623E-03
SNU 24 hours	hsa-miR-12136	3.334	4.93E-16	0.0000	NRARP	ENSG00000198435.4	-1.763	3.743E-12
SNU 24 hours	hsa-miR-12136	3.334	4.93E-16	0.0000	MAF	ENSG00000178573.7	-1.748	5.311E-04
SNU 24 hours	hsa-miR-4492	8.427	4.75E-12	0.0000	IGFBP5	ENSG00000115461.5	-1.073	2.379E-18
SNU 24 hours	hsa-miR-4488	13.091	1.08E-32	0.0000	PDGFB	ENSG00000100311.17	-1.101	1.086E-51

The next step in the study consisted of predicting plausible interactions of lncRNAs with selected miRNAs in Figure 4C using DIANA-LncBase v3 (DIANA tools) (Figure 5A and B). From total predictions, only significantly deregulated lncRNA with coherent expression with target miRNAs were selected for further analysis (log₂FC > 1, P-adj < 0.05). No downregulated lncRNAs by Sorafenib were found to be coherently altered. Nineteen upregulated lncRNAs were selected, preferentially expressed at 24 hours of Sorafenib treatment (Figure 5C). CYTOR, DANT2, MALAT1, MIR4435-2HG, PCBP2-OT1, SNHG17, THAP9-AS1 and ZFAS1 were specifically regulated in HepG2 cells, whereas AZIN1-AS1, CARMN, EPB41 L4A-AS1, GAS5, LINC01578, LMCD1-AS1, PVT1, SNHG12, SNHG20, SNHG5 and SNHG7 were specifically related to SNU499 cells (Table 2). The role of lncRNA regulating mRNA instability was next explored. To complete the lncRNA-miRNA-mRNA network, we performed our in silico analysis using as targets those mRNAs that had been shown to be part of the miRNA-mRNA axes described in previous lines and their corresponding lncRNA extracted from miRNA-lncRNA signatures. The number of interactions at 24 hours of Sorafenib treatment is summarized in heatmaps in Figure 5D (HepG2 cells) and Figure 5E (SNU449 cells) (Supplementary Figure 3). In HepG2 cells, the most intense interaction was developed by the couple MALAT1-TFAP4, followed by MALAT- SMAD7 (Figure 5D), whereas in SNU499 cells PVT1-VEGFA or LINC01578-VEGFA were the most enriched pairs of lncRNA-mRNA (Figure 5E). The resulting lncRNA-miRNA-mRNA interaction networks were constructed accordingly (Figure 6), showing VEGA as a hotspot for ncRNA mediated regulation in SNU449 cells (Figure 6A).

Figure 5 Identification of lncRNA candidates. A and B: Gene regulatory networks showing miRNA-lncRNA interaction predictions in HepG2 and SNU449 cells treated with Sorafenib for 24 hours; C: Heatmap of selected lncRNAs in HepG2 and SNU449 cells for lncRNA-mRNA interactions; D and E: Heatmaps representing the number of interactions between lncRNA and mRNA targets of their corresponding miRNA interactions in HepG2 and SNU449 cells treated with Sorafenib (24 hours).

Figure 6 Gene regulatory networks of interactions among miRNAs (large green dots), lncRNAs (red squares) and mRNAs (small light blue dots) in HepG2 and SNU449 cells treated with Sorafenib (24 hours). A: HepG2 cells; B: SNU449 cells.

Table 2 Predicted interactions between miRNAs and mRNAs, and between their corresponding lncRNAs in HepG2 and SNU449 cells treated with Sorafenib for 6 and 24 hours.

Condition	miRNA	log2 (Fold-change)	Target	Target type	Transcript ID	log2 (Fold-change)	P-adj
HepG2 6 hours	hsa-miR-19b-3p	-1.277	RORA	mRNA	ENSG00000069667.16	1.048	1.51E-03
HepG2 6 hours	hsa-miR-19b-3p	-1.277	TNFAIP3	mRNA	ENSG00000118503.15	1.186	4.13E-07
HepG2 6 hours	hsa-miR-19b-3p	-1.277	ZNF385B	mRNA	ENSG00000144331.20	1.159	1.23E-02
HepG2 6 hours	hsa-miR-29b-3p	-1.374	HRK	mRNA	ENSG00000135116.9	1.967	8.93E-08
HepG2 6 hours	hsa-miR-29b-3p	-1.374	JMY	mRNA	ENSG00000152409.9	1.105	2.25E-13
HepG2 6 hours	hsa-miR-29b-3p	-1.374	TNFAIP3	mRNA	ENSG00000118503.15	1.186	4.13E-07
HepG2 24 hours	hsa-miR-193b-3p	-1.193	CYTOR	lncRNA	ENSG00000222041.11	1.421	7.45E-21
HepG2 24 hours	hsa-miR-193b-3p	-1.193	MIR4435-2HG	lncRNA	ENSG00000172965.16	1.204	7.19E-10
HepG2 24 hours	hsa-miR-193b-3p	-1.193	PCBP2-OT1	lncRNA	ENSG00000282977.1	1.086	8.35E-03
HepG2 24 hours	hsa-miR-193b-3p	-1.193	ZFAS1	lncRNA	ENSG00000177410.13	1.413	4.96E-35
HepG2 24 hours	hsa-miR-193b-3p	-1.193	SIAH1	mRNA	ENSG00000196470.12	1.581	1.14E-38
HepG2 24 hours	hsa-miR-193b-3p	-1.193	LAMC2	mRNA	ENSG00000058085.15	1.327	2.15E-03
HepG2 24 hours	hsa-miR-3187-3p	-1.169	MIR4435-2HG	lncRNA	ENSG00000172965.16	1.204	7.19E-10
HepG2 24 hours	hsa-miR-3187-3p	-1.169	DANT2	lncRNA	ENSG00000235244.4	1.014	1.37E-02
HepG2 24 hours	hsa-miR-3187-3p	-1.169	PLK2	mRNA	ENSG00000145632.15	1.634	3.01E-08
HepG2 24 hours	hsa-miR-3187-3p	-1.169	ATOH8	mRNA	ENSG00000168874.13	1.691	6.50E-08
HepG2 24 hours	hsa-miR-374a-3p	-1.327	MIR4435-2HG	lncRNA	ENSG00000172965.16	1.204	7.19E-10
HepG2 24 hours	hsa-miR-374a-3p	-1.327	OXR1	mRNA	ENSG00000164830.18	1.050	2.09E-08
HepG2 24 hours	hsa-miR-374a-3p	-1.327	GOLGA4	mRNA	ENSG00000144674.16	1.230	8.42E-03
HepG2 24 hours	hsa-miR-374a-3p	-1.327	ARRDC3	mRNA	ENSG00000113369.9	2.849	2.20E-41
HepG2 24 hours	hsa-miR-374b-3p	-1.064	THAP9-AS1	lncRNA	ENSG00000251022.6	1.954	2.61E-55
HepG2 24 hours	hsa-miR-374b-3p	-1.064	MALAT1	lncRNA	ENSG00000251562.8	1.024	2.38E-10
HepG2 24 hours	hsa-miR-374b-3p	-1.064	SMAD7	mRNA	ENSG00000101665.9	1.999	1.67E-42
HepG2 24 hours	hsa-miR-374b-3p	-1.064	TIRARP	mRNA	ENSG00000163659.13	1.070	1.55E-18
HepG2 24 hours	hsa-miR-769-5p	-1.058	MALAT1	lncRNA	ENSG00000251562.8	1.024	2.38E-10
HepG2 24 hours	hsa-miR-769-5p	-1.058	SNHG17	lncRNA	ENSG00000196756.13	1.093	1.46E-10
HepG2 24 hours	hsa-miR-769-5p	-1.058	TFAP4	mRNA	ENSG00000090447.12	1.010	1.16E-11
HepG2 24 hours	hsa-miR-769-5p	-1.058	FAXDC2	mRNA	ENSG00000170271.11	1.565	1.25E-13
HepG2 24 hours	hsa-miR-4521	-4.496	MBNL1	mRNA	ENSG00000152601.17	1.242	2.60E-16
HepG2 24 hours	hsa-miR-7974	-3.432	CREB3 L3	mRNA	ENSG00000060566.14	1.062	4.00E-06
HepG2 24 hours	hsa-miR-194-3p	-1.463	ARHGAP24	mRNA	ENSG00000138639.18	1.281	1.19E-08
HepG2 24 hours	hsa-miR-194-3p	-1.463	HOXA1	mRNA	ENSG00000105991.9	1.192	2.48E-06
HepG2 24 hours	hsa-miR-27b-5p	-1.236	OPTN	mRNA	ENSG00000123240.17	1.531	6.35E-20
HepG2 24 hours	hsa-miR-4488	9.059	MELTF	mRNA	ENSG00000163975.12	-1.085	2.72E-05
SNU 6 hours	hsa-miR-15a-5p	-1.526	CHAC1	mRNA	ENSG00000128965.13	2.262	5.79E-28
SNU 6 hours	hsa-miR-12135	-1.298	FOXJ1	mRNA	ENSG00000129654.8	1.240	2.10E-06
SNU 24 hours	hsa-miR-29b-1-5p	-2.126	HEY1	mRNA	ENSG00000164683.17	1.692	2.30E-07
SNU 24 hours	hsa-miR-374c-5p	-1.884	CEBPB	mRNA	ENSG00000172216.6	1.002	5.31E-17
SNU 24 hours	hsa-miR-374c-5p	-1.884	VEGFA	mRNA	ENSG00000112715.22	1.448	2.13E-45
SNU 24 hours	hsa-miR-374c-5p	-1.884	ADM	mRNA	ENSG00000148926.10	2.340	1.78E-143
SNU 24 hours	hsa-miR-374c-5p	-1.884	GAS5	lncRNA	ENSG00000234741.8	1.795	2.73E-202
SNU 24 hours	hsa-miR-374c-5p	-1.884	SNHG5	lncRNA	ENSG00000203875.12	1.370	1.13E-50
SNU 24 hours	hsa-miR-374c-5p	-1.884	SNHG20	lncRNA	ENSG00000234912.12	1.176	1.34E-23
SNU 24 hours	hsa-miR-32-3p	-1.821	YOD1	mRNA	ENSG00000180667.10	1.658	1.06E-111
SNU 24 hours	hsa-miR-32-3p	-1.821	TNFAIP3	mRNA	ENSG00000118503.15	2.316	3.47E-10
SNU 24 hours	hsa-let-7b-3p	-1.541	KLF4	mRNA	ENSG00000136826.15	1.572	1.17E-31
SNU 24 hours	hsa-let-7b-3p	-1.541	ADRB2	mRNA	ENSG00000169252.5	1.743	9.09E-48
SNU 24 hours	hsa-let-7b-5p	-1.541	SNHG12	lncRNA	ENSG00000197989.14	1.092	1.50E-21
SNU 24 hours	hsa-let-7b-5p	-1.541	LMCD1-AS1	lncRNA	ENSG00000227110.7	1.189	5.83E-03
SNU 24 hours	hsa-let-7b-5p	-1.541	CARMN	lncRNA	ENSG00000249669.10	1.637	1.11E-04
SNU 24 hours	hsa-miR-195-5p	-1.456	VEGFA	mRNA	ENSG00000112715.22	1.448	2.13E-45
SNU 24 hours	hsa-miR-195-5p	-1.456	ADRB2	mRNA	ENSG00000169252.5	1.743	9.09E-48
SNU 24 hours	hsa-miR-195-5p	-1.456	SNHG12	lncRNA	ENSG00000197989.14	1.092	1.50E-21
SNU 24 hours	hsa-miR-195-5p	-1.456	EPB41 L4A-AS1	lncRNA	ENSG00000224032.7	2.170	4.36E-34
SNU 24 hours	hsa-miR-195-5p	-1.456	AZIN1-AS1	lncRNA	ENSG00000253320.7	1.070	2.01E-03
SNU 24 hours	hsa-miR-195-5p	-1.456	PVT1	lncRNA	ENSG00000249859.11	1.264	4.13E-79
SNU 24 hours	hsa-miR-195-5p	-1.456	LINC01578	lncRNA	ENSG00000272888.7	1.446	4.59E-72
SNU 24 hours	hsa-miR-34a-5p	-1.316	KLF4	mRNA	ENSG00000136826.15	1.572	1.17E-31
SNU 24 hours	hsa-miR-34a-5p	-1.316	HSPA1B	mRNA	ENSG00000204388.7	1.676	2.82E-82
SNU 24 hours	hsa-miR-34a-5p	-1.316	SNHG7	lncRNA	ENSG00000233016.7	1.116	1.01E-40
SNU 24 hours	hsa-miR-5701	1.535	SYNGR3	mRNA	ENSG00000127561.15	-1.505	1.63E-08
SNU 24 hours	hsa-miR-5701	1.535	RFLNB	mRNA	ENSG00000183688.4	-1.340	1.86E-27
SNU 24 hours	hsa-miR-12136	3.334	ACVR2B	mRNA	ENSG00000114739.14	-1.075	3.21E-13
SNU 24 hours	hsa-miR-12136	3.334	NR4A2	mRNA	ENSG00000153234.14	-1.754	9.17E-06
SNU 24 hours	hsa-miR-12136	3.334	CSPG5	mRNA	ENSG00000114646.10	-1.089	1.71E-08
SNU 24 hours	hsa-miR-12136	3.334	AREG	mRNA	ENSG00000109321.11	-1.016	2.62E-03
SNU 24 hours	hsa-miR-12136	3.334	NRARP	mRNA	ENSG00000198435.4	-1.763	3.74E-12
SNU 24 hours	hsa-miR-12136	3.334	MAF	mRNA	ENSG00000178573.7	-1.748	5.31E-04
SNU 24 hours	hsa-miR-4492	8.427	IGFBP5	mRNA	ENSG00000115461.5	-1.073	2.38E-18
SNU 24 hours	hsa-miR-4488	13.091	PDGFB	mRNA	ENSG00000100311.17	-1.101	1.09E-51

Finally, potential targets between lncRNA and proteins under mRNA regulation were analyzed, and an integrated analysis was performed including deregulated mRNAs from Table 2 to study connected pathways by performing a STRING analysis (Figure 7A and B) and subsequent enrichment analysis (Figure 7C and D). In HepG2 cells (24 hours), SMAD dependent pathways were enriched (Figure 7A and C), whereas in SNU499 cells (24 hours), pathways related to growth factor receptors and angiogenesis were among the significant terms (Figure 7B and D). At 6 hours of treatment, although no lncRNA-miRNA-mRNA networks were highlighted, analysis was performed using mRNAs and targeted proteins. Enriched terms were related to apoptosis and regulation of transcription in HepG2 and SNU449 cells, respectively (Supplementary Figure 4).

Figure 7 Identification of protein interactions at 24 hours. A and B: Identification of protein interactions at 24 hours. STRING pathways of targets regulated by lncRNAs or miRNAs with no more than 10 interactions in HepG2 and SNU499 cells treated with Sorafenib (24 hours); C and D: Top 3 enriched Gene Ontology and Kyoto Encyclopedia of Genes and Genome terms in the STRING pathway analysis in HepG2 and SNU499 cells treated with Sorafenib (24 hours); E-G: Validation of proposed interaction networks using GEPIA data between MALAT1 and TFAP4, and MALAT and SMAD7 in HepG2 cells, and PVT1 and VEGFA in SNU449 cells.

Validation of the most prominent lncRNA-miRNA-mRNA was explored using RNAseq data by accessing GEPIA[29]. According to these in silico analysis, the axes MALAT1/hsa-mir-374b-3p/SMAD7 and MALAT1/hsa-mir-769-5p/TFAP4 were proposed to be of high relevance for Sorafenib response in HepG2 cells, whereas PVT1/hsa-miR-195-5p/VEGFA was related to the SNU449 cell line. A positive significant correlation was found between the expression levels of MALA and TFAP4 (Figure 7E), MALAT and SMAD7 (Figure 7F) (HepG2 cells), and PVT1 and VEGFA (Figure 7G) (SNU449 cells).

DISCUSSION

HCC is the most common form of primary liver cancer being globally recognized as the fourth cause of cancer-related death, and will lead to more than a million deaths annually in 2030[34]. In this sense, it is extremely important to identify biomarkers of disease. In particular, in the context of Sorafenib treatment, several lncRNAs and miRNAs have been related to the initiation, progression and metastasis of HCC[15]. miR-30e-3p, miR-30d, miR-31-5p and miR-221 are related to Sorafenib resistance in cultured renal and liver cancer cells[35-38]. Furthermore, Sorafenib increases the expression of miR-27a-3p, miR-122-5p, miR-193b-3p, miR-375 and miR-505-5p, but reduces the expression of miR-148b-3p, miR-194-5p, miR-200c-3p, miR-222-5p, miR-512-3p and miR-551a in HepG2 cells[18]. In this study, miR-200c-3p, miR-222-5p, and miR-512-3p appeared as prognostic miRNA markers of disease progression in two independent cohorts of patients under Sorafenib treatment[18]. Different lncRNAs such as NIFK-AS1, SNHG1, SNHG3, SNHG16, and NEAT1 have been related to Sorafenib resistance in liver cancer cells[39-43]. For instance, lncRNA PRR34-AS1 acts as a sponge of miR-296-5p and stimulates Wnt/beta-catenin signaling that leads to increased cell proliferation, migration, and invasion in cultured HCC cells, and facilitated tumor growth in vivo[44]. Similarly, the expression of RUSC1-AS1 negatively correlated with miR-340-5p expression in HCC cells that is tightly associated with the prevention of tumor cell proliferation and tumor progression[45].

In this study, we have identified lncRNA-miRNA-mRNA regulatory axes in HepG2 and SNU449 cells under Sorafenib treatment. We conducted an in silico analysis to identify interactions between RNA molecules and construct interaction networks based on bulk expression data that could help in deciphering treatment response in two different settings of Sorafenib response. However, it should be taken into account that HCC is a highly heterogeneous tumor that presents several layers of complexity, both inter and intra tumoral (differential genomic, pathology and molecular landscapes), that might be partially influenced by the wide diversity of etiologies[46]. In this sense, more clinically relevant studies as well as single-cell based experiments would be needed to confirm the relevance of predicted networks[47]. In HepG2 cells treated with Sorafenib we observed two key signatures based on the interaction of metastasis-associated lung adenocarcinoma transcript 1 (MALAT1)/hsa-miR-374b-3p/SMAD7 and MALAT1/hsa-miR-769-5p/TFAP4 (Figures 5 and 6; Table 3). Therefore, the lncRNA MALAT1 was of high importance for HepG2 cells during their response to Sorafenib. In fact, this lncRNA plays a very important role in normal physiology, as well as in hepatic fibrosis, hepatic carcinoma, liver regeneration, and fatty liver diseases[48]. In the context of Sorafenib response, it is involved in treatment resistance by targeting the miR-140-5p/Aurora-A signaling pathway[49]. Here, we identified mothers against SMAD7, TIRARP and TFAP4 as targets of MALAT1. The axis MALAT1/miR-769-5p/TFAP4 was the signature with the highest degree of lncRNA-mRNA interactions (Figure 5). TFAP4 belongs to the basic helix-loop-helix leucine-zipper (bHLH-LZ) family that has been found to be widely involved in the proliferation, differentiation, metastasis, angiogenesis, and other biological regulatory functions of malignancies[50]. However, it has also been shown to be a potent suppressor of c-MYC driven lymphoma[51]. In this regard, miR-769-5p should have anti-tumor properties. Several studies support its anti-oncogenic profile, by regulating transforming growth factor-β (TGF- β)[52] or growth factor receptors, among other pathways[53]. Moreover, the interaction between MALAT1 and miR-769-5p has also been experimentally tested[54]. MALAT1 was also able to interact with SMAD7 in this study (Table 3). Interestingly, SMAD7 inhibits TGF-β signaling that might impact various human diseases, such as tissue fibrosis, inflammatory disease as well as carcinogenesis[55]. In concordance with Figure 7, current literature supports the relevance of SMAD signaling in treatment response and within RNA networks. SMAD7 and MALAT1 were also able to interact with miR-374b-3p, another anti-tumor miRNA, that has been shown to regulate renal cell carcinoma[56]. In this complex network of interactions, MALAT1 was also shown to modulate TIRARP. TIRARP is a protein mono-ADP-ribosyl transferase that catalyzes the transfer of ADP-ribose from NAD⁺ to acceptor proteins involved in signal transduction, intracellular trafficking, cell metabolism and death[57]. Therefore, their role in cancer is in accordance with the results obtained in this study. Finally, SMAD7 and miR-374b-3p might also interact with the lncRNA THAP9-A1 (Table 3). Increased expression of THAP9-AS1 in HCC tissues and cells has been associated with tumor size, TNM stage and poor survival prognosis, and promoted HCC cell proliferation[58].

Table 3 Summary of lncRNA-miRNA-mRNA interactions in HepG2 and SNU449 cells at 24 hours of Sorafenib treatment.

lncRNA	miRNA	mRNA	Condition
MIR4435-2HG	hsa-miR-3187-3p	ATOH8	HepG2 24 hours
MALAT1	hsa-miR-769-5p	FAXDC2
ZFAS1	hsa-miR-193b-3p	SIAH1
THAP9-AS1	hsa-miR-374b-3p	SMAD7
MALAT1	hsa-miR-374b-3p
MALAT1	hsa-miR-769-5p	TFAP4
SNHG17	hsa-miR-769-5p
MALAT1	hsa-miR-374b-3p	TIPARP
SNHG12	hsa-let-7b-5p	ADRB2	SNU449 24 hours
EPB41 L4A-AS1	hsa-miR-195-5p
PVT1	hsa-miR-195-5p
LINC01578	hsa-miR-195-5p
GAS5	hsa-miR-374c-5p	CEBPB
SNHG5	hsa-miR-374c-5p
SNHG7	hsa-miR-34a-5p	KLF4
EPB41 L4A-AS1	hsa-miR-195-5p	VEGFA
PVT1	hsa-miR-195-5p
LINC01578	hsa-miR-195-5p
GAS5	hsa-miR-374c-5p
SNHG5	hsa-miR-374c-5p

The expression of ADRB2, CCAAT/CEBP, KLF4, and VEGFA were crucial for the response to Sorafenib in SNU449 cells. Of these, it was the regulation of VEGFA via the lncRNA and hsa-miR-195-5p axis which resulted in the higher degree of regulation (Figures 5 and 6; Table 3). PVT1 is known to be an oncogene that promotes tumorigenesis and treatment resistance by regulating YAP[59] or Wnt/beta-catenin[60] in different cancer contexts. In our setting, it was related to downregulation of hsa-miR-195-5p and upregulation of VEGFA. While this miRNA acts as an anti-tumor ncRNA in HCC[61], VEGFA is a strong immune suppressor that is used as a target in therapies such as Bevacizumab[62].

Regarding the other key players of Sorafenib response in SNU449 cells, KLF4 is a zinc-finger transcription factor expressed in differentiated epithelial cells in which it plays a relevant role in cell proliferation, differentiation, and apoptosis[63]. KLF4 upregulates CD9 and CD81 in exosomes that suppresses cell proliferation by negatively regulating MAPK/JNK signaling in HCC[64]. Curiously, KLF4 is regulated by has-let-7b-5p, a family of miRNAs already described to play a relevant role in differentiation and stemness[65]. No reports evaluating the role of ADRB2 in liver cancer exist. However, the expression of ADRB2 is related to a reduction of the invasiveness and motility of prostatic epithelial cancer cells[66].

CONCLUSION

In conclusion, this study identified differential lncRNA-miRNA regulatory axes in HepG2 and SNU449 liver cancer cells. They lead to increased expression of mRNAs with positive (SMAD7, TIRARP, TFAP4, FAXDC2 and ADRB2) and negative (VEGFA) therapeutic effects in HepG2 and SNU449 under Sorafenib treatment.