Skip to main content

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

Abstract

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.

Peer Review reports

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.

Accumulating evidence suggests that traditional Chinese medicine formulae are promising agents for treating NASH [8]. The Qianggan formula was designed according to the theory of traditional Chinese medicine, contains 16 ingredients: Artemisia scoparia Waldst. & Kitam., Isatis tinctoria L., Angelica sinensis(Oliv.)Diels., Paeonia lactiflora Pall., Salvia miltiorrhiza Bunge., Curcuma wenyujin Y.H.Chen et C.Ling, Astragalus membranaceus(Fisch.)Bunge., Codonopsis pilosula(Franch.)Nannf., Alisma orientale(Sam.)Juz., Polygonatum kingianum Collett & Hemsl., Rehmannia glutinosa (Gaertn.) DC., Dioscorea oppositifolia L., Crataegus pinnatifida Bunge., Medicated Leaven Massa Medicata Fermentata, Gentiana macrophylla Pall., and Glycyrrhiza uralensis Fisch., and has been reported as an effective and safe remedy for the treatment of NAFLD/NASH [9]. Protecting hepatocytes is proposed; however, the molecular mechanisms underlying the efficacy are still unclear.

Coding ribonucleic acids (RNAs) coupled with non-coding (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.

Methods

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.

Table 1 The composition of Qianggan formula

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 temperature-and 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-and-choline-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 plug-in 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 miRNA-mediated 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.

Table 2 Sequences of the primers used for PCR

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.

Results

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).

Fig. 1
figure 1

Phenotypic characteristics of the mice. a Body weight. b Liver weight. c Liver/body weight ratio. d Serum AST. e Serum ALT. f Serum ALP. g Serum LDH. h Liver sections stained with H&E. Data were represented as mean ± SD, n = 12 per group, a P < 0.05, b P < 0.01

Differentially expressed mRNAs, lncRNAs and circRNAs

By the cutoff of a 1.5-fold change and P < 0.05, 6193 differentially expressed mRNAs (differentially expressed genes, DEGs), 5692 differentially expressed lncRNAs (DElncRNAs) and 4843 differentially expressed circRNAs (DEcircRNAs) were identified in the MCD vs normal group. There were 514 DEGs, 1182 DElncRNAs and 443DEcicRNAs in QG vs MCD group. By Venn diagram, overlapping 244 DEGs (Fig. 2a), 259 DElncRNAs (Fig. 2b) and 98DEcircRNAs (Fig. 2c) were filtered out and chosen for further analysis. Representative overlapping DEGs, DElncRNAs and DEcircRNAs are listed in Tables 3, 4 and 5, respectively. Detailed information of differentially expressed RNAs was listed in Additional file 1: Table S1, Additional file 2: Table S2 and Additional file 3: Table S3.

Fig. 2
figure 2

Hierarchical cluster of overlapped DEGs, DElncRNAs and DEcircRNA. a Venn diagram of overlapping DEGs. b Venn diagram of overlapping DElncRNAs. c Venn diagram of overlapping DEcircRNAs. d Heatmap of overlapping DEGs. e Heatmap of overlapping DElncRNAs. f Heatmap of overlapping DEcircRNAs. Gene expression is represented by colors, with brighter red for higher values and brighter green for lower values

Table 3 Representative overlapped differentially expressed genes
Table 4 Representative differentially expressed lncRNAs
Table 5 Representative overlapped differentially expressed circRNAs

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. DEcircRNAscirc_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).

Fig. 3
figure 3

Functional enrichment of overlapping DEGs. a Top 30 enriched GO terms, y-axis represents GO terms, and x-axis represents enrichment factor. Size and color of the bubble represents amount of DEGs enriched in GO terms and enrichment significance, respectively. b Top 30 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

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.

Fig. 4
figure 4

PPI networks. a PPI network with 118 nodes. Pink nodes represent the hub nodes with more neighbors (≥5). b A significant module with 7 nodes are identified by MCODE (score = 5). c A significant module with 7 nodes are identified by MCODE (score = 4.3)

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.

