lncRNA–mRNA competing endogenous RNA network in IR-hepG2 cells ameliorated by APBBR decreasing ROS levels: a systematic analysis

Introduction

Transcriptome sequencing in various organisms has revealed that abundant transcription derived from human transcriptome generates a large proportion of non-coding RNAs (ncRNAs). Increasing evidence have shown the critical regulatory potential of ncRNAs in biological processes. The ncRNAs are classified into two categories according to their length: long non-coding RNAs (lncRNAs) with a length of more than 200 nucleotides and short noncoding RNAs (18 to 200 nucleotides), respectively (Li et al., 2019). lncRNA accounts for 68% of ncRNAs (Iyer et al., 2015), but they were considered “junk RNA” in the genome because they could not encode proteins (Zhang et al., 2019). However, recent studies have shown that lncRNA exert critical roles in regulating biological processes through comprehensive mechanisms, including genetic imprinting, immune response, tumorigenesis, cellular development and metabolism (Liyanage & Ganegoda, 2017). LncRNAs are also involved in the regulation of glucose and lipid metabolism in the liver, and their abnormal expression may be involved in the development of type 2 diabetes mellitus (T2DM) (Eslam, Valenti & Romeo, 2018).

Insulin resistance (IR)refers to a physiological condition in which target tissue such as fat, liver and muscle are unable to use insulin effectively. The metabolism of substances in the body is almost all regulated by insulin (Williams & Olefsky, 1990). Once IR occurs, the metabolism of sugars, lipids and proteins will be disturbed. The physiological regulation of insulin is mainly carried out by combining with the insulin receptors on the cells of various tissues. Any defect in the signal transmission in the process may affect the normal physiological regulation process of insulin and thus cause the occurrence of IR. IR degree in different tissues is specific, which makes its mechanism complex.

Reactive oxygen species (ROS) is a collective term to include oxygen-containing reactive species (Li, Jia & Trush, 2016). It is a significant contributor to the cell and tissue dysfunction in diabetes and the progression of complex diabetic complications (Newsholme et al., 2007). The increased ROS formation may cause oxidative stress, which plays a significant role in the development of hepatic IR (Liu et al., 2017). HepG2 cells are ideal model for studying hepatic IR pathogenesis and action mechanism of hypoglycemic drugs. The number of insulin receptors on the surface of HepG2 cells decreases with high level insulin intervention (Williams & Olefsky, 1990). At present, the underlying molecular mechanisms and regulation of ROS change in IR-HepG2 cells remains poorly understood.

Astragalus Polysaccharide (AP) is the main effective ingredient extracted from the Chinese herb Radix Astragali (Astragalus membranaceus var. mongholicus (Bunge)), which has the potential activity of antioxidant effects (Fu et al., 2018). Berberine (BBR) is an isoquinoline extracted from Coptis chinensis (Coptis chinensis var. angustiloba), which has many pharmacological activities such as anti-microbial, anti-diarrheal, reducing glucose and cholesterol, anti-tumor and immunosuppression. The plant names have been checked with http://www.theplantlist.org. Both of these two herbs have been used to treat diabetes for thousands of years in Traditional Chinese Medicine, study showed that the combination of Coptidis rhizoma, Astragali radix and Lonicera japonica(Lonicera alba L.) has a synergistic effect by increasing insulin sensitivity and ameliorating IR (Yang, 2017). In addition, it has been verified that AP and BBR could promote glucose up-take and improve carbohydrate metabolism respectively with different characteristics (Wang et al., 2019). Our previous study results showed that AP and BBR used together at a ratio of 1:1 could promote the basic secretion ability of IR-INS-1 cells model significantly, and showed a synergistic effect in diabetes treating (Mao, Shou & Chai, 2017). The study found has been applied for a patent, the application number: CN201810266347.9, publication/announcement number: CN108478591A. However, the cellular and molecular compatibility mechanism responsible for APBBR in attenuating IR remains unclear.

