This article is mentioned in:

Introduction

Excessive fat deposition in hepatocytes is a major characteristic of NAFLD, and its pathological types generally include four stages of simple steatosis, non-alcoholic steatohepatitis (NASH), hepatic fibrosis and hepatocellular carcinoma (1). The worldwide prevalence of NAFLD is 20-30%, occurring particularly in wealthy or developed countries (2), and it is estimated that NAFLD prevalence is 30-40% in the US, with 3-12% of adults having NASH (3). Previous studies have shown that NASH is a key stage in the progression of NAFLD to liver cirrhosis, hepatocellular carcinoma and other liver diseases (4-7). Moreover, the ‘two-hit’ theory is the most widely accepted theory for the pathogenesis of NAFLD. The hypothesis states that simple hepatocyte steatosis caused by the deposition of triglycerides and fatty acids in liver cells represents the ‘first hit’; and increased levels of oxidative stress, insulin resistance and inflammatory cytokines induced by fat deposition in the liver represents the ‘second hit’ (8). However, the ‘multiple hits’ hypothesis has also been proposed, in which it is speculated that the progression of NAFLD is driven by multiple factors, including insulin resistance (IR), hormones, intestinal microbiota and genetic susceptibility (9). However, there is currently no effective treatment for NAFLD (10).

P-element-induced wimpy testis (PIWI)-interacting RNAs (piRNAs) were first identified in 2006 in the ovarian germline cells of Drosophila. These RNAs are short non-coding RNAs measuring 24-32 nucleotides and perform their biological functions by interacting with Piwi protein to form piRNA-induced silencing complex (piRISC) complexes (11-13). It has been shown that piRNAs can have biological roles in transposon silencing, gene regulation, protein expression, genome rearrangement and reproductive stem-cell maintenance (12). Previous studies have also reported that piRNAs are differentially expressed in various tumors, rheumatoid arthritis, human ageing, neuronal axon regeneration and chemotherapy resistance (14-18). Furthermore, as key proteins in assisting the biological role of piRNAs, the Piwi protein associated with piRNAs shows the same characteristics. For example, the expression of Piwi-like protein 1 (PiwiL1) in tumors is significantly correlated with histological tumor grade, clinical stage and poor clinical outcome (11).

However, to the best of our knowledge, no direct correlation between NAFLD and piRNA has been previously reported. Therefore, in the present study, piRNA expression in NAFLD was analyzed using a piRNA gene chip and subsequent bioinformatics techniques. Thus, the present results may provide a basis for future research on the role piRNAs in the pathogenesis of NAFLD.

Materials and methods

Animals and materials

C57BL/6 mice were purchased from Chongqing TengXing Biotechnology Co. Ltd. (n=20; weight, 20 g; age, 8 weeks). Methionine- and choline-deficient (MCD) diet and methionine- and choline-sufficient (MCS) diet model feeds were purchased from the Trophic Animal Feed High-tech Co., Ltd. Alanine aminotransferase (ALT; cat. no. C009-2), aspartate aminotransferase (AST; cat. no. C010-2), total cholesterol (TC; cat. no. A111-1) and triglyceride (TG) assay kits (cat. no. A110-1) were purchased from Nanjing Jian Cheng Bioengineering Institute. Arraystar MM9 piRNA array and microarray experiments were performed by Kang Chen Bio-Tech, Inc. TRIpure® Total RNA extraction reagent, EntiLink™ 1st Strand cDNA synthesis kit and EnTurbo™ SYBR Green PCR SuperMix were purchased from ELK (Wuhan) Biotechnology Co., Ltd.

Mouse NAFLD model

