circRNomics Balance is Disrupted in HSCs of NASH
circRNAs have become a novel research hotspot in RNA biology. They localize in specific subcellular compartments, and play different biological roles[34]. However, the functions of circRNAs in NASH is still elusive. As activation of HSCs is well regarded as the central driver of fibrosis in NASH, we decided to explore the differentially expressed circRNAs in HSCs of NASH patients[8].
Thus, we downloaded the microarray gene profiling dataset (GSE134146, including 4 patients with NASH cirrhosis and 4 patients without NASH) from Gene Expression Omnibus (GEO) for bioinformatics analysis[17]. The circRNAs with P-value<0.05 and |Log2FC|≥ 1.5 from t-test were identified as differentially expressed circRNAs (Figure. 1A-1B). In total, 20 circRNAs were upregulated, while 8 circRNAs were downregulated in NASH (Figure. 1C). Among them, 12 circRNAs are less than 500nt in length (inclined to function as molecular scaffolds), 4 are between 500-1000nt, and 12 are longer than 1000nt (inclined to function as miRNA sponges) (Figure. 1D). Consistent with previous reports, most of them (20 out of 28) are encoded by exons (Figure. 1E). These results indicated that circRNomics homeostasis in HSCs is disrupted in NASH patients, and these changes in circRNA profiles may play a significant role in the occurrence and progression of NASH.
mecciRNAs Account for Half of Downregulated circRNAs in HSCs of NASH
Surprisingly, we found that half (4 out of 8) of downregulated circRNAs were encoded by the mitochondrial genome, namely hsa_circ_0089761, hsa_circ_0089762 (also known as circSCAR[17]), hsa_circ_0089763 and hsa_circ_0008882 (Figure. 2A-2B). Given the majority of circRNomics are occupied by nuclear-encoded circRNAs (99.93%), it is noteworthy that mecciRNAs accounted for 14.3% of the differentially expressed circRNAs in NASH (Figure. 2C)[35]. Among them, hsa_circ_0008882 are encoded by the heavy strand of mitochondrial genome, while hsa_circ_0089761, hsa_circ_0089762 and hsa_circ_0089763 are generated from the light strand (Figure. 2D). Additionally, their sizes span a wide range, in which lengths of hsa_circ_0089761 and hsa_circ_0089763 are over 5000 nt, while hsa_circ_0089762 and hsa_circ_0008882 are less than 400 nt (Figure. 2D). Then, we used mitochondrial-cytoplasmic separation assays to determine their cellular localization. The results demonstrated that these mecciRNAs were mainly located in mitochondria, while a part of them were also presented in cytoplasm (Figure. 2E).
mecciRNAs Regulate Fibrosis-related Signaling Pathways in HSCs
During the development of NASH, there is a constant dysfunction of mitochondria, including alterations in enzyme activities, protein expression, and signaling networks[36]. Thus, we planned to investigate the role of mecciRNAs imbalance in NASH progression.
circRNAs could function as miRNA decoys, which are defined as miRNA sponges that bind miRNAs and prevent them from suppressing their target mRNAs[37]. Noteworthy, circRNAs are generally with relatively few binding sites for miRNAs, the idea of controlling the stability and quantity of miRNAs by circRNAs and achieving measurable effects should be considered with caution[38, 39]. However, hsa_circ_0089761 and hsa_circ_0089763 are highly possible to possess powerful miRNA-sponge capacity, resulted from their long sequence lengths (over 5000 nt). Interestingly, as shown in Figure. 2D, hsa_circ_0089761's sequence completely contains hsa_circ_0089763's sequence, indicating they might have double sponge capacity for specific miRNAs.
First, utilizing the circMINE database, we predicted the miRNAs targeted by hsa_circ_0089761 and hsa_circ_0089763[40]. According to the filter criteria (a. Score cutoff>150; b. Energy cutoff<-7; c. Amount of specific miRNA-binding sites≥3), we identified 52 and 26 mecciRNA-miRNA pairs respectively (Figure. 3A). As expected, hsa_circ_0089761 and hsa_circ_0089763 shared 26 identical miRNA targets (hsa-mir-136-5p, hsa-mir-154-3p, hsa-mir-155-3p, hsa-mir-195-3p, hsa-mir-345-5p, hsa-mir-365a-3p, hsa-mir-365b-3p, hsa-mir-374a-5p, hsa-mir-384, hsa-mir-3943, hsa-mir-4705, hsa-mir-4732-3p, hsa-mir-5003-5p, hsa-mir-5193, hsa-mir-548aa, hsa-mir-548as-3p, hsa-mir-548t-3p, hsa-mir-5583-5p, hsa-mir-6511a-3p, hsa-mir-6511b-3p, hsa-mir-664a-3p, hsa-mir-6729-3p, hsa-mir-6731-3p, hsa-mir-6758-3p, hsa-mir-7844-5p, hsa-mir-487a-3p). Then, we identified 2000 mRNAs targeted by the 26 miRNAs in miRNet database (Supplementary Table. 1)[41]. It is worth mentioning that we also took the Protein-Protein Interaction (PPI) database into account when identifying targeted mRNAs. Therefore, among these targets, there may not only be changes in mRNA or protein expression levels, but also changes in protein modification status. Next, To further enhance the reliability, we only selected the part with degree frequency>1 for subsequent analysis (Supplementary Table. 2, Supplementary Figure. 1A). Finally, we successfully established a novel mecciRNA-miRNA-mRNA network in HSCs, including 2 mecciRNA nodes, 26 miRNA nodes and 1343 mRNA nodes (Supplementary Figure. 1B).
To achieve a better understanding of the functions associated with these mecciRNAs-related genes, we performed GO Term and KEGG pathway enrichment analysis. GO analysis showed an enrichment of GO terms indicative of transcription regulation, such as transcription coactivator activity, enhancer binding, and histone modification (Figure. 3B). Meanwhile, KEGG analysis results revealed the enrichment of several pro-fibrotic signaling pathways, such as Hippo signaling pathway and TGF-β signaling pathway (Figure. 3C)[10, 42]. Strikingly, the top 20 most enriched genes contained many famous NASH-related genes, such as c-MYC, SMAD2 and SMAD3 (Figure. 3D-3E) [43-46]. Taken together, these results indicated that mecciRNAs might have the potential to regulate fibrosis-related signaling pathways in HSCs.
Exosomes-mediated Crosstalk Between HSCs and Hepatic Cells in NASH
Death of hepatocytes could induce liver inflammation and then directly or indirectly promote the activation of HSCs[47]. Conversely, activated HSCs also further aggravated the damage and death of hepatocytes[48]. Hence, the crosstalk between HSCs and hepatic cells performs important functions in NASH progression, in which exosome-mediated intercellular communication cannot be ignored[49]. The exosomal cargo consists of proteins, lipids and ncRNAs, while many of the biologic effects of exosomes have been attributed to miRNAs[50]. Therefore, we were particularly interested in whether the changes in mecciRNA-miRNA profile of HSCs could influence the signaling pathways of hepatic cells.
First, using the exosomal miRNAs database in miRNet, we screened out 7 candidates among the 26 mecciRNA-targeted miRNAs (hsa-mir-485-3p, hsa-mir-618, hsa-mir-642a-5p, hsa-mir-1224-3p, hsa-mir-1248, hsa-mir-3529-3p and hsa-mir-3679-3p), which had been validated that can be secreted to extracellular microenvironment through exosomes to function (Figure. 4A)[41]. Then, utilizing the same analytical methods and filter criteria, we finally identified 419 mRNAs in hepatic cells, which might be regulated by the mecciRNA-miRNA network (Supplementary Table.3, Figure. 4B-4C).
Meanwhile, we downloaded the raw sequencing data of GSE46300 from GEO, which includes liver tissues with or without NASH, and obtained a set of differentially expressed genes between NASH and normal liver (Figure. 4D-4E) [51]. In total, 21 overlapped genes were identified (STAT3, POU5F1, TRAF3, UBE2K, CXCL5, SPAG9, MBD2, THBS1, TMOD3, NPHP1, NEDD9, COL5A3, ABCC5, FRMD6, FBXL18, TXNIP, FYTTD1, PSG4, DKK3, CYP2B6, FOXK2), in which STAT3 and THBS1 were localized at the core of protein interaction network (Figure. 4F-4H). Despite several reports about their functions in liver steatosis, in-depth molecular mechanisms about STAT3 and THBS1 regulating NASH progression still deserve further study[52-54].
A Novel Immunotyping of NASH Based on mecciRNA-related Network
Classification of specific disease based on different cellular and molecular features could deepen our understanding of this disease and even provide valuable information for treatment. Hence, in order to evaluate the effects of these 21 genes on NASH, we extracted the expression profiles of these targets from GEO46300 to construct a consensus cluster (Figure. 5A). Two subtypes (Cluster 1 and Cluster 2) were obtained with the minimum variance within the group and the maximum variance across the groups (Figure. 5B, Supplementary Figure. 2A-2D). Meanwhile, the results of principal component analysis (PCA) confirmed that there are significant differences between the two clusters (Figure. 5C).
We then performed differential gene expression analysis across the transcriptomes of these two subtypes. In total, we identified 273 differentially expressed genes (Supplementary Table. 4, Figure. 5D). GO analysis revealed that these genes were highly enriched for functions related to chemokine activity and immune cell chemotaxis (Table. 1, Figure. 5E). And KEGG analysis demonstrated that, in addition to the upregulation of immune cell chemotaxis-related pathways, IL-17 signaling pathway and TNF signaling pathway were also significantly activated in Cluster 1, which could facilitate the development of NAHS by exerting proinflammatory effects (Table. 2, Figure. 5F-5G) [55, 56].
Considering the differences between 2 subtypes were mainly focused on immune-related pathways, we further utilized ssGSEA scores to characterize the immune cell components in the 2 subtypes[33]. Despite most of the differences did not reach statistical significance, immune cell infiltration was higher in Cluster 1 than in Cluster 2 generally, which was consistent with the GO and KEGG analysis (Figure. 5H). Noteworthy, infiltration of T helper 1 (Th1) cells and NK CD56dim cells was significantly higher in Cluster 1, while central memory T (Tcm) cells infiltration was reduced (Figure. 5I). It is universally agreed that Th1 cells are proinflammatory cells and aggravate liver fibrosis development[57-59]. However, the roles of NK CD56dim cells and Tcm cells in NASH remain elusive and require further study.
In general, we established a novel immunotyping of NASH based on mecciRNA-related network, in which Cluster 1 presents increased immune cell infiltration and activation of proinflammatory signaling pathways.
Validation of mecciRNAs Networks in vitro and in vivo
First, to validate the bioinformatic conclusions, we successfully established HSC-activation model by exposing LX-2 cells (human hepatic stellate cell line) to lipopolysaccharide (LPS) at a concentration of 100 ng/ml for 48h (Figure. 6A). We confirmed that both hsa_circ_0089761 and hsa_circ_0089763 was detected in a lower amount in activated HSCs compared to wild-type HSCs (Figure. 6B). To further validate our computational predictions, we detected the expression levels of 5 core miRNAs (as test cases) in the mecciRNA-miRNA network. The results largely conformed to our predictions that miR-642a-5p, miR-1248, miR-670-3p and miR-1224-3p was upregulated in activated HSCs, while only miR-4667-3p level was comparable between the two groups (Figure. 6C). As shown in Figure. 3D, SMAD2, SMAD3 and c-MYC are core targets in HSCs mecciRNA-miRNA network, which have been reported to play vital roles in NASH progression[43-46]. The western blot results demonstrated the significant upregulation in expression levels of these proteins in activated HSCs as expected (Figure. 6D).
Meanwhile, a set of candidates were considered as the downstream effectors of mecciRNAs-miRNAs network in hepatic cells, in which THBS1 and STAT3 occupied key positions. As feeding mice an MCD diet is a widely accepted model to study relevant mechanisms in NASH, we established this animal model, and detected the expression levels of THBS1 and STAT3 in control group and MCD-diet group (Figure. 6E) [60, 61]. The results revealed that the expression of THBS1 was increased in MCD-diet group, despite the similar expression levels of STAT3 between the two groups (Figure. 6F). Considering PPI database was also included in mecciRNA network analysis, we further detected the expression level of phosphorylated STAT3. As shown in Figure. 6F, the results demonstrated the elevated level of phosphorylated STAT3 in MCD-diet group.
Based on these experimental results, we confirmed the reliability of mecciRNA-miRNA networks preliminarily. Additionally, our study also brought a distinct perspective on the understanding of relationship between mecciRNAs and NASH.