Additionally, there were few studies to identify key genes and pathways involved in ROS change in IR-HepG2 cells. Spurred by the development of deep sequencing technology, we are allowed to investigate the mechanisms of hepatic IR at the whole-transcriptome level. In the present study, we used high-throughput RNA-seq technology to explore the potential mechanism underlying hepatic IR at whole-transcriptome level by analyzing changes in the global gene expression profile in normal, IR and APBBR-treated HepG2 cells for the first time. First, differentially expressed genes (DEGs), related gene ontology (GO) and pathways were determined. Subsequently, gene act, pathway act and co-expression network were constructed to further explore the role of the ROS and IR related genes and pathways, which could shed new light to the underlying mechanisms of ROS and IR and therapeutic applications of APBBR in treating diabetes.

Materials & Methods

Results

High-dose insulin induced decrease of glucose content in IR-HepG2 cells

The results of intracellular glucose content determination in HepG2 cells are shown in Fig. 1A. Compared to control group, the glucose content in model group decreased significantly (^ΔΔP < 0.01.)

Expression profiles of mRNAs in experimental groups and ROS related mRNAs

Total RNA and whole transcriptome sequencing data (mRNA, miRNA, lncRNA) were acquired from qualified biospecimens (see “Methods”). The Illumina Hiseq2500 platform was used to obtain mRNA, lncRNA and miRNA data. As a result of this analysis, we detect thousands of genes. To determine the regulation of mRNA expression, we performed an unsupervised clustering analysis of the significantly regulated genes in IR and APBBR cells. The EBseq algorithm was applied to filter the differentially expressed genes(DEGs) according to the criteria of a—log2 (fold change)— ≥ 1and a P-value <0.05. There were 157 genes whose changes met the criteria. 50 genes were significantly and differentially expressed in Model group relative to Control group, with 24 up-regulated genes and 26 down-regulated genes. 107 genes were significantly altered in APBBR group compared with that in Model group, with 56 genes up-regulated and 51 genes down-regulated Fig. 2A. The Volcano plots for mRNAs are showed as Figs. 2B, 2C. The heatmaps for DEGs are showed as Figs. 2D, 2E. To analyze the function of these DEGs, we next conducted GO Figs. 2F, 2G, and KEGG pathway analysis on these genes Fig. 2H. DEGs in Model vs. Control were largely involved in lysine degradation, RNA degradation, taurine and hypotaurine metabolism, sulfur relay system, ascorbate and aldarate metabolism, AMPK signaling pathway. DEGs in APBBR vs. Model were involved in Protein processing in endoplasmic reticulum, spliceosome, systemic lupus erythematosus, steroid hormone biosynthesis, tryptophan metabolism, MAPK signaling pathway. According to GO and KEGG pathway results as well as previous reports of IR, we notice that genes such as SCD, DUSP5 and PLK3were quite important for ROS development, which may result in IR Fig. 2I.

Differentially expressed lncRNAs in experimental groups and lncRNAs that participates in APBBR treatment

A total of 66,980 lncRNAs were annotated in this study. In this study, 12,962 lncRNAs were first annotated as known lncRNAs (database: GRCh38 Ensembl v88).The 32,433 newly identified lncRNAs included 17,508 lncRNAs (category code u), 85 antisense lncRNAs, 526 intron lncRNAs, and 3466 other lncRNAs. Almost all known lncRNAs are distributed on each chromosome, but those universal exon lncRNAs and reference lncRNA transcriptors (category code o) do not show significant chromosome location preference Fig. 3A.

From the RNA-seq results, a total of 521 lncRNAs were detected with P < 0.05.Compared to Control, 384 lncRNAs were detected to be differentially expressed in Model (89 up regulated and 295 down regulated). Compared to Model, 196 lncRNAs were identified to be differentially expressed in APBBR (113 up regulated and 83 downregulated) Fig. 3B. The Volcano plot for lncRNAs are showed respectively as Figs. 3C, 3D. Heatmap for top 100 differentially expressed lncRNAs in different groups are showed as Figs. 3E, 3F.