In total, 20 C57BL/6 male mice were randomly divided into two groups. Then, 10 mice in the control group were fed the MCS diet, while the other 10 mice in the model group were fed the MCD diet. All mice were fed for 8 weeks in a clean room with free access to water and food (temperature, 20±5˚C; humidity, 50%; 12/12 h light/dark cycle). The mice were weighed regularly and the liver index was calculated using the following formula: (Liver wet weight/Body weight of mice) x100%. In order to alleviate pain at the time of sacrifice, the mice were first weighed and then given deep anesthesia by intraperitoneal injection of 10% chloral hydrate (dose, 500 mg/kg). Mice were judged to be in a state of deep anesthesia after complete muscle relaxation, disappearance of corneal reflex and no response to external stimuli. Then, the mice were sacrificed by cervical dislocation, and mortality was assessed by the inability to observe respiratory movements in the thorax and the inability to feel the heartbeat.

Liver histopathological examination

After the mice were sacrificed, a longitudinal incision was made in the middle of the abdomen to expose the abdominal cavity, and the abdominal skin, muscles and fascia were gradually separated. The liver tissue was fully exposed and isolated, and the wet weight of the liver tissue was measured. Fresh liver tissue was placed into a fixed solution for 24 h at 4˚C (4% neutral formaldehyde) followed by alcohol dehydration and xylene replacement of the alcohol in the tissue. Then, paraffin embedding was performed and liver tissue was cut into 5-µm sections. These sections were then incubated at 45˚C for 2 h, and stained with 0.4% hematoxylin for 2-3 min at 25˚C and 0.5% eosin (H&E) for 1 min at 25˚C. Subsequently, the liver tissue was assessed by NAFLD activity score (NAS) (19). The NAS (0-8 points) was assessed by i) hepatocyte steatosis: 0 points (<5%); 1 point, 5-33%); 2 points, 34-66%; 3 points >66%; ii) inflammation in the hepatic lobule (count necrotic foci at x20 magnification): 0 points, none; 1 point, <2; 2 points, 2-4; 3 points, >4; and iii) hepatocyte ballooning: 0 points, none; 1 point, rare; 2 points, many. NASH was excluded if the NAS was <3, and NASH was diagnosed if the NAS was >4.

Determination of serum biochemical indexes

After the mice were sacrificed, blood samples (2 ml) of the mice were obtained by orbital collection. Subsequently, the blood samples were centrifuged at 110 x g for 10 min and 1 ml plasma was collected. Serum levels of ALT, AST, TG and TC were determ ined by using AST, ALT, TC and TG assay kits (Nanjing Jian Cheng Bioengineering Institute), according to the manufacturer's instructions.

Quality control of total RNA in the liver tissues

Total RNA was extracted with TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. RNA integrity was evaluated by electrophoresis on 1% denatured agarose gels. Total RNA of eukaryotic samples run on the denaturing gel had distinct bands of 28S and 18S rRNA. If the 28S rRNA band was ~2 times as strong as the 18S rRNA band, then the RNA of this sample could be considered complete (20). RNA concentration [optical density (OD) 260], organic compound contamination (ratio OD260/OD230) and protein contamination (OD260 /OD280 ratio) were measured using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc.). Total RNA with an OD260/OD280 ratio >1.8 was then selected (20).

Microarray analysis of piRNAs

