Qianggan extract improved nonalcoholic steatohepatitis by modulating lncRNA/circRNA immune ceRNA networks

Background The traditional Chinese medicine prescription, Qianggan formula have been confirmed to be effective on non-alcoholic steatohepatitis (NASH), however, the underlying molecular mechanisms remain obscure. Methods Thirty-six male C57BL/6 mice were randomly divided into three groups: normal chow diet group; methionine-and-choline-deficient diet (MCD) group, and Qianggan extract (QG) intervention group (0.4 g/kg daily) that fed with MCD. The efficacy of QG was biochemically and histologically evaluated. The expression profiles of messenger ribonucleic acids (mRNAs), long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) were examined using microarray and verified by RT-qPCR. Results QG significantly improved the phenotypic characteristics of NASH, as serum alanine aminotransferase (ALT), aspartate aminotransferase (AST), alkaline phosphatase (ALP), and lactate dehydrogenase (LDH) levels and liver inflammatory cytokines were significantly decreased. By the cutoff of a 1.5-fold change and P < 0.05, 6193 mRNAs, 5692 lncRNAs and 4843 circRNAs were identified as differentially expressed between the MCD and normal groups, and 514 mRNAs, 1182 lncRNAs and 443 circRNAs were identified as differentially expressed between the QG and MCD groups. The intersections (244 mRNAs, 259 lncRNAs and 98 circRNAs) among the three groups were chosen for analysis. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment revealed that most overlapping mRNAs were related to immune functions such as natural-killer-cell-mediated cytotoxicity, intestinal immune network for IgA production, and T cell receptor signaling pathway. Pathway interactions, protein-protein interactions and molecular complex detection (MCODE) analysis identified numerous immune-related hub genes e.g. natural cytotoxicity triggering receptor 1(Ncr1), C-X-C motif chemokine ligand 9 (Cxcl9), Klra1, and Cd28. Finally, two lncRNAs (Sngh1 and Slc36a3os) and four circRNAs (circ_0009029, circ_0004572, circ_0009212 and circ_0009453) in competing endogenous RNA (ceRNA) networks were constructed by Cytoscape, and immune-related mRNAs (e.g., Cd28, Cd8a, Il15, and Klrk1) were involved in the ceRNA networks. Conclusions LncRNA and circRNA-associated immune ceRNA networks might be the targets of QG in alleviating NASH, and our work may provide valuable clues for exploring the mechanisms underlying the effect of QG. Electronic supplementary material The online version of this article (10.1186/s12906-019-2577-6) contains supplementary material, which is available to authorized users.


Background
Nonalcoholic fatty liver disease (NAFLD) refers to a broad spectrum of diseases characterized by fat infiltration in the liver [1]. Non-alcoholic steatohepatitis (NASH) is a progressive form of NAFLD characterized by lobular inflammation, hepatocellular ballooning, and fibrosis with an inherent risk of progression to end-stage hepatocellular carcinoma [2,3]. Twenty percent of NASH patients are reported to develop cirrhosis, and 30-40% of patients with NASH cirrhosis die from liver-related diseases, NASH has become the third most common indication for liver transplantation in the United States [4].
The pathogenesis of NASH has not been fully elucidated, and the mechanisms appear to be multifactorial. "Two hits" and "multiple hits" hypotheses have been proposed to describe the pathogenesis of NASH [5,6]. Insulin resistance, inflammatory cytokines, and oxidative stress are thought to be important in the progression of NASH. Pharmacological agents such as insulin sensitizers, antioxidants, and lipid-lowering drugs have been implicated to improve NASH characteristics and progression [7]. However, there is no consensus regarding the effective and appropriate drug therapy for NASH.
Coding ribonucleic acids (RNAs) coupled with noncoding (nc) RNAs are involved in NASH pathogenesis [10]. MicroRNAs (miRNAs), long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs) are different species of ncRNAs [11]. MiRNAs are 18-25 nucleotides long, evolutionarily conserved, single-stranded RNA, that negatively modulate the expression of their target mRNAs [12]. Numerous miRNAs have been reported to be correlated with NASH [13,14]. LncRNAs are RNA transcripts > 200 nucleotides, that can regulate gene expression through many patterns, including chromosome remodeling, and transcriptional and post-transcriptional processing. LncRNAs are involved in the pathogenesis of liver diseases and have potential diagnostic, prognostic, and therapeutic importance [15], and contribute to the development of NASH and fibrosis [10,16]. LncRNAs can act as competing endogenous RNAs (ceRNAs) or miRNA sponges, interacting with miRNAs in a manner that sequesters these molecules and reduce their regulatory effect on target mRNAs [17]. Similar to lncRNAs, circRNAs (a class of RNA molecules that lack 5′-3′ ends and poly A tails and covalently form closed loops) can also function as ceRNAs or miRNA sponges [18]. The importance of ceRNAs has been addressed in various diseases. A deep RNA sequencing study constructed comprehensive circRNA-associated ceRNA networks, and inferred the diagnosis and therapy of Alzheimer's disease [19]. Another study provided a lncRNA-related ceRNA network and suggested two-lncRNA signatures as independent biomarkers for the prognosis of lung squamous cell carcinoma [20].
To systemically understand the efficacy of Qianggan extract (QG) on NASH, a comprehensive microarray analysis of mRNAs, lncRNAs and circRNAs was conducted to obtain expression profiles. By identifying differentially expressed mRNAs, lncRNAs and circRNAs and constructing ceRNA networks, our data shed new light on potential novel biomarkers and drug targets of NASH.