If APBBR can improve IR, it is necessary to ensure that the lncRNA expression trend is exactly opposite in both Model vs. Control and APBBR vs. Model. That is to say, we should find out genes decreased in Model group compared to Control group but increase in APBBR group contrasted with Model group (decrease-increase type) and genes increased in Model group compared to Control group but decreased in APBBR group in contrast to Model group (increase-decrease type). According to this principle, Venn analysis was additionally applied to learn the lncRNAs that participates in the exertion of APBBR treatment effect, including increasing/decreasing type and decreasing/increasing type. We found two of them(CTD-2600O9.2, MIR4435-2HG) can be reversed by APBBR in the 295 downregulated genes following IR. Similarly, APBBR treatment tend to reversed five (RP5-1057I20.5, CTD-3014M21.1, HOXB-AS3, LAMA5-AS1, CTD-2517M22.14) of the 89 upregulated lncRNAs following IR (Figs. 3G, 3H).

Comparative characteristics analysis of mRNAs and lncRNAs

To investigate the difference between differentially expressed mRNAs and lncRNAs in genomic characteristics, transcript lengths, the gene lengths, exons number, ORF lengths, isoforms number were analyzed. FPKM data shows that the abundance of lncRNA is lower than that of mRNA in RNA-seq samples, indicating that lncRNAs exhibited a notable bias toward fewer transcript isoforms Figs. 4A, 4B. In addition, we found that the length of lncRNA is usually shorter than that of mRNA Fig. 4C. For example, most lncRNAs have less than six exons, while mRNA have more exons and exon numbers distributed over a wider range. Some mRNA have as many as 30 exons Fig. 4D. LncRNAs also have shorter (60–90 nucleotide) ORF than mRNA, while most mRNAs have more than 500 nucleotide ORF Figs. 4E, 4F.

Construction of a lncRNAs–miRNA–mRNA network

To systematically explore the influence of dynamic changes in ceRNA regulation on gene expression related to IR and ROS in Normal, Model and APBBR-treated hepG2 cells. we constructed a ceRNA network, integrating matched expression profiles of lncRNAs, miRNAs and mRNAs. Target genes were selected according to the results of GO or KEGG pathway analysis or according to previous reports of ROS-associated genes that may lead to IR, and then the global miRNA-mediated ceRNA network was constructed (Fig. 5). In the lncRNA-miRNA-mRNA predicted interaction network, seven lncRNAs and three mRNAs had the same expression trend (belonging to the increase-decrease type/decrease-increase type). Twenty-five miRNAs were predicted to be possible target genes for regulation.

Discussion

Systemic or local IR occurs when insulin target organs or tissues become less sensitive and responsive to endogenous or exogenous insulin (Semenkovich, 2016). IR refers to the weakened response of tissues to circulating insulin (Archer, Von Schulze & Geiger, 2018), it mainly take place in skeletal muscle cells, fat cells and liver cells, it can also occur in vascular endothelium and islet beta cells (Dontsov & Vasil’eva, 2016). IR involving multiple molecules and signaling mechanism, including insulin and its antagonist, insulin receptor substrates, phosphatidyl inositol 3 kinase (PI3K), glucose transporter gene and protein, mitogen-activated protein kinase (MAPK), etc (Seino et al., 2010). It has been proved that IR is the common pathological basis of abnormal glucose and lipid metabolism of essential hypertension, T2DM,coronary heart disease, hyperlipidemia and other diseases (Tan, Sasagawa & Mori, 2017; Otto et al., 2019).

ROS include superoxide ionon (O²⁻), hydroxyl radical (OH⁻) and hydrogen peroxide (H₂O₂), among which O²⁻ has high activity and cytotoxicity and is mainly produced by aerobic microorganisms (Newsholme et al., 2016). If the cell antioxidant system fails to suppress ROS, ROS would react with cell macromolecules, leading to lipid peroxidation, causing cell DNA damage, which may result in oxidative stress and oxidative damage to the cell (Chattopadhyay et al., 2015). Oxidative stress is also known as ROS imbalance, where reactive oxygen species produced by the body exceeds their antioxidant capacity, leading to excessive generation of free ROS or dysfunction of antioxidant system (Meex, Blaak & Van Loon, 2019). Oxidative stress can damage critical cellular macromolecules and stimulate various stress sensitive intracellular pathways such as c-Jun, N-terminal kinase, ERK1/2 and the NF-κB (Keane et al., 2015). It can also produce chronic low-grade inflammatory reaction (Klaunig et al., 2011). It is increasingly accepted that oxidative stress is a significant contributor to the progression of IR.

