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.

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.

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.

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.

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

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.

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

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 ddH₂O 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.

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.

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

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

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.

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

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

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

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.

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.