In total, three liver tissues from control and model groups were randomly selected for detection by Arraystar MM9 piRNA array. Then, 1 µg RNA was extracted from each sample, and cyanine 3 (Cy3) was fluorescently labelled at the 3' end using T4 RNA ligase (cat. no. EL0021; Thermo Fisher Scientific, Inc.) via a RNA ligase method. Microarray analysis was performed using the Agilent Array platform (Agilent Technologies, Inc.). RNA labelled with Cy3 was hybridized with Arraystar piRNA Array in Agilent's SureHyb hybridization chambers for 17 h at 65˚C. Then, an Agilent DNA Microarray Scanner (Agilent Technologies, Inc) was used for the scanning of the arrays. The array images were analyzed by Agilent Feature Extraction software (version 11.0.1.1; Agilent Technologies, Inc). The GeneSpring GX v11.5.1 software package (Agilent Technologies, Inc.) was used to normalize and process data for quantiles. Filtering fold change and volcano plots was performed to identify piRNAs with significant differential expression between the two groups. Moreover, R software (version 3.1.2; https://www.r-project.org/) was used to create hierarchical clustering.

Bioinformatics analysis

miRanda (http://www.microrna.org/microrna/home.do) was used to search for potential target genes of the differentially expressed piRNAs (21,22). The function and signaling pathways of target genes were predicted by Gene Ontology (GO; http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses (KEGG; http://www.genome.jp/kegg/ko.html).

Detection of piRNA expression by reverse transcription-quantitative PCR (RT-qPCR)

Liver samples were stored at -80˚C, and ~100 mg of tissue was collected from the ice with sterilized tools and ground in 1 ml pre-cooled TRIpure®. The homogenate was poured into a 1.5-ml EP tube, and 250 µl trichloromethane was added, mixed and left on ice for 5 min. The mixture was centrifuged at 10,000 x g at 4˚C for 10 min. Then, 500 µl supernatant was absorbed into a 1.5-ml EP tube in an ultra-clean working platform, and an equal volume isopropanol pre-cooled at 4˚C was added, mixed upside down and incubated at -20˚C for 15 min. The solution was centrifuged at 4˚C and 10,000 x g for 10 min, the liquid was carefully poured out, 1 ml 4˚C pre-cooled 75% ethanol was added and inverted several times. The RNA precipitate was washed using 75% ethanol, centrifugated at 4˚C and 10,000 x g for 5 min, the liquid was poured out and RNA was dried for 5 min at 4˚C on a clean work table. Then, the RNA was fully volatilized with ethanol and 10 µl RNase-Free Water was added to fully dissolve the RNA.

Subsequently, an EntiLink™ 1st Strand cDNA Synthesis kit [ELK (Wuhan) Biotechnology Co., Ltd.; cat. no. EQ003] was used to perform first strand cDNA. The mixture of 1.0 µl internal reference gene-specific reverse transcription primer, 1.0 µl target gene-specific reverse transcription primer, 1.0 µl dNTPs and 15.0 µl RNA was placed on a PCR instrument at 70˚C for 5 min, and then cooled quickly on ice for 2 min. Then, the reaction solution was placed on ice, followed by the addition of 4.0 µl 5X RT buffer, 1.0 µl M-MLV reverse transcriptase and 1.0 µl RNase inhibitor, and the mixture was placed on a PCR instrument at 42˚C for 60 min and 85˚C for 5 min.

RT-qPCR was performed using the StepOne™ RT PCR instrument (Thermo Fisher Scientific, Inc.) using the Turbo™ SYBR Green PCR SuperMix kit [ELK (Wuhan) Biotechnology Co., Ltd.; cat. no. EQ001]. The reaction system was as follows: 1.0 µl primer working solution, 5.0 µl 2X Master Mix, 1.0 µl template, 2.0 µl ddH2O and 1.0 µl Rox. The following thermocycling conditions were used: Initial denaturation at 95˚C for 3 min, followed by 40 cycles at 95˚C for 10 sec, 58˚C for 30 sec and 72˚C for 30 sec. Then, the dissolution curve was drawn according to the default settings of the instrument. Relative gene expression was calculated using the 2-ΔΔCq method (23). All experiments were repeated three times with U6 as the internal reference gene. The primers used in this assay are shown in Table I.

Table I Oligonucleotide sequences of the primers.

Table I

Oligonucleotide sequences of the primers.

Primer	Primer location	Sequence (5'-3')
U6	RT	AACGCTTCACGAATTTGCGT
	Forward	CTCGCTTCGGCAGCACAT
	Reverse	AACGCTTCACGAATTTGCGT
DQ566704	RT	CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGCTCCCGCC
	Forward	CTTAGGACAGTGCCGAGGG
	Reverse	CTCAACTGGTGTCGTGGAGTC
DQ723301	RT	CTCAACTGGTGTCGTGGAGTCGGCAATTCAGTTGAGGGGCGTCT
	Forward	GTGGCATTTAAAGATACCGGA
	Reverse	CTCAACTGGTGTCGTGGAGTC

[i] RT, reverse transcription.

Statistical analysis

SPSS for Windows version 19.0 was used for statistical analysis (SPSS, Inc.). Pairwise comparisons were made using Student's t-test. Data are presented as the mean ± SD. P<0.05 was considered to indicate a statistically significant difference.

Results

Body weight, liver wet weight and liver index of each group of mice

It was found that the post-modeling body weight and liver wet weight of mice in the model group were significantly lower compared with the control group. However, the liver index of the model group was significantly increased in the model group (Table II).

Table II Statistical analysis of liver index and other indexes.

Table II

Statistical analysis of liver index and other indexes.

Group	Weight pre-modeling, g	Weight post-modeling, g	Liver wet weight, g	Liver index, %
Control	20.9300±0.5696	25.4200±1.0336	1.7500±0.1957	0.0687±0.0061
Model	20.8900±0.4817	16.6600±0.6040a	1.4800±0.1032a	0.088±0.0042a

[i] Data are presented as the mean ± SD.

[ii] ^aP<0.05 vs. control group.

Serum biochemical index

Serum levels of AST and ALT in mice of the model group were significantly higher compared with the control group (Fig. 1A and B); however, the levels of TG and TC were significantly decreased compared with the control group (Fig. 1C and D).

Figure 1 Serum biochemical index of mice. (A) Level of AST. (B) Level of ALT. (C) Level of TG. (D) Level of TC. *P<0.05, ***P<0.001. ALT, alanine aminotransferase; AST, aspartate aminotransferase; TC, total cholesterol; TG, triglyceride.

H&E staining of liver tissue and pathological assessment

The pathological results identified no obvious abnormalities in the liver tissue of the control group; the hepatic lobules were orderly and clear, and the hepatic cells were arranged radially around the central vein. Furthermore, in the control group the hepatocytes were normal in size with uniform cytoplasm, there was no lipid droplet deposition in liver cells and the morphology structure was normal (Fig. 2A and B). However, the liver tissue of mice in the model group was highly different from the healthy liver tissue, as the hepatic lobule and hepatic cord were disorganized in the model group. Moreover, liver cells were swollen and poorly demarcated. In addition, with varying degrees of steatosis and balloon degeneration, numerous hepatocytes had punctate or focal necrosis, and inflammatory cells were found around the lobules and portal veins. The results also identified lipid droplets of different sizes in cells, and nuclei were extruded to the edges of the cytoplasm (Fig. 2C and D). Moreover, the NAS of the liver tissue of the control group was 0-1, thus NASH diagnosis could be excluded. However, the NAS of the model group was 6-7 points, which indicated a potential NASH diagnosis.

Figure 2 Hematoxylin and eosin staining of liver tissues obtained from mice in the (A) control group at magnification, x100 and (B) control group at magnification, x200. In the control group, hepatic lobules structures were regular and clearly visible, and hepatocytes were radially arranged around the central vein. The hepatocytes were normal in size and there was no lipid droplet deposition in liver cells. Hematoxylin and eosin staining of liver tissues obtained from mice in the (C) model group at magnification, x100 and (D) model group at magnification, x200. The hepatic lobular structures and the hepatic cord structure in the model group were disordered, and liver cells were swollen and poorly demarcated. There was varying degrees of fatty degeneration and ballooning degeneration (blue arrow), numerous liver cells that demonstrated punctate or focal necrosis (red arrow) and there were also numerous inflammatory cells (black arrow).

piRNA microarray data analysis

A box plot was used to visualize the dataset distributions and show the intensity of normalized gene data. In the piRNA array experiment, the absolute optical density of each array is different. It is necessary to normalize the experimental data before comparing them. The purpose of normalizing the gene data is to eliminate the gene intensity errors caused by experimental techniques, rather than to adjust the differences in biological RNA samples, so that each experimental sample and parallel experimental data are at the same level, so as to obtain the gene intensity with biological significance. The results suggested that the distributions of log2 ratios among the three paired samples were highly similar, indicating that the normalization results of gene expression data are well and can be used for further comparative analysis (Fig. 3A). A hierarchical clustering was given according to the ‘All Targets Value piRNAs’ (Fig. 3B), and the result of hierarchical clustering demonstrated distinguishable gene expression profiling among the samples. A total of 1,285 differentially expressed piRNAs were identified in the model group, of which 641 were upregulated and 644 were downregulated. Furthermore, a scatter plot was used to assess the variation between gene chips (Fig. 3C).

Figure 3 piRNA microarray data analysis. Control 1, Control 2, Control 3: Three liver tissues from control group. Model 1, Model 2, Model 3: Three liver tissues from model group. (A) Box plot was used to compare the intensity distribution from the three paired samples. (B) Hierarchical clustering was performed based on ‘All Targets Value-piRNAs’. The result of Hierarchical Clustering identified a distinguishable piRNA expression profile among samples. Red indicates high relative expression, and green indicates low relative expression. (C) A scatter-plot was used for assessing piRNA expression variation between arrays. The values on the x- and y-axes in the scatter-plot are the normalized signal values of the samples (log2 scaled). Black represents fold change lines (the default fold change value given was 2.0). The piRNAs above the top black part and below the bottom black part indicated >2.0 fold change of piRNAs between the two samples. piRNA, PIWI-interacting RNA.

A volcano plot was performed to conveniently visualize the differentially expressed piRNAs between the two groups. (Fig. 4). According to the P-value, fold change and raw intensity of samples in the chip analysis of differentially expressed piRNAs, 10 piRNAs with potential research value were found in the upregulated group and downregulated group, respectively (Table III).

Figure 4 A volcano plot was used to visualize the differentially expressed piRNAs. The vertical lines correspond to 2.0-fold up and down, respectively, and the horizontal line represents a P-value of 0.05. Therefore, the red points and green points in the plot represent the differentially expressed piRNAs with statistical significance; piRNA, PIWI-interacting RNA.

Table III Fold change and length of the 10 upregulated piRNAs and 10 downregulated piRNAs of the piRNA microarray data.

Table III

Fold change and length of the 10 upregulated piRNAs and 10 downregulated piRNAs of the piRNA microarray data.

Regulation	piRNA	P-value	Fold change	Sequence (5'-3')
Up	DQ606826	0.000674055	12.6207222	TGAACCAATGAGTTTATTGGGGTTACTTAG
	DQ710371	0.008904193	11.2274599	TAGAACAGAAAGTATATGGTAATGGAGAC
	Muu12120	0.011419083	8.7360271	GTTGAAACAGTAAAGGAGACCACCAG
	DQ545927	0.010209107	8.1960626	TGAATACGATAAGGAAACTCTGTGGAGAC
	DQ714914	0.002075764	7.0829049	TGAGAGAATAATGAAACCAAAAGGAGACTC
	DQ697977	0.037974174	6.3357343	TAATAAAAAGAAGGAACTCAGGGGGAGA
	DQ566704	0.001872297	4.8126516	TAGGACAGTGCCGAGGGCAATGGCGGGAG
	DQ692555	0.001940126	4.1488898	TACGTAGAATATTGAATGTAGGTGGAGAC
	DQ702021	0.013655089	3.6915243	TAAAAGATTAAAACTAGGGAGAGCTAAGA
	DQ723301	0.004944453	3.428635	TGGCATTTAAAGATACCGGAAAGACGCCC
Down	DQ723951	0.000805666	8.1575768	TCAGTTTTTTTTTGGTGGCCCCTCCCCC
	DQ705842	0.000190539	6.6600453	TGTTGTTTTTTTTCTTGAGTCTTCCCT
	DQ692800	0.020804572	6.0003043	TGTGGGAGAAGTGACCGCAGTCTTCTGCC
	DQ555235	0.027459524	5.9141076	TGGGTGACTTGTGTTGTAGAATTATGTA
	DQ706232	0.040900385	5.782601	TGATTATTCCTGGGGGAATTCCAATTTC
	DQ707539	0.045159786	5.6721689	TGTCCTAAGAACAAAAGCAAAACTCAAGGAA
	DQ719597	0.005348353	5.11935298	GGTCGATGATGAGAGCTTTGTTCTGAGC
	DQ708554	0.03392006	5.1573021	TGGCACGATGATCTCTGCGGCCGGC
	DQ553898	0.008277426	5.1545273	TGGCATGGCTCTGAGTGGTATATGTGGTC
	DQ701588	0.21921817	5.1502508	TGGAGTCCAGACTGGTTGGACTGGGTC

[i] piRNA, PIWI-interacting RNA.

Bioinformatics assessments

It has previously been shown that piRNAs regulate protein-coding genes (24). It was speculated that the differentially expressed piRNAs in NAFLD may also affect the expression of downstream genes. The bioinformatics assessment results demonstrated that the differentially expressed piRNAs had numerous target genes, and GO analysis was then carried out on the predicted target gene. The results indicated that the predicted target genes were associated with numerous biological processes, most of which were enriched in ‘cellular metabolic processes’ and ‘metabolic processes’ (Fig. 5A). Furthermore, the target genes were found to be associated with several cellular components (Fig. 5B).

Figure 5 Bioinformatics assessments. (A) GO analysis was performed to investigate biological processes associated with target genes. Target genes were found to be associated with numerous BPs, most of which were concentrated in the process of cellular metabolism and metabolism. (B) GO analysis was performed to investigate the CCs associated with target genes. Target genes were identified to be associated with numerous CCs. GO, Gene Ontology; BP, biological process; CC, cellular components.

The results also suggested that the majority of target genes were associated with ‘protein binding’ in the domain of molecular function (Fig. 6A). It was found that several major signaling pathways were associated with differentially expressed piRNAs, such as ‘Pathways of cancer’, the ‘Hippo signaling pathway’, the ‘Wnt signaling pathway’, ‘Cushing syndrome’, the ‘Forkhead box protein O (FoxO) signaling pathway’ and the ‘Mitogen-activated protein kinase (MAPK) signaling pathway’ (Fig. 6B). Moreover, several of these pathways have been shown to play a role in the development and progression of NAFLD (25-27). Therefore, the results suggest that the identified piRNAs may represent a novel class of regulators in NAFLD.

Figure 6 Bioinformatics assessments. (A) GO analysis was performed to investigate the molecular functions of target genes. The majority of target genes were indicated to be associated with protein binding. (B) Kyoto Encyclopedia of Genes and Genomes pathway cluster analysis of target genes. Differentially expressed piRNAs were found to be associated with ‘Cancer’, ‘Hippo signaling pathway’, ‘Wnt signaling pathway’, ‘Cushing syndrome’, ‘FoxO signaling pathway’ and ‘MAPK signaling pathway’. GO, Gene Ontology; MAPK, mitogen-activated protein kinase; FoxO, Forkhead box protein O; MF, molecular functions; piRNA, P-element-induced wimpy testis interacting RNA.

Validation of the expression levels of piRNAs by RT-qPCR

The expression levels of piRNAs (piRs) in the control group and model group were examined by RT-qPCR. The results demonstrated that piR-DQ566704 and piR-DQ723301 were significantly upregulated in the model group compared with the control group (Fig. 7). Furthermore, the results of RT-qPCR were mostly consistent with those of gene chip analyses, which suggested that the results of gene chip analyses were reliable.

Figure 7 Validation of the expression levels of piRNAs by RT-qPCR. RT-qPCR was performed to investigate the expression levels of piR-DQ566704 and piR-DQ723301 in the model and control groups. Quantification of expression levels was determined using the 2-ΔΔCq method. The results demonstrated that the expression levels of piR-DQ566704 and piR-DQ723301 were increased in the model group compared with the control group, which were largely consistent with the results obtained from piRNA gene chip analyses. Data are from three independent experiments and are presented as the mean ± SEM. *P<0.05 vs. control group. RT-qPCR, reverse transcription-quantitative PCR; piRNA, PIWI-interacting RNA.

Discussion

Currently, there are two main methods to construct a mouse non-alcoholic fatty liver disease (NAFLD) model; the first method is to establish a mouse NAFLD model is via genetic modification, and the second type is by dietary intervention. In the present study, the mouse NAFLD model was constructed by adjusting diet, which was hypothesized to be more similar to the natural course of patients with NAFLD and minimizes the influence on normal gene expression.

Common NAFLD model diets include: i) MCD (28); ii) choline deficiency and L-amino acid-defined diet (29); iii) cholesterol and cholate diet (30); iv) high-fat diet (31); v) high-fructose diet (32); vi) high-fat and gold-glucose-glucose diet (33); and vii) high-fructose high-fat diet (34). Moreover, the various model diets have certain differences in the modeling of mouse NAFLD, including different modeling time, pathological effects of models, degrees of inflammatory fibrosis, metabolic patterns of models and experimental repeatability (28-34). A methionine- and choline-deficient (MCD) diet-induced NAFLD is a commonly used animal model for studying NAFLD and its associated inflammation and fibrosis (28,35). Furthermore, the MCD diet can induce hepatocyte steatosis and inflammation within a few weeks and even lead to liver fibrosis in severe cases, providing an effective mouse pathological model to study NAFLD and its complications (28,35). In addition, C57BL/6 mice are a strain of mice that are sensitive to metabolic syndrome (male C57BL/6 mice are more sensitive) and are widely used to replicate metabolic syndrome models, including IR, hyperlipidemia, diabetes, atherosclerosis and other metabolic diseases (36-38). Therefore, the present study chose a MCD diet combined with C57BL/6 male mice to construct the NAFLD mouse model. However, the single NAFLD modeling method and the single mouse sex may lead to less comprehensive results in the present study. Therefore, the differentially expressed piRNA screened in the present study require examination in a variety of NAFLD mouse models. Moreover, the expression levels of differentially expressed piRNA require investigation in NAFLD models constructed in female mice.