It is noteworthy that IR also produces excessive free radicals and superoxides through proton electrochemical gradient, which destroys the antioxidant defense ability of various tissues (Valko et al., 2007), and also affects the expression of glucose transporter-GLUT1 and antioxidant defense function, thus increasing oxidative stress (Gonzalez-Menendez et al., 2018).

Astragalus Polysaccharide (AP) and berberine (BBR) have been reported to have antidiabetic properties (Meng et al., 2018; Jiang et al., 2015). AP can modulate the insulin-initiated phosphorylation cascades in a similar manner as metformin in HepG2/IR cell model (Sun et al., 2019). BBR is reported to demonstrate an excellent activity with oxidative stress and inflammation through a variety of signaling pathways including NF-κB, AMPK, Nrf2/HO, JNK, c-Jun and MAPKs pathways and various kinases in cell (Ma et al., 2018) which indicates its IR-improving potential. However, the molecular mechanisms for this action are not fully clarified. Our study showed that APBBR can improve the proliferation activity and decrease the level of ROS in IR-HepG2 cell model, suggesting the antioxidant mechanism of APBBR and APBBR may be a stable and reliable means to attenuate IR and its adverse metabolic consequences in HepG2 cells.

RNA-seq analysis is a powerful approach for investigating underlying molecular mechanisms of complex biological processes (Han et al., 2014). NcRNAs include miRNAs and LncRNAs. LncRNA is a hotspot in the field of biological medicine, its important function in the cell has been widely accepted as the successor of another wave of research followed miRNAs. It was found that the expression of lncRNA in human pancreatic island cells had exceeded 1,100 species (Moran et al., 2012). Various miRNAs and lncRNAs have been found to be associated with IR and intensive studies on ncRNAs had shed light on understanding the pathogenesis of IR (Esguerra et al., 2011; Davalos et al., 2011; Zeggini et al., 2007). Recently, competing endogenous RNA networks have been reported as a novel mechanism to set up an extensive regulatory network in a post-transcriptional layer. According to the “ceRNA hypothesis”, lncRNAs competitively interacted with the miRNAs and thereby “talk” to mRNAs through microRNA response elements (MREs) (Phelps et al., 2016). Based on this reasoning, we constructed a global triple network ,in which lncRNA and mRNA from a triple sharing the same miRNA.

To the best of our knowledge, this was the first study to detect genome-wide transcriptome of IR-HepG2 cell to predict potential lncRNA-mediated networks in APBBR attenuating IR at the molecular level. In the present study, we found that compared with Control group, 384 DElncRNAs and 50 DEGs were present in Model group and there were 196 DElncRNAs and 107 DEGs in APBBR vs. Model. GO and KEGG pathway analyses showed the potential functions of the DEGs involved in the mechanism of APBBR treatment; and the ceRNA analysis identified several potential lncRNA/miRNA/mRNA interaction networks that might ameliorate the development of IR. The sequencing depth involves mRNA, miRNA and lncRNA, which allowed us to detect differences in gene expression at the transcriptome level and predict the regulatory mechanisms that may mediate ROS changes by ceRNA. These results suggest more comprehensive perspectives to understand IR progress and explore mechanisms of APBBR treatment. Firstly, we select the candidate mRNA according to GO and KEGG pathway analyses and previous relevant studies of ROS and IR. Secondly, for validation, we chose the lncRNAs that exhibited opposite variation trend before and after APBBR treatment. Venn analysis was additionally applied to learn the possible marker that participates in the process of APBBR’s improving IR.