LncRNA Snhg1 was a ceRNA of eight miRNAs (e.g., miR-199a-3p, miR-490-3p, and miR-412-3p) that targeted 74 mRNAs (e.g., Cd28, Cd8a, and Cxcl9) (Fig. 5a). LncRNA Slc36a3os was a ceRNA of four miRNAs (miR-324-3p, miR-680, miR-3473b and miR-3102-5p.2-5p) that targeted 66 mRNAs (e.g., Cd28, Il18rap, and Il15) (Fig. 5b). LncRNA Snhg1 and Slc36a3os shared 53 mRNAs such as Cd28, Cd8a, and Klrk1 (Fig. 5c &f). Moreover, the GO terms revealed that mRNAs involved in lncRNASngh1 or Slc36a3os ceRNA networks were enriched in immune processes (Fig. 5d, e&g).

Fig. 5
figure 5

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

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).

Fig. 6
figure 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 circRNA-associated ceRNAs. Blue represents circ_0009029-associated mRNAs, red represents circ_0004572-associated mRNAs, green represents circ_0009453-associated mRNAs and yellow represents circ_0009212-associated mRNAs

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.

Fig. 7
figure 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

Discussion

Here we confirmed the effect of QG on NASH mice, and comprehensive microarray analysis revealed a batch of differentially expressed mRNAs, lncRNAs and circRNAs. 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 high-throughput analysis could facilitate the exploration of the underlying mechanisms, and the ceRNA networks we have reported may provide valuable information.

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- cell-mediated 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]. 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 co-stimulatory 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 co-receptor with the T-cell receptor [38]. Here, we noticed that Cd8a was up-regulated in NASH mice but down-regulated 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 circRNAs 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.

Conclusions

Our comprehensive expression profiles of mRNAs, lncRNAs and circRNAs revealed that QG might ameliorate NASH by modulating lncRNA and circRNA associated ceRNA networks. This may broaden our knowledge of QG targets and have promising therapeutic implications for NASH.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Abbreviations

ALP:

Alkaline phosphatase

ALT:

Serum alanine aminotransferase

AST:

Aspartate aminotransferase

ceRNAs:

Competing endogenous RNAs

circRNAs:

Circular RNAs

DEcircRNAs:

Differentially expressed genes circRNAs

DEGs:

Differentially expressed genes

DElncRNAs:

Differentially expressed lncRNAs

GO:

Gene ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

LDH:

Lactate dehydrogenase

lncRNAs:

Long non-coding RNAs

MCD:

Methionine-and choline-deficient diet

MCODE:

Molecular Complex Detection

miRNAs:

MicroRNAs

mRNAs:

Ribonucleic acids

NAFLD:

Nonalcoholic fatty liver disease

NASH:

Non-alcoholic steatohepatitis

ncRNAs:

Non-coding RNAs

PPI:

Protein-protein interaction

QG:

Qiang-Gan extracts