Identified as nuclear proteins, there are four human Piwis: PiwiL1, PiwiL2, PiwiL3 and PiwiL4(39). These proteins form the piRISC complex with piRNA to regulate germline transposable elements and the protein-coding gene transposable elements in soma (13,40). Inflammation is the key point of the progression of NAFLD, and it has gained increasing attention. Previous studies have reported high expression of inflammatory factors in the progression of NAFLD, such as NF-κB, interleukin (IL)-6 and tumor necrosis factor (TNF)-α, and their expression levels vary with the severity of inflammation (41,42). Piwi proteins, as a key proteins associated with piRNA, have been shown to be correlated with inflammatory cytokines. For example, PiwiL2 protein, an incomplete Piwi protein, is associated with the upregulation of the Bcl-2 gene and the activation of NF-κB (43). The synovial fibroblasts of rheumatoid and osteoarthritis, with the stimulation of TNFα + IL1β/Toll-like receptor-ligands, the mRNA expression levels of PiwiL2 and PiwiL4 are significantly increased, and piRNA-16659 is also significantly induced by the stimulation of Poly (I:C) (14). Therefore, results suggest that piRNAs may also be associated with NAFLD inflammation.

In addition, in the present study, the related pathways of differential expression piRNA in NAFLD were predicted by bioinformatics; these pathways included Hippo, Wnt, FoxO and MAPK. Furthermore, previous studies have also revealed that these signaling pathways play an important role in liver fibrosis (44-46), metabolic syndrome and lipid metabolism (47-50).

