Adult male specific pathogen free (SPF) c57BL/6 mice (Jihui Laboratory Animal Care Co., Shanghai, China) were randomized using the lottery method into the normal control (NC) group (n = 6), which was fed a normal diet, and the NAFLD group (n = 6), which was fed an HFHC diet (2% cholesterol, 10% lard, and 88% normal diet). The sample size was estimated by the prospective difference in circRNAs between the NC and NAFLD groups using Gpower v3.1.9.7 software (Faul et al., 2009). The mice were measured for body weight weekly and raised in specific pathogen free (SPF) environment, with a 12 h light/12 h dark cycle, 24 ± 2°C ambient temperature, 50 ± 5% humidity and free access to water and pelleted chow. The HFHC diet was administrated for 16 weeks to establish a rodent model of NAFLD (Liu et al., 2020). And liver samples were collected at the end of the experiments and measured for liver weight. The liver index of the mice was calculated as follow: liver index = liver weight/bodyweight. All animal experiments were conducted under the Guideline for the Care and Use of Laboratory Animals (1996), with the approval of the ethical committee of Xinhua hospital.

Liver samples were collected from each mouse and embedded with paraffin or optimal cutting temperature (OCT) compound for pathological evaluation after the 16-weeks treatment. Frozen sections of OCT-embedded liver samples were stained with Oil Red O to assess the accumulation of neutral lipids (Mehlem et al., 2013). Hematoxylin-eosin (HE) staining of paraffin sections was used to evaluate the NAFLD-specific characteristics according to the NAS method (Kleiner et al., 2005). The steatosis was ranked by the area of hepatic steatosis according to S0 < 5%, S1>5%–33%, S2>33%–66%, S3>66%. The pathological alterations in each sample were assessed by two pathologists in a blinded manner.

Liver tissue was ground into a homogenate and centrifuged for the supernatant. The total protein concentration of the supernatant was determined using a BCA protein concentration determination kit G2026 (Servicebio, Wuhan, China). The BCA working solution was prepared by Bicinchonic acid and copper sulfate solution in 50: 1 ratio of volume. The supernatant (2 μl) was diluted in 18 μl phosphate buffer saline (PBS) buffer, mixed with 200 μl working solution, and incubated at 37°C for half an hour. After the incubation, the optical density (OD) was measured using BioTek Epoch Microplate Spectrophotometer (Agilent technologies, Santa Clara, CA, United States) at 562 nm wave length. The standard curve was calculated by the protein concentrations of 0, 25, 50, 100, 200, 300, 400, 500 μg/ml against the optical density obtained under the same condition. The concentration of total protein was calculated using the optical density according to the standard curve. TG contents in the liver of mice in both NC and NAFLD groups were measured using the TG Determination Kit (GPO-PAP method) (C061, Huili, Changchun, China) against total hepatic protein contents (Zhang et al., 2015). The 250 μl GPO-PAP working solution was mixed with 2.5 μl PBS, standard glycerol (1.7 mmol/L), or supernatant, respectively. The mixture was then subjected to automatic biochemistry analyzer (Chemray 800, Rayto, Shenzhen, China) and analyzed the OD at 300 s and 510 nm wave length. The TG concentration was determined using colorimetric method with standard glycerol after zero set with PBS.