References

  1. Chalasani N, Younossi Z, Lavine JE, Diehl AM, Brunt EM, Cusi K, Charlton M, Sanyal AJ. The diagnosis and management of non-alcoholic fatty liver disease: practice guideline by the American Association for the Study of Liver Diseases, American College of Gastroenterology, and the American Gastroenterological Association. Hepatology. 2012;55(6):2005–23.

    Article  Google Scholar 

  2. Baffy G, Brunt EM, Caldwell SH. Hepatocellular carcinoma in non-alcoholic fatty liver disease: an emerging menace. J Hepatol. 2012;56(6):1384–91.

    Article  Google Scholar 

  3. Karbalaei R, Piran M, Rezaei-Tavirani M, Asadzadeh-Aghdaei H, Heidari MH. A systems biology analysis protein-protein interaction of NASH and IBD based on comprehensive gene information. Gastroenterology and hepatology from bed to bench. 2017;10(3):194–201.

    PubMed  PubMed Central  Google Scholar 

  4. Takahashi Y, Sugimoto K, Inui H, Fukusato T. Current pharmacological therapies for nonalcoholic fatty liver disease/nonalcoholic steatohepatitis. World J Gastroenterol. 2015;21(13):3777–85.

    Article  CAS  Google Scholar 

  5. Day CP, James OF. Hepatic steatosis: innocent bystander or guilty party? Hepatology. 1998;27(6):1463–6.

    Article  CAS  Google Scholar 

  6. Tilg H, Moschen AR. Evolution of inflammation in nonalcoholic fatty liver disease: the multiple parallel hits hypothesis. Hepatology. 2010;52(5):1836–46.

    Article  CAS  Google Scholar 

  7. Ratziu V, Goodman Z, Sanyal A. Current efforts and trends in the treatment of NASH. J Hepatol. 2015;62(1 Suppl):S65–75.

    Article  CAS  Google Scholar 

  8. Jadeja R, Devkar RV, Nammi S. Herbal medicines for the treatment of nonalcoholic steatohepatitis: current scenario and future prospects. Evidence-based complementary and alternative medicine : eCAM. 2014;2014:648308.

    Article  Google Scholar 

  9. Li L, Zhang XJ, Lan Y, Xu L, Zhang XZ, Wang HH. Treatment of non-alcoholic fatty liver disease by Qianggan capsule. Chinese journal of integrative medicine. 2010;16(1):23–7.

    Article  Google Scholar 

  10. Leti F, Legendre C, Still CD, Chu X, Petrick A, Gerhard GS, DiStefano JK. Altered expression of MALAT1 lncRNA in nonalcoholic steatohepatitis fibrosis regulates CXCL5 in hepatic stellate cells. Transl Res. 2017;190:25–39 e21.

    Article  CAS  Google Scholar 

  11. Ragusa M, Barbagallo C, Brex D, Caponnetto A, Cirnigliaro M, Battaglia R, Barbagallo D, Di Pietro C, Purrello M. Molecular Crosstalking among noncoding RNAs: a new network layer of genome regulation in Cancer. Int J Genomics. 2017;2017:4723193.

    Article  Google Scholar 

  12. Friedman RC, Farh KK, Burge CB, Bartel DP. Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 2009;19(1):92–105.

    Article  CAS  Google Scholar 

  13. Lakshman R, Shah R, Reyes-Gordillo K, Varatharajalu R. Synergy between NAFLD and AFLD and potential biomarkers. Clin Res Hepatol Gastroenterol. 2015;39(Suppl 1):S29–34.

    Article  CAS  Google Scholar 

  14. Szabo G, Bala S. MicroRNAs in liver disease. Nat Rev Gastroenterol Hepatol. 2013;10(9):542–52.

    Article  CAS  Google Scholar 

  15. Takahashi K, Yan I, Haga H, Patel T. Long noncoding RNA in liver diseases. Hepatology. 2014;60(2):744–53.

    Article  CAS  Google Scholar 

  16. Atanasovska B, Rensen SS, van der Sijde MR, Marsman G, Kumar V, Jonkers I, Withoff S, Shiri-Sverdlov R, Greve JWM, Faber KN, et al. A liver-specific long noncoding RNA with a role in cell viability is elevated in human nonalcoholic steatohepatitis. Hepatology. 2017;66(3):794–808.

    Article  CAS  Google Scholar 

  17. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell. 2011;146(3):353–8.

    Article  CAS  Google Scholar 

  18. Cortes-Lopez M, Miura P. Emerging functions of circular RNAs. Yale J Biol Med. 2016;89(4):527–37.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Zhang S, Zhu D, Li H, Li H, Feng C, Zhang W. Characterization of circRNA-associated-ceRNA networks in a senescence-accelerated mouse prone 8 brain. Mol Ther. 2017;25(9):2053–61.

    Article  CAS  Google Scholar 

  20. Sui J, Xu SY, Han J, Yang SR, Li CY, Yin LH, Pu YP, Liang GY. Integrated analysis of competing endogenous RNA network revealing lncRNAs as potential prognostic biomarkers in human lung squamous cell carcinoma. Oncotarget. 2017;8(39):65997–6018.

    Article  Google Scholar 

  21. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, et al. Gene ontology: tool for the unification of biology. The gene ontology consortium. Nat Genet. 2000;25(1):25–9.

    Article  CAS  Google Scholar 

  22. Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44(D1):D457–62.

    Article  CAS  Google Scholar 

  23. Wei C, Luo T, Zou S, Zhou X, Shen W, Ji X, Li Q, Wu A. Differentially expressed lncRNAs and miRNAs with associated ceRNA networks in aged mice with postoperative cognitive dysfunction. Oncotarget. 2017;8(34):55901–14.

    PubMed  PubMed Central  Google Scholar 

  24. Zhu M, Wang Q, Zhou W, Liu T, Yang L, Zheng P, Zhang L, Ji G. Integrated analysis of hepatic mRNA and miRNA profiles identified molecular networks and potential biomarkers of NAFLD. Sci Rep. 2018;8(1):7628.

    Article  Google Scholar 

  25. Wang H, Zhao YL, Xu KC. Clinical and pathological study on effects of Qianggan capsule combined lamivudine on hepatic fibrosis in patients with chronic hepatitis B. Zhongguo Zhong Xi Yi Jie He Za Zhi. 2006;26(11):978–80.

    PubMed  Google Scholar 

  26. Liu L, Mao D, Lv J, Huang R. Meta analysis of efficacy and safety of Qianggan capsules on patients with non-alcoholic fatty liver disease. Chinese Traditional Patent Medicine. 2018;40(8):6.

    Google Scholar 

  27. Zheng P, Wang L, Zhang L, Liu T, Xing L, Ji G. Therapeutic effect of Qiangganjiaonang on nonalcoholic rat fatty liver by increasing the expression of liver leptin receptor and P-JAK2/P-STAT3. Chinese Journal of Integrated Traditional and Western Medicine on Digestion. 2009;17(3):141–5.

    Google Scholar 

  28. Sutti S, Bruzzi S, Albano E. The role of immune mechanisms in alcoholic and nonalcoholic steatohepatitis: a 2015 update. Expert Rev Gastroenterol Hepatol. 2016;10(2):243–53.

    Article  CAS  Google Scholar 

  29. Byun JS, Yi HS. Hepatic immune microenvironment in alcoholic and nonalcoholic liver disease. Biomed Res Int. 2017;2017:6862439.

    Article  Google Scholar 

  30. Sun C, Liu X, Yi Z, Xiao X, Yang M, Hu G, Liu H, Liao L, Huang F. Genome-wide analysis of long noncoding RNA expression profiles in patients with non-alcoholic fatty liver disease. IUBMB Life. 2015;67(11):847–52.

    Article  CAS  Google Scholar 

  31. Sahu D, Hsu CL, Lin CC, Yang TW, Hsu WM, Ho SY, Juan HF, Huang HC. Co-expression analysis identifies long noncoding RNA SNHG1 as a novel predictor for event-free survival in neuroblastoma. Oncotarget. 2016;7(36):58022–37.

    Article  Google Scholar 

  32. Sun X, Wang Z, Yuan W. Down-regulated long non-coding RNA SNHG1 inhibits tumor genesis of colorectal carcinoma. Cancer Biomark. 2017;20(1):67–73.

    Article  CAS  Google Scholar 

  33. Zhang M, Wang W, Li T, Yu X, Zhu Y, Ding F, Li D, Yang T. Long noncoding RNA SNHG1 predicts a poor prognosis and promotes hepatocellular carcinoma tumorigenesis. Biomedicine & pharmacotherapy = Biomedecine & pharmacotherapie. 2016;80:73–9.

    Article  CAS  Google Scholar 

  34. Fang YL, Chen H, Wang CL, Liang L. Pathogenesis of non-alcoholic fatty liver disease in children and adolescence: from "two hit theory" to "multiple hit model". World J Gastroenterol. 2018;24(27):2974–83.

    Article  CAS  Google Scholar 

  35. Li J, Zhang Z, Xiong L, Guo C, Jiang T, Zeng L, Li G, Wang J. SNHG1 lncRNA negatively regulates miR-199a-3p to enhance CDK7 expression and promote cell proliferation in prostate cancer. Biochem Biophys Res Commun. 2017;487(1):146–52.

    Article  CAS  Google Scholar 

  36. Papotto PH, Marengo EB, Sardinha LR, Goldberg AC, Rizzo LV. Immunotherapeutic strategies in autoimmune uveitis. Autoimmun Rev. 2014;13(9):909–16.

    Article  CAS  Google Scholar 

  37. Poggi M, Morin SO, Bastelica D, Govers R, Canault M, Bernot D, Georgelin O, Verdier M, Burcelin R, Olive D, et al. CD28 deletion improves obesity-induced liver steatosis but increases adiposity in mice. Int J Obes. 2015;39(6):977–85.

    Article  CAS  Google Scholar 

  38. Biros E, Gabel G, Moran CS, Schreurs C, Lindeman JH, Walker PJ, Nataatmadja M, West M, Holdt LM, Hinterseher I, et al. Differential gene expression in human abdominal aortic aneurysm and aortic occlusive disease. Oncotarget. 2015;6(15):12984–96.

    Article  Google Scholar 

  39. Czaja AJ. Review article: chemokines as orchestrators of autoimmune hepatitis and potential therapeutic targets. Aliment Pharmacol Ther. 2014;40(3):261–79.

    Article  CAS  Google Scholar 

  40. Saeki C, Nakano M, Takahashi H, Saito S, Homma S, Tajiri H, Zeniya M. Accumulation of functional regulatory T cells in actively inflamed liver in mouse dendritic cell-based autoimmune hepatic inflammation. Clin Immunol. 2010;135(1):156–66.

    Article  CAS  Google Scholar 

  41. Semba T, Nishimura M, Nishimura S, Ohara O, Ishige T, Ohno S, Nonaka K, Sogawa K, Satoh M, Sawai S, et al. The FLS (fatty liver Shionogi) mouse reveals local expressions of lipocalin-2, CXCL1 and CXCL9 in the liver with non-alcoholic steatohepatitis. BMC Gastroenterol. 2013;13:120.

    Article  CAS  Google Scholar 

  42. Tuo H, Wang Y, Wang L, Yao B, Li Q, Wang C, Liu Z, Han S, Yin G, Tu K, et al. MiR-324-3p promotes tumor growth through targeting DACT1 and activation of Wnt/beta-catenin pathway in hepatocellular carcinoma. Oncotarget. 2017;8(39):65687–98.

    Article  Google Scholar 

  43. Cepero-Donates Y, Lacraz G, Ghobadi F, Rakotoarivelo V, Orkhis S, Mayhue M, Chen YG, Rola-Pleszczynski M, Menendez A, Ilangumaran S, et al. Interleukin-15-mediated inflammation promotes non-alcoholic fatty liver disease. Cytokine. 2016;82:102–11.

    Article  CAS  Google Scholar 

  44. Yonekawa C, Nakae H, Zheng YJ, Wada H, Tanaka H, Tajimi K. IL-15 levels in patients with acute hepatic failure. Res Commun Mol Pathol Pharmacol. 2004;115-116:5–14.

    CAS  PubMed  Google Scholar 

  45. Pandey A, Sahu AR, Wani SA, Saxena S, Kanchan S, Sah V, Rajak KK, Khanduri A, Sahoo AP, Tiwari AK, et al. Modulation of host miRNAs transcriptome in lung and spleen of Peste des Petits ruminants virus infected sheep and goats. Front Microbiol. 2017;8:1146.

    Article  Google Scholar 

  46. Yang R, Yang F, Huang Z, Jin Y, Sheng Y, Ji L. Serum microRNA-122-3p, microRNA-194-5p and microRNA-5099 are potential toxicological biomarkers for the hepatotoxicity induced by Airpotato yam. Toxicol Lett. 2017;280:125–32.

    Article  CAS  Google Scholar 

  47. Heymann F, Hammerich L, Storch D, Bartneck M, Huss S, Russeler V, Gassler N, Lira SA, Luedde T, Trautwein C, et al. Hepatic macrophage migration and differentiation critical for liver fibrosis is mediated by the chemokine receptor C-C motif chemokine receptor 8 in mice. Hepatology. 2012;55(3):898–909.

    Article  CAS  Google Scholar 

  48. Markowitz GJ, Yang P, Fu J, Michelotti GA, Chen R, Sui J, Yang B, Qin WH, Zhang Z, Wang FS, et al. Inflammation-dependent IL18 signaling restricts hepatocellular carcinoma growth by enhancing the accumulation and activity of tumor-infiltrating lymphocytes. Cancer Res. 2016;76(8):2394–405.

    Article  CAS  Google Scholar 

  49. Lanier LL. NKG2D receptor and its ligands in host defense. Cancer immunology research. 2015;3(6):575–82.

    Article  CAS  Google Scholar 

  50. Kahraman A, Schlattjan M, Kocabayoglu P, Yildiz-Meziletoglu S, Schlensak M, Fingas CD, Wedemeyer I, Marquitan G, Gieseler RK, Baba HA, et al. Major histocompatibility complex class I-related chains a and B (MIC a/B): a novel role in nonalcoholic steatohepatitis. Hepatology. 2010;51(1):92–102.

    Article  CAS  Google Scholar 

  51. Steinke FC, Xue HH. From inception to output, Tcf1 and Lef1 safeguard development of T cells and innate immune cells. Immunol Res. 2014;59(1–3):45–55.

    Article  CAS  Google Scholar 

  52. Desert R, Mebarki S, Desille M, Sicard M, Lavergne E, Renaud S, Bergeat D, Sulpice L, Perret C, Turlin B, et al. "Fibrous nests" in human hepatocellular carcinoma express a Wnt-induced gene signature associated with poor clinical outcome. Int J Biochem Cell Biol. 2016;81(Pt A):195–207.

    Article  CAS  Google Scholar 

  53. Ni F, Guo C, Sun R, Fu B, Yang Y, Wu L, Ren S, Tian Z, Wei H. MicroRNA transcriptomes of distinct human NK cell populations identify miR-362-5p as an essential regulator of NK cell function. Sci Rep. 2015;5:9993.

    Article  CAS  Google Scholar 

  54. Ni F, Zhao H, Cui H, Wu Z, Chen L, Hu Z, Guo C, Liu Y, Chen Z, Wang X, et al. MicroRNA-362-5p promotes tumor growth and metastasis by targeting CYLD in hepatocellular carcinoma. Cancer Lett. 2015;356(2 Pt B):809–18.

    Article  CAS  Google Scholar 

  55. Cazanave SC, Mott JL, Elmi NA, Bronk SF, Masuoka HC, Charlton MR, Gores GJ. A role for miR-296 in the regulation of lipoapoptosis by targeting PUMA. J Lipid Res. 2011;52(8):1517–25.

    Article  CAS  Google Scholar 