There are four subtypes of MAPK: ERK1/2, p38 MAPK, JNK and ERK5(51). Previous studies have shown that ERK1/2 can stimulate the expression of fatty acid synthase by regulating the nuclear maturation of the Sterol regulatory element-binding transcription factor 1 C (SREBP-1C), and an increase in ERK1/2 phosphorylation may increase the expression of SREBP-1C (47). Moreover, SREBP-1C can regulate the synthesis and storage of triglycerides and phospholipids, and its excessive expression can increase the synthesis of fatty acid synthase (48). Fatty acid synthase is a key enzyme in the synthesis of fatty acid, which can inhibit the phosphorylation of insulin receptor (IRS-1), downregulate the transcription level of IRS-2 and inhibit the normal transduction of insulin signal, leading to IR (49), thus forming NAFLD. The P38 MAPK signaling pathway is primarily involved in IR induction of abnormal expression of intrahepatic lipid regulatory protein, leading to NAFLD (52). When the phosphorylation level of p38 MAPK increases in vivo, the phosphorylation level of IRS-1/2 decreases and leads to IR, which promotes the level of free fatty acid (FFA) in blood. In addition, it has been shown that FFA can induce the upregulation of aquaporin 9 and increase the uptake and utilization of glycerol in hepatocytes (50).