Preparation of Qianggan extract
The Qianggan formula was composed of 16 ingredients (Table 1), and the extract was prepared using the following procedures: all of the ingredients were soaked in 10 times the volume of water and boiled for 2 h, the supernatant were collected, 8 times volume of water was added and the mixture was boiled another 1.5 h. Then, the supernatant was collected, 6 times the water volume was added, and the mixture was and boiled 1.5 h for the third time. All the supernatants were filtered, the PH was adjusted to 8. The extract powder of was then obtained by concentrating the solution to a density ratio of 1.35.

Animal treatment and sample collection
Six-week-old male C57BL/6 mice were obtained from Shanghai Laboratory Animal Co. Ltd., and maintained in a specific pathogen-free environment with temperatureand humidity-controlled conditions. Thirty-six mice were allocated to six per cage, and randomly divided into three groups (n = 12 per group) based on bodyweight: the normal group, fed a chow diet; the methionine-andcholine-deficient diet (MCD) group, fed a MCD diet; and the QG group, fed a MCD diet and supplemented with the extract (Orient Pharmaceutical Co., Ltd., 0.4 g /kg daily). The extract was dissolved in distilled water and intragastrically administered (1 ml/100 g body weight) for 4 consecutive weeks, and the other mice were given an equal volume of distilled water. Body weight was recorded once per week. At the end of the experiment, mice were weighed at 8 am the next morning and anesthetized by intraperitoneal injection with 2% pentobarbital sodium (1.5 ml/kg bodyweight). Then, blood was collected by cardiac puncture, centrifuged at 3000 rpm for 10 min at 4°C, and serum was separated. Liver tissue was quickly removed, rinsed with 0.9% sodium chloride solution and weighed. Liver tissues from the identical lobe and position were cut and then fixed in 10% neutral-buffered formalin. The remainder was stored at − 80°C after snap freezing in liquid nitrogen. After sample collection, the animals were physically euthanized by quickly breaking the spines to avoid insufficient anesthetization, and bodies were disposed of by the central department of animal care.
All animal procedures were approved by the Animal Experiment Ethics Committee of Shanghai University of Traditional Chinese Medicine (approval number: 201703014).

Determination of the serum enzymes levels
Serum alanine aminotransferase (ALT), aspartate aminotransferase (AST), alkaline phosphatase (ALP), and lactate dehydrogenase (LDH) were analyzed using the Hitachi full-automatic system with the corresponding kits (Wako, Richmond, VA, USA).

Histology staining
Liver samples were fixed in 10% formalin for 48 h, embedded in paraffin, sectioned to 4 μm thickness, stained with hematoxylin and eosin (H&E), and examined under a light microscope at 200 × magnification.