The ceRNA network offers considerable clues for understanding the key roles of ceRNA-mediated gene regulatory networks in ROS genesis and APBBR regulation mechanism as well as further detail analysis of biological function of lncRNA. From the ceRNA network, we identified a substantial amount of cross-talk within the non-coding RNAs, in which mRNAs were co-expressed with candidate lncRNAs mediated by miRNAs and subsequently formed a complex network. According to this analysis, we discovered 7 lncRNAs involving in the process of APBBR’s improving IR. Of which, the lncRNA MIR4435-2HG decreased expression in Model group but increased in APBBR group, connecting many mRNAs that may ameliorate IR according to GO enrichment analysis and pathway analysis, and it could directly interact with aberrant miRNA expression contributes to IR. For instance, miR-320 is reported to promote IR in high glucose treated adipocytes and anti-miR-320 oligo was found to regulate IR in adipocytes by improving insulin–PI3-K signalling pathways (Ling et al., 2009). In our study, miR-320 increased expression in Model group but decreased after APBBR treatment. According to the ceRNA theory, there may be a triplet that contained MIR4435-2HG, miR-320 and the downstream target mRNAs, the lncRNA could sponge the miRNA and change the expression of target mRNAs. We found that PLK3 and DUSP5 tend to be the target mRNAs in the network. Polo-like kinase 3 (Plk3) is a serine/threonine protein kinase of the Polo-like kinase family, which is critical for tumor suppression and stress responses (Barr, Sillje & Nigg, 2004). Moreover, up-regulation of Dusp5 (a tumor suppressor) is able to negatively regulate the MAPK signaling pathway (Zhang et al., 2014). MAPK plays a critical role during the process of inflammation and IR (Zhou et al., 2018), the negative regulation of MAPK might ameliorate inflammation and activate insulin signaling in liver. That is to say, miR-320 may cause IR by activating the MAPK pathway to increase hepatic oxidative stress and APBBR may decrease ROS level in IR-HepG2 cells by negatively regulating miR-320.

IR can play a role in lipid metabolism disorders in the body and promotes lipid deposition in liver. We notice that the mRNA related to lipid metabolism SCD increase significantly in IR vs. Control and it tend to be decreased in APBBR vs. IR. SCD encodes an enzyme involved in fatty acid biosynthesis, primarily the synthesis of oleic acid. The enzyme is also important in the conversion of polyunsaturated fatty acids (PUFAS) to monounsaturated fatty acids (MUFAS). The PUFAs are the major cellular target of ROS and readily provide an extractable proton to oxygen-free radicals (Ushio-Fukai & Alexander, 2004). The decreased expression of SCD after APBBR treatment may decrease ROS level in IR-HepG2 cells and thereby attenuating IR. These results are consistent with previous studies that SCD has been shown to be over expressed in fat tissues with obesity and metabolic disorders and absence of SCD can improve metabolic syndrome in mice (MacDonald et al., 2008). Using ceRNA, we narrowed the targets of mRNAs related to IR and ROS with a lncRNA_miRNA_mRNA network, pointing the way for further studies to verify the genes and describe their function.

Generally, we believe that biochemical experiments combined with high-throughput technologies could shed new light to explore complex biological mechanisms underlying the efficacy of APBBR. At the same time, the lncRNA played an important role in the IR and could offer new therapeutic targets for IR mechanism research. It is noteworthy that many lncRNAs in the ceRNA network are not annotated, which is quite worth futher study. We hope to inspire researchers to study the role of non-coding RNAs in IR and APBBR treatment.

Conclusions

In summary, our investigation indicates that attenuated IR in hepG2 cells manifestation by increased cell survival rate and decreased ROS level. Furthermore, we explored the target genes highly enriched in pathways related to ROS and IR and predicted the network regulated by lncRNA according to ceRNA theory, revealing the potential biofunctional roles for lncRNAs as drivers of IR. Transcriptome profiling and regulatory network provides a rationale for exploiting the insulin-sensitizing potential of APBBR by decreasing ROS level via MAPK pathway, which enriches the therapeutic options in the treatment of IR and enhanced our understanding of APBBR. This study provided the expression profile of lncRNAs in IR-hepG2 cells for the first time and laid a foundation for the study on the molecular mechanism of lncRNA regulating IR as well as the follow-up study on the lncRNA regulating mechanism of APBBR in treating IR.

Supplemental Information