Download references

Acknowledgments

Not applicable

Funding

This work was supported by grants from the National Natural Science Foundation of China (No. 81620108030, 81774084), the Shanghai Rising-Star Program (No.17QA1404000). The funding body has no role in the design of the study, collection, analysis, and interpretation of data, and in writing the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

MZ and ML contributed equally to this work; ML, LZ, YY and FL performed the majority of the experiments; MZ and LZ wrote the manuscript; GJ and LZ designed the project, guided the whole process. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Li Zhang or Guang Ji.

Ethics declarations

Ethics approval

All animal procedures were approved by the Animal Experiment Ethics Committee of Shanghai University of Traditional Chinese Medicine.

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Detailed information of DEGs among groups. Sheet1: The 6193 DEGs in MCD vs normal group. Sheet2: The 514 DEGs in in QG vs MCD group. Sheet3: The overlapped 244 DEGs between MCD vs normal and QG vs MCD. (XLSX 545 kb)

Additional file 2:

Detailed information of DElncRNAs among groups. Sheet1: The 5692 DElncRNAs in MCD vs normal group. Sheet2: The 1182 DElncRNAs in in QG vs MCD group. Sheet3: The overlapped 259 DElncRNAs between MCD vs normal and QG vs MCD. (XLSX 453 kb)

Additional file 3:

Detailed information of DEcircRNAs among groups. Sheet1: The 4843 DEcircRNAs in MCD vs normal group. Sheet2: The 443 DEcircRNAs in in QG vs MCD group. Sheet3: The overlapped 98 DEcircRNAs between MCD vs normal and QG vs MCD. (XLSX 349 kb)

Additional file 4:

Detailed information of GO germs and KEGG pathways. Sheet1: GO terms. Sheet2: KEGG pathways. (XLSX 64 kb)

Additional file 5:

Detailed information of predicted miRNAs for DEGs, DElncRNAs and DEcircRNAs. Sheet1: The predicted miRNAs for DEGs. Sheet2: The predicted miRNAs for DElncRNAs. Sheet3: The predicted miRNAs for DEcircRNAs. (XLSX 771 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhu, M., Li, M., Zhou, W. et al. Qianggan extract improved nonalcoholic steatohepatitis by modulating lncRNA/circRNA immune ceRNA networks. BMC Complement Altern Med 19, 156 (2019). https://doi.org/10.1186/s12906-019-2577-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12906-019-2577-6

Keywords