Microarray detection and analysis RNA preparation and microarray detection
Liver total RNA was extracted using TaKaRa RNAiso Plus regent and RNA integrity was assessed using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, US). RNA was further purified by an RNeasy mini kit (QIAGEN, Hilden, Germany) and RNase-Free DNase Set (QIAGEN), amplified and labeled by Low Input Quick Amp WT Labeling Kit (Agilent Technologies), purified and hybridized with a SBC-Mouse ceRNA microarray (SBC, Shanghai, China). Nine slides were scanned by an Agilent Microarray Scanner with default settings: dye channel, green; scan resolution, 3 μm; PMT 100%; 20 bit. Data were extracted with Feature Extraction version 10.7 (Agilent Technologies). Raw data were normalized by the Quantile algorithm, limma packages in R.

Microarray data processing
The statistical significance of differentially expressed mRNAs, lncRNAs and circRNAs among groups (MCD vs normal and QG vs MCD) was identified by the cutoff of 1.5-fold change and P < 0.05. The intersections of differentially expressed mRNAs, lncRNAs and circRNAs between MCD vs normal and QG vs MCD were determined by Venn diagrams. The overlapped mRNAs, lncRNAs and circRNAs were subjected to hierarchical clustering by Cluster version 3.0 to generate an overview of the characteristics of expression profiles.