Activation of hepatic stellate cells (HSCs) is a key link in liver fibrosis (53). A DNA microarray was used in a previous study to detect the differentially expressed genes in quiescent and activated rat HSCs, which revealed that the ligands Wnt4 and Wnt5 and the Wnt receptor Frizzled-2 are upregulated, while no alterations are shown in the phosphorylation state or nuclear translocation of β-catenin, thus suggesting that the non-canonical Wnt pathway functions in the process of HSC activation (44). Furthermore, previous research has shown that exogenous Wnt5a can promote the proliferation of HSCs, and the Wnt signaling pathway promotes liver fibrosis by activating HSCs (45). In addition, after treatment with a Wnt signaling inhibitor, HSCs can be restored to a resting state and this inhibitor can also induce apoptosis in cultured HSCs, indicating that blocking Wnt signaling may be a potential treatment for liver fibrosis (46). Thus, it was speculated that piRNAs selected in the present study have potential research value in metabolism, IR, liver fibrosis and other pathways of NAFLD.

To the best of our knowledge, the present study is the first to demonstrate that there were numerous differential expressed piRNAs in NAFLD using piRNA chip detection. Thus, the present study established a foundation for further research on NAFLD. However, there are limitations to the present study. Firstly, this study was focused on the piRNA genes that were positively expressed in NAFLD, and did not examine the importance of negatively expressed piRNA genes. Secondly, as a key protein involved in the biological role of piRNA, the expression of piwi protein in NAFLD was not assessed in the present study. Therefore, if negative piRNA genes and Piwi protein expression levels were examined, the present results would be more comparative. Thus, future research is warranted to further identify the potential mechanisms of NAFLD.

Acknowledgements

Arraystar MM9 piRNA array and microarray experiments were performed by Kang Chen Bio-Tech, Inc.

Funding

No funding was received.

Availability of data and materials

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

Authors' contributions

XM conceived the study protocol, processed the samples and wrote manuscript. YH, YD and LS detected the biochemical indexes and performed RT-qPCR. XZ and MK analyzed the data and performed the statistical analysis. CL conceived the study protocol, critically revised the manuscript and approved its final version. All authors read and approved the final manuscript.

Ethics approval and consent to participate

The ethical review of animal experiments was approved by the Ethics Committee of Affiliated Hospital of The Southwest Medical University (Luzhou; approval no. 201906-6).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References