Total RNA was extracted from the liver of mice in NC (n = 3) and NAFLD (n = 3) groups using RNAiso Plus reagent (Takara Bio, Shiga, Japan). The concentration and purity of hepatic RNA were evaluated using the NanoDrop ND-100 (Thermo Fisher Scientific, Waltham, MA, United States) and its integrity was assessed using the Agilent Bioanalyzer 2,100 (Agilent technologies). After excluding ribosomal and linear RNA, circRNA samples were subjected to cDNA synthesis, end repair, 3′ end adenylation, and adapter ligation. The cDNA library was then constructed, and its concentration and size were assessed using the Qubit^® 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, United States) and Agilent Bioanalyzer 2,100, respectively. Paired end sequencing was performed using Illumina HiSeq 2,500 (Illumina, San Diego, CA, United States) following the manufacturer’s instruction. Clean reads were further filtered using Seqtk (https://github.com/lh3/seqtk) and mapped to reference genome GRCm38 using the BWA-MEM alignment algorithm (http://bio-bwa.sourceforge.net/). Finally, the hepatic circRNA profile was obtained from clean reads using the CIRI algorithm (Gao et al., 2015) against circBase (http://circrna.org/) (Glažar et al., 2014). After normalizing data using spliced reads per billion mapping (SRPBM) (Jeck et al., 2013), differentially expressed circRNAs in the NAFLD group were obtained using the edgeR R package (Robinson et al., 2010) based on the criteria of fold-change ≥ 2 and p value < 0.05.

First, circRNA sets identified from the liver of NC and NAFLD groups were intersected to shed light on the scale and difference of their expression profile. The sequences of these circRNAs were mapped to the reference genome GRCm38 by CIRI software to uncover their chromosomal and genic distribution. The differentially expressed circRNAs were further subjected to unsupervised hierarchical clustering using TreeView and depicted using a volcano plot.

Second, the interactions between differentially expressed circRNAs and target miRNAs were investigated by analyzing the base complementation between the MREs of circRNA and the seed sequence of miRNA using the miRanda algorithm (Version:3.3a, https://omictools.com/miranda-tool) (Enright et al., 2003). The key circRNAs involved in miRNA regulation were filtered using the following parameters: score >180, energy < −30, and degree of circRNA-miRNA network ≥2. Moreover, downstream mRNAs of the circRNA-targeted miRNAs were investigated using integrated analysis of miRNA databases (miRDB (Chen and Wang, 2020), Targetscan (Agarwal et al., 2015), miRanda (Enright et al., 2003), miRWalk (Dweep and Gretz, 2015), miRMap (Vejnar and Zdobnov, 2012), miRNAMap (Hsu et al., 2008), RNA22 (Miranda et al., 2006)). Those mRNAs found in ≥6 databases, together with the miRDB score >80, were considered to be crucial in miRNA-based gene regulation (Chen and Wang, 2020). The circRNA-miRNA-mRNA regulatory network was finally established using Cytoscape v 3.8.2 (https://cytoscape.org/) (Shannon et al., 2003).

Third, target miRNA sets of circRNAs were intersected to reveal the characteristics of the circRNA-miRNA-mRNA regulatory network. To highlight the pathophysiological effects of the intersecting miRNAs, their target mRNA sets were subjected to KEGG enrichment analysis using the clusterProfiler R package (Yu et al., 2012; Kanehisa et al., 2021). Top-ranking signaling pathways with significant p-values were visualized using the ggplot2 3.3.0 R package (https://ggplot2.tidyverse.org/). The protein-protein interaction network of the key signaling pathway was constructed using the STRING database (https://cn.string-db.org/) (Szklarczyk et al., 2021) and modified using Cytoscape v3.8.2.

The circRNAs and mRNAs of each liver sample were extracted and subjected to reverse transcription with random primers and oligo dT primers using PrimeScript™ RT Reagent Kit (Takara Bio, Shiga, Japan) according to the manufacturer’s instructions, while the RT of miRNAs and U6 was performed with specific stem-loop RT primers and U6 reverse primers using the same method respectively. Quantitative PCR was performed on circRNAs, miRNAs, and mRNAs using Hieff UNICON^® qPCR SYBR^® Green Master Mix (Low Rox) (Yeasen Biotechnology, Shanghai, China) on QuantStudio 3 (Thermo Fisher Scientific, Waltham, MA, United States) with three replications. Relative expression levels were normalized to GAPDH (circRNAs and mRNAs) or U6 (miRNAs) housekeeping genes using the 2^−ΔΔct method. The divergent primers (circRNAs), convergent primers (miRNAs and mRNAs), and specific stem-loop primers used in the experiments are shown in Table 1.

Protein samples were harvested from mouse liver using RIPA lysis buffer and quantified by BCA protein concentration determination kit G2026 (Servicebio) with the same method. Total protein extracts were incubated with loading buffer at 100°C for 10 min and subjected to electrophoresis on 8% sodium dodecyl sulfate-polyacrylamide gel electrophoresis. Subsequently, proteins were transferred onto polyvinylidene difluoride membranes. The membranes were incubated with protein-free rapid blocking buffer for 20 min, and subsequently, with primary antibodies against a-tubulin (1:1,000, 66031, Proteintech, Rosemont, IL, United States), AMPK-α1 (1.5:1,000, AF6422, Affinity, Liyang, China), SREBP1 (1.5:1,000, AF6283, Affinity), and FASN (1:1,000, 10624, Proteintech) overnight at 4 °C, and with horseradish peroxidase (HRP) -conjugated secondary antibodies at room temperature for 1 h. Protein bands were visualized by chemiluminescence using AI680 (Cytiva, Marlborough, MA, United States) and assessed using ImageJ software (https://imagej.nih.gov/ij/) with integrated density/area.

The paraffin sections of liver samples were exposed to deparaffinization in xylene twice and ethanol three times sequentially for 5 min each time, rehydration with PBS, and antigen retrieval in boiled citrate buffer for 3 min. Then they were successively incubated with 3% hydrogen peroxide buffer in dark for 25 min, 3% bovine serum albumin for 30 min, and primary antibodies against AMPK-α1 (AF6422, 1:50, Affinity), SREBP1 (AF6283, 1:50, Affinity), and FASN (1:100, 10624, Proteintech) overnight at 4 °C in humidor. After reaction with HRP conjugated secondary antibody for 1 h, immunohistochemistry signals were visualized using the DAB method and assessed using ImageJ software with integrated density/area. The sections were washed in PBS 5 min for three times between every step mentioned above.

The numerical variables of this study were presented as mean ± standard error (SEM). Kolmogorov-Smirnov test and F-test were applied for assessing the normality and homogeneity of variance of the data, respectively. To investigate the difference between NC and NAFLD groups, data following normal distribution were analyzed using Student’s t-test (homogenous variance) or Welch’s t-test (inhomogeneous variance). Data following a non-normal distribution were analyzed using the Mann-Whitney test. In addition, a two-way repeated-measures ANOVA was performed to assess the effect of diet on bodyweight. SPSS version 19.0 (SPSS Inc., Chicago, IL, United States) was used for the statistical analyses. Two-side significance threshold was set at p < 0.05.

Mice that were fed an HFHC diet for 16 weeks presented a bodyweight significantly higher than that of mice in the NC group fed a normal diet (p < 0.01, Figure 1A). The liver index in the HFHC diet group was much higher than that in the NC group (p < 0.001, Figure 1B). The pathological assessment showed that all mice in the NC group were graded S0 according to the steatosis score, while those in the NAFLD group were graded ≥ S1, which also indicated that hepatic steatosis was induced by the HFHC diet (Mann Whitney test: NC vs. NAFLD, 0 vs. 2.667 ± 0.211, p = 0.0022, Figure 1C). Additionally, mice fed an HFHC diet presented a significant upregulation of hepatic TG concentration compared to mice fed a normal diet (p < 0.001, Figure 1E). Moreover, HE and Oil Red O staining highlighted widespread, hepatocyte-enriched lipid droplets in mice fed an HFHC diet (Figures 1D,F), indicating the presence of hepatocyte steatosis, and hence, the successful establishment of the animal model of NAFLD.

The sequence length, annotation, chromosomal location, and gene region distribution of the liver circRNAs in the NC and NAFLD groups were investigated. Interestingly, the sequence length of circRNAs in the NC group and NAFLD group showed normal distribution and power-law distribution, respectively (Figure 2A). The annotation ratio of circRNAs was similar in both groups (Figure 2B). There were more circRNAs encoded by mitochondrial DNA than in chromosome 1–19 and chromosome X in the NAFLD group than in the NC group (Figure 2C). In addition, the natural logarithm of circRNA count revealed that they derived from a much more reduced number of exons and introns than those of the NC group (p < 0.01–0.001, Figure 2D).

Based on the fold-change and p-value of their expressive difference, 28 annotated circRNAs were found to be differentially expressed between the NC and NAFLD groups (Figures 2E,F). These NAFLD-related circRNAs were significantly downregulated in the NAFLD compared to the NC group (Figures 2E,F). Thus, both the number and expression level of circRNAs were reduced by HFHC diet-induced rodent NAFLD.

To gain insight into the roles of NAFLD-related circRNAs, their sponging ability was assessed using the miRanda algorithm, including the criteria of the score, energy, and degree of circRNA-miRNA network. Databases screening and miRDB scoring further revealed crucial mRNAs downstream of circRNA-regulated miRNAs. The circRNA-miRNA and miRNA-mRNA interactions formed a circRNA-miRNA-mRNA regulatory network comprised by five circRNAs (circ_0001452, circ_0001453, circ_0001454, circ_0000503, circ_0001499), 10 miRNAs (miR-1943-5p, miR-466m-3p, miR-3470b, miR-669c-3p, miR-706, miR-466d-5p, miR-466i-5p, miR-466i-3p, miR-1187, miR-466k), and 282 mRNAs (Figure 3).

Most of the key circRNAs (4/5) of the circRNA-miRNA-mRNA regulatory network shared the target miRNAs miR-466i-3p and miR-669c-3p (Figure 4A). circ_0001452, circ_0001453, and circ_0001454 were derived from LNCPINT (Gene ID: 232685). Additionally, miR-466i-3p and miR-669c-3p shared 77.84% (74/95) and 57.81% (74/128) of the target mRNAs, respectively (Figure 4B). Members of the AMPK signaling pathway (such as AMPK-α1, Pfkfb2, and Irs4) featured these miRNA targets (Figure 4B). KEGG pathway analysis further revealed the signaling pathways with top-ranking fold enrichment and p < 0.05, which included signaling pathways associated with the glycolipid metabolism, such as the AMPK, insulin, mTOR, and insulin resistance signaling pathways (Figure 4C). The gene encoding AMPK-α1 was the target of both miR-466i-3p and miR-669c-3p. Additionally, AMPK-α1 was a central hub in the protein-protein interaction network of the AMPK signaling pathway (Figure 4D). Therefore, LNCPINT-derived circRNAs were suggested to play an important role in NAFLD through miR-466i-3p- and miR-669c-3p-dependent regulation of the AMPK signaling pathway.

In the NAFLD group, the hepatic levels of LNCPINT-derived circRNAs (circ_0001452, circ_0001453, circ_0001454) and circ_0000503 were lower compared to those in the NC group (Figures 5A,C,E,G). Deficient circRNAs generated by back-spliced junctions contained multiple MREs specific to miR-466i-3p and miR-669c-3p (Figures 5B,D,F,H). The complementary binding between the MREs of deficient circRNAs and the seed sequences of the corresponding miRNAs (Figures 5I,K) was abolished, leading to the reactivation of miR-466i-3p and miR-669c-3p. The hepatic expression level of miR-466i-3p and miR-669c-3p was also significantly upregulated in the NAFLD group compared to that in the NC group (Figures 5J,L).

Both miR-466i-3p and miR-669c-3p were recognized to repressively interact with AMPK-α1, one of the most important members of the AMPK family, through the seed sequence (located in the miRNAs) and the 3’ UTR (in the mRNA) (Figures 6A,B). The loss of LNCPINT-derived circRNAs in the NAFLD group lead to an increased level and functional activation of these miRNAs, and subsequent downregulation of AMPK-α1 expression at the transcriptional and translational levels (Figures 6C,F,G) in hepatocytes (Figure 6J). AMPK-α1 has been shown to inhibit the expression of downstream lipogenic genes, such as SREBP1 and FASN. Hence, as expected, miR-466i-3p- and miR-669c-3p-mediated AMPK-α1 downregulation markedly increased the mRNA and protein levels of both SREBP1 and FASN (Figures 6D–F,H,I). Immunohistochemical analysis revealed that hepatocytes with high-level of SREBP1 and FASN were more abundant in the NAFLD group than in the NC group (Figure 6J). Semi-quantitative analysis directly exhibited the expressive difference of AMPK-α1, SREBP1 and FASN between NC and NAFLD groups (Supplementary Figure S1). Finally, hepatic steatosis, the most characteristic pathological feature of NAFLD, also occurred in these hepatocytes (Figure 6J).

Consistent with their putative roles in lipogenesis, the expression level of circ_0001452 (Figure 7A, r = −0.671, p = 0.017), circ_0001453 (Figure 7B, r = −0.692, p = 0.013), circ_0001454 (Figure 7C, r = −0.686, p = 0.014), and circ_0000503 (Figure 7D, r = −0.686, p = 0.014) was negatively associated with the TG concentration in the liver of rodents. Mice in the NAFLD group presented lower levels of LNCPINT-derived circRNAs and higher hepatic TG content compared with those of the NC group. The loss of LNCPINT-derived circRNAs lead to miR-466i-3p and miR-669c-3p activation, which was remarkably associated with the significant upregulation of hepatic TG levels (Figures 7E,F). In similar, both miR-466i-3p and miR-669c-3p levels were negatively correlated with the expression levels of LNCPINT-derived circRNAs, respectively (Supplementary Figures S2A,B). But the correlation was not observed between hepatic levels of both miRNAs and circ_0000503 (Supplementary Figures S2A,B).