Gene ontology (GO) and pathway analysis
To explore the main biological function of differentially expressed genes, functional enrichment was performed based on GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) database [21,22]. GO classifies gene functions according to molecular function, biological process and cellular component. Pathway enrichment analysis revealed the main metabolic and signaling pathways. The enriched pathways interaction networks were constructed using the Cytoscape plug-in ClueGO and CluePedia.
Protein-protein interaction (PPI) network analysis PPI network analysis was performed based on the STRING database (https://string-db.org/). The network was visualized by Cytoscape, and hub genes were obtained by screening the degree of connectivity of each node in the network. In the network, nodes represent genes and lines represent the interactions. To identify the significant modules in the network, Cytoscape plugin MCODE (Molecular Complex Detection) was conducted with a score > 4.

LncRNA and circRNA associated ceRNA network construction
According to the ceRNA hypothesis, ceRNA members can compete for the same miRNA response elements (MREs) to regulate each other. Various types of RNA transcripts compete with one another for binding to a shared miRNA, thereby negatively regulating miRNAmediated gene silencing [23]. Hence, the expression pattern of mRNAs should be in the same direction as lncRNAs or circRNAs. The potential miRNA-binding sites were searched by the sequences of lncRNAs, circRNAs and mRNAs based on the TargetScan database (http:// www.targetscan.org/). Then, lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks were constructed based on lncRNA/circRNA-miRNA and miRNA-mRNA regulation pairs. In addition to overlapping miRNA binding, mRNAs with an expression pattern in the same direction as lncRNAs or circRNAs were filtered out to construct lncRNA/circRNA -miRNA-mRNA ceRNA networks by Cytoscape.

Quantitative real-time PCR assay
To validate our microarray data, quantitative Real-time experiments were performed as previously described [24]. We selected 13 representative RNAs to verify the expression, and primers sequences were listed in Table 2.

Statistical analysis
Data were expressed as mean ± SD and were analyzed by one-way analysis of variance by SPSS version 18.0. P < 0.05 was considered significant different.

Phenotypic characteristics
Four weeks of MCD feeding induced typical NASH features, and H&E staining showed significant hepatic steatosis and inflammation in MCD mice compared to normal mice (Fig. 1h). Serum AST (Fig. 1d), ALT (Fig.  1e), ALP (Fig. 1f ) and LDH (Fig. 1g) levels were all markedly increased in MCD mice. QG extract supplementation improved hepatic steatosis and inflammation, and restored increased serum parameters (Fig. 1d-g).
The body weight and liver weight were decreased in the MCD mice in comparison to normal mice, while QG had no effect on these parameters (Fig. 1a&b). The ratio of liver weight and body weight among the groups revealed no significant difference (Fig. 1c).
To obtain an overview of the expression pattern of DEGs, DElncRNAs and DEcircRNAs among groups, hierarchical clustering was performed. Most of the DEGs, DElncRNAs and DEcircRNAs exhibited opposite expression patterns between MCD vs normal and QG vs MCD (Fig. 2d-f ). For instance, the DEGs Cxcl9, Cd28, Cd8a, Btla, Klra1, Klrk1 and Ncr1were up-regulated in the MCD vs normal group, but down-regulated in the QG vs MCD group. The DElncRNAs Snhg1 and Slc36a3os were increased in the MCD vs normal group, but decreased in the QG vs MCD group. DEcircRNAs-circ_0009029 and circ_0004572 were induced in the MCD vs normal group, but reduced in the QG vs MCD group. These results indicated that Qianggan extract could affect the gene expression of NASH mice.

GO and pathway analysis
To further explore the main biological functions of DEGs, GO and pathway enrichments were performed (Additional file 4: Table S4). The 244 overlapping DEGs were enriched in immune-related GO terms and KEGG pathways such as natural killer cell-mediated cytotoxicity, intestinal immune network for IgA production, and T-cell receptor signaling pathways. Figure 3a and b lists the top 30 enriched GO terms and KEGG pathways. To identify the key genes in the pathways, enriched KEGG pathway interactions were constructed, and hub genes such as Cd28, Klra1, and Klra7 were indicated (Fig. 3c).

PPI network analysis
PPI network analysis and MCODE were conducted to depict interactions and significant modules among the overlapped 244 DEGs. A PPI network with 118 nodes (average of 2.66 neighbors) was obtained (Fig. 4a). The important hub genes with no less than five neighbors were identified. By MCODE, two significant modules were filtered out with a score > 4. Most hub genes, such as Cd28, Cd8a, Cxcl9, Klka1, Klrk1and Ncr1c, and the modules were immune related.

LncRNA and circRNA associated immune ceRNA networks
LncRNAs and circRNAs interacting as ceRNAs play important roles in regulating gene expression. However, the biological functions of most lncRNAs and circRNAs are unknown. To unravel the functions of identified DElncRNAs and DEcircRNAs, we constructed ceRNA networks. Interactions of miRNAs with mRNAs, lncRNAs and circRNA were predicted based on the Targetscan database (Additional file 5: Table S5). By combining of literature retrieval and the ceRNA hypothesis, we constructed two lncRNA-and four circRNA-associated immune ceRNA networks that were implied in NASH.
In addition, circ_0004572 was a ceRNA of 2 miRNAs that targeted 32 mRNAs (Fig. 6a). Circ_0009453 was a ceRNA of four miRNAs that targeted 54 mRNAs (Fig. 6c). Circ_0009212 was a ceRNA of 4 miRNAs that targeted 68 mRNAs (Fig. 6e). Circ_0009029 was a ceRNA of miR-362-5p that targeted 24 mRNAs (Fig. 6g). GO analysis demonstrated that mRNAs involved in the networks were enriched in immune-related terms (Fig. 6b, d, f and h). We obtained circRNA-associated immune ceRNA networks that comprised 4 circRNAs, 11 miRNAs and 100 m RNAs, and we also noticed that the four circRNAs shared four mRNAs (Fig. 6i and j).

Validation of the RNA expression
To validate our findings based on microarray data, we selected 13 RNAs for verification by RT-qPCR. Fold changes of RNAs (MCD vs normal and QG vs MCD) were calculated based on the average expression values. As shown in Fig. 7, most of the results from RT-qPCR were in consistence with microarray data, which might suggest the reliability of our microarray analysis and findings.

Discussion
Here we confirmed the effect of QG on NASH mice, and comprehensive microarray analysis revealed a batch of differentially expressed mRNAs, lncRNAs and cir-cRNAs. Through bioinformatics analysis, we constructed lncRNA-circRNA-associated immune ceRNA networks, and obtained certain immune-related biological processes, pathways and hub genes that responded QG intervention.
Qianggan formula was initially designed for treat hepatitis B virus in China, and it showed obvious effect in reducing plasma ALT, AST and inflammatory parameters [25]. In clinical trials of NAFLD, Qianggan formula treatment significantly reduced hepatic lipid content in patients, as evidenced by ultrasound test. In addition, plasma TG and ALT levels were decreased upon treatment [26]. In high-fat diet-induced NAFLD rats, 4 weeks of QG treatment prevented leptin resistance, and alleviated lipid accumulation and inflammation in the liver [27]. Because the formula is a mixture of 16 ingredients and undoubtedly has multiple targets, integrative highthroughput analysis could facilitate the exploration of   Immune responses are involved in the development of NASH [28]. Many types of innate (NK cells, NKT cells, and kupffer cells/macrophages) and adaptive (T cells and B cells) immune cells are enriched in the liver and play pivotal roles in NASH progression [29]. NK-cellmediated cytotoxicity and the T-cell receptor signaling pathway are important immune processes, and studies have shown that the dysfunction of these pathways is involved in liver inflammation. Accordingly, we identified changes in several immune hub genes, such as Cd28, Cd8a, Cxcl9, Klra1, Klrk1, Klra9, Il15, and Ncr1that related to the effect of Qianggan extract.
LncRNAs might have important roles in NAFLD/ NASH [30]. LncRNA SNHG1 is up-regulated by MYCN amplification and could be a potential prognostic biomarker for neuroblastoma intervention in patients [31]. enriched KEGG pathways, y-axis label represents pathway names, and x-axis represents enrichment factor. Size and color of the bubble represents amount of DEGs enriched in the pathway and enrichment significance, respectively. c Pathway interaction, each node is a KEGG pathway item; the node size reflects pathway significance, edge between nodes reflected common genes; and different node colors represented different functional groups Another study reported that decreased lncRNASNHG1 inhibits colorectal carcinoma tumorigenesis and might act as an important potential therapeutic target [32]. Of note, lncRNA Snhg1 is reported to be increased in hepatocellular carcinoma [33]. In the present study, we observed that the level of lncRNASnhg1was increased in NASH, whereas it was decreased by Qianggan extract. Considering the association between NASH and hepatocellular carcinoma, lncRNASnhg1 might be a therapeutic target of QGE and a potential biomarker of NASH.
According to the ceRNA hypothesis, lncRNAs act as miRNA sponges to negatively regulate miRNAs and positively regulate target mRNAs. A recent study revealed that lncRNA SNHG1 acted as a ceRNA of miR-326 targeted NOB1 [34]. A previous report showed that lncRNA SNHG1 serves as a ceRNA to negatively regulate miR-199a-3p and enhance CDK7 expression, which might be a potential therapeutic target for prostate cancer [35]. We noticed that lncRNA Snhg1 serves as a ceRNA of miR-199a-3p, and 7 other miRNAs that targeted a series of immune-related mRNAs. CD28 is a membrane protein belonging to the family of costimulatory molecules that synergistically acts as a second signal with the T-cell receptor complex [36]. Researchers have revealed that deletion of CD28 establishes a new pro/anti-inflammatory balance, and protects the liver against steatosis [37]. In accordance with previous studies, we observed that Cd28 mRNA expression was elevated in NASH and reduced after QGE treatment in mice. CD8A is a constituent of the CD8 antigen found on most cytotoxic T cells and acts as a coreceptor with the T-cell receptor [38]. Here, we noticed that Cd8a was up-regulated in NASH mice but downregulated after Qianggan treatment, suggesting Cd8a as a potential therapeutic target for NASH. CXCL9 is secreted by monocyte-derived dendritic cells within the liver and is a pro-inflammatory cytokine that attracts Th1 and Th17 cells to the liver [39]. CXCL9 was elevated in patients with autoimmune hepatitis and mice with fatty liver [40,41], indicating that CXCL9 inhibition is a potential strategy in NAFLD treatment.
MiR-324-3p is associated with malignant clinicopathologic features and poor prognosis of hepatocellular carcinoma [42]. In the present study, we noticed that miR-324-3p-targeted mRNAs were differentially expressed among the groups. IL-15 has been implicated in high-fat-diet-induced hepatic lipid accumulation and inflammation through immune regulation [43]. Researchers have analyzed the level of IL-15 in patients with acute hepatic failure and have suggested that overexpression of IL-15 might cause liver injury in humans [44]. Similar to lncRNAs, circRNAs can also act as ceRNAs to inhibit miRNAs and positively regulate associated mRNAs. We identified several cir-cRNAs acting as ceRNAs in NASH. MiR-760-3p has been implicated in immunological processes and miR-5099 in liver injury [45,46]. Intriguingly, we observed that 32 targeted mRNAs were enriched in immune related functions and several mRNAs (e.g., Ccl1 and Il18r1) were related to liver inflammation and fibrosis. Inhibition of the CC chemokine ligand CCL1 has been suggested as a successful therapeutic way to limit liver inflammation and fibrosis [47]. In addition, IL18R1 −/− mice displayed a less severe liver fibrotic phenotype than wild-type animals [48].
Klrk1, also known as Nkg2d, is an activating receptor expressed on the surface of NK cells, CD8 + T cells, and subsets of CD4 + T cells and invariant NKT cells [49]. KLRK1 mRNA increased in patients with NASH [50]. Lymphoid enhancer binding factor Lef1 is essential for the early stages of thymocyte maturation, and is involved in liver fibrosis [51,52]. The physiological significance of miR-362-5p in the regulation of NK cell function has been studied [53]. Up-regulation of miR-362-5p is associated with hepatocellular carcinoma and inhibition of miR-362-5p dramatically decreases tumor growth and metastasis [54]. Regulation of miR-362-5p could be inferred from our results. MiR-296-5p expression is reduced in liver samples from NASH patients [55]. Here, we observed that cir_0009212 and miR-296-5ptargeted mRNAs were increased in NASH mice and that QG partially restored these alterations.
Although we obtained abundant information through ceRNA networks, and hinted at the molecular mechanisms of QG on NASH, the limitations of our study should also be noted. The MCD diet-induced mice manifested obvious hepatic steatosis and inflammation, but lacked the human metabolic profile, such as insulin resistance, and dyslipidemia. The substantial weight loss of the mice was not observed in NASH patients. In addition, since we only investigated the levels of mRNAs and not their coded proteins, a thorough and comprehensive human investigation needs to be performed to confirm the results. LncRNAs associated ceRNA networks. Red rectangles represent lncRNAs, green ellipses represent mRNAs, yellow triangles representsmiRNAs, and pink nodes representoverlapped RNAs. a lncRNAsnhg1 associated ceRNA networks. b lncRNASlc36a3os-associatedceRNA networks. c Merged ceRNA networks of lncRNAsnhg1and Slc36a3os. d Enriched GO terms of mRNAs involved in lncRNAsnhg1ceRNA network. e Enriched GO terms of mRNAs involved in lncRNASlc36a3osceRNA network. f Venn diagram of overlapped mRNAs between lncRNAsnhg1and Slc36a3os associated ceRNAs. g Enriched GO terms of overlapped mRNAs Fig. 6 CircRNA-associated ceRNA networks. Red rectangles represent circRNAs, green ellipses represent mRNAs, yellow triangles represent miRNAs, and pink nodes represent overlapped RNAs. a circ_0004572-associated ceRNA networks. b Enriched GO terms of mRNAs involved in circ_0004572 networks. c circ_0009453 associated ceRNA networks. d Enriched GO terms of mRNAs involved in circ_0009453 network (e) circ_0009212-associated ceRNA networks. f Enriched GO terms of mRNAs involved in circ_0009212 network. g circ_0009029-associated ceRNA networks. h Enriched GO terms of mRNAs involved in circ_0009029 network. i Merged ceRNA networks of four circRNAs. j Venn diagram of overlapped mRNAs among four circRNAassociated ceRNAs. Blue represents circ_0009029-associated mRNAs, red represents circ_0004572-associated mRNAs, green represents circ_0009453associated mRNAs and yellow represents circ_0009212-associated mRNAs Fig. 7 Validation of representative RNAs. The x-axis represents RNA names, and the y-axis represents RNA fold changes based on average expression values. Blue bars represent data yielded by real-time qPCR, and red line with points represent data obtained by microarray. a MCD vs normal groups, (b) QG vs MCD groups