Evaluation of the severity of nonalcoholic fatty liver disease by the analysis of serum exosomal miRNA expression

Noninvasive modalities exhibit limited diagnostic performance for evaluating the severity of nonalcoholic fatty liver disease (NAFLD). MicroRNAs (miRNAs) serve as useful biomarkers for diagnosing and monitoring disease progression and treatment response. Here, we evaluated whether serum exosomal miRNAs can be used to clinically predict the severity of NAFLD. Exosomal miRNAs were isolated from the sera of 41 NAFLD patients (diagnosed using biopsy) for array proling. To compare differentially expressed miRNAs, the degree of severity was determined using inammation, steatosis, ballooning, and NAFLD activity score (NAS). The correlation of miRNAs with clinical and biochemical parameters was analyzed. The correlation between miRNA and mRNA expression was analyzed by comparing our miRNA expression data with publicly available mRNA expression data. Twenty-ve, eleven, thirteen, and fourteen miRNAs correlated with inammation score, steatosis score, ballooning score, and NAS, respectively. Thirty-three signicant correlations were observed between twenty-seven miRNAs and six clinical variables (age, ALP, AST, ALT, GGT, and NAS). In brosis, 52 and 30 interactions corresponding to high miRNA-low mRNA and low miRNA-high mRNA expression, respectively, were observed. We demonstrated that serum exosomal miRNAs can be used to evaluate NAFLD severity and also delineate potential targets for NAFLD treatment.


Introduction
Nonalcoholic fatty liver disease (NAFLD) comprises a broad spectrum of diseases, including nonalcoholic fatty liver, nonalcoholic steatohepatitis (NASH), cirrhosis, and hepatocellular carcinoma (HCC) 1 . NASH is a progressive disease resulting in in ammation and hepatocyte injury, whereas nonalcoholic fatty liver or simple steatosis are benign conditions. To evaluate the severity of NAFLD, NAFLD activity score (NAS) was developed by the NASH Clinical Research Network, and is a composite score of steatosis, in ammation, and hepatocyte ballooning 2 . Fibrosis is another important histologic feature, which is the most relevant factor associated with prognosis 3 . To analyze the histological ndings, liver biopsy is the gold standard for the diagnosis of NASH and evaluation of NAFLD disease severity 4 . However, liver biopsy has certain disadvantages, such as sampling error, inter-/intra-observer variation; there is a risk for complications such as pain, infection, and bleeding 5 . As its global prevalence has increased remarkably in recent years, it is impossible to perform this invasive procedure in all patients 6 . Therefore, noninvasive evaluation of disease severity of NAFLD has been developed using serologic and imaging biomarkers [7][8][9] . However, the clinical application of such noninvasive biomarkers has several disadvantages, such as high cost, di culty in availability, low accuracy and reliability, and limited validation. In addition to the conventional radiological, serological, and pathological markers, exosomes and microRNAs (miRNAs) are promising biomarkers for the diagnosis and prognosis of NAFLD 10,11 . Exosomes contain several types of cargo molecules for cell signaling, such as proteins, lipids, and nucleic acids including mRNA, miRNAs, and other non-coding RNAs 12 . Among the non-coding RNAs, miRNAs are small-sized ones that regulate epigenetic gene expression by binding and suppressing their target mRNAs 13 . Therefore, miRNAs are involved in the progression of several diseases and have been studied as potential therapeutic targets 14 . Many miRNAs have been demonstrated to be involved in the progression of NAFLD and are emerging as novel biomarkers for distinguishing the different stages of NAFLD disease severity 15 . A recent meta-analysis identi ed several miRNAs, including miR-34a, miR-122, and miR-192, as potential biomarkers of NAFLD and NASH 16 .
In this study, we conducted a microarray analysis of serum exosomal miRNAs in patients with NAFLD diagnosed by a biopsy. Following bioinformatic analysis, we observed important implications of exosomal miRNAs, which re ected disease severity, such as the grade of NAS score or degree of liver brosis. We also demonstrated the correlation between miRNAs and clinical and biochemical parameters. With the help of publicly available datasets, we observed miRNA-mRNA interactions that may help elucidate future targets for improvements in diagnostic and therapeutic performance.

Baseline characteristics
Demographic and laboratory ndings and histopathological observations of 41 patients who were diagnosed with NAFLD based on the histological examination are presented in Table 1. The prevalence of diabetes/impaired fasting glucose, hypertension, and dyslipidemia was 68.3%, 46.3%, and 31.7%, respectively. The median age was 55 years, and the patients were predominantly female (68.3%). The median body mass index (BMI) was 29.94 kg/m 2 . We observed that 22 patients had NASH [demonstrated hepatocyte ballooning (53.7%)] and 16 patients had advanced brosis. The NAS score was distributed between 2 and 8 points, and 17 patients had a score of 5 or higher. The brosis stages of 0, 1, 2, 3, and 4 were demonstrated by 8, 10, 7, 12, and 4 patients, respectively. Sixteen patients (39.0%) demonstrated an advanced stage of brosis. The patients were divided into groups with low, medium, and high NAS scores. Analysis of expression patterns based on the degree of in ammation, steatosis, ballooning, and total

NAS
We analyzed the three groups and observed that 25, 11, 13, and 14 miRNAs correlated with the in ammation score, steatosis score, ballooning score, and total NAS, respectively. The expression of the miRNAs, which were observed in this study, was expressed as a heatmap, and a column annotation bar describing the scores of each patient was added. (Fig. 1)

Patterns of correlation of miRNA with clinical parameters in patients with NAFLD
Since correlations between the expression of certain miRNAs and disease severity were observed, we decided to evaluate the correlation of the miRNAs with clinical parameters. We selected certain clinical demographic parameters and laboratory ndings to evaluate their correlation with the miRNAs. Sixteen continuous variables were selected (age, BMI, hemoglobin, WBC, platelet, albumin, AST, ALT, ALP, GGT, bilirubin, PT, BUN, Cr, CRP, and NAS score). miRNAs with a statistical power of 0.9 or more were selected. All correlations between the 147 expression patterns of 117 miRNAs and 15 variables are listed in Supplementary Table 1 online. We found that negatively correlated with miRNA expression, whereas ALP, ALT, AST, GGT, and NAS scores positively correlated with miRNA expression. When we ltered the r 2 based on a cut-off of 0.35 or higher, 33 correlation patterns between the expression of miRNA with six variables were detected. (Table 2  We constructed a Circos plot to show a comprehensive correlation pattern of miRNA expression in each chromosome. Six variables including brosis, NAS, presence or absence of diabetes, pathologic parameters (in ammation, steatosis, and ballooning score), and biochemical parameters (ALT and platelet counts) were selected. Each variable was grouped into the categories of low or high scores or values ( brosis stage, 0-2 versus 3, 4; NAS score, < 5 versus ≥ 5; diabetics, yes versus no; in ammation score, score 1 versus 2 versus 3; steatosis score, score 1 versus 2 versus 3; ballooning score, score 0 versus 1 versus 2), and the correlation between miRNA expression and two biochemical parameters was analyzed (ALT and PLT). (Fig. 3) Two hundred miRNA expression patterns with low p-values were selected based on the t-test performed by comparison between two groups ( brosis, NAS, and diabetes) of each variable or ANOVA performed by comparison between two or three groups of each variable. (Supplementary Fig. 1) In each group, the miRNA expression was altered between two or three groups. Among the clinical parameters, ALT and platelet levels were correlated with each group as shown by two internal peaks. When the miRNA expression of ALT and platelet according to the group was positively correlated, a yellow or red peak appeared, and a negative correlation was indicated by green or sky blue. The correlation of miRNA expression at each genomic location was con rmed by the height of the peak.

miRNA-mRNA interaction networks
Co-expression networks between the correlated mRNAs and miRNAs were analyzed based on the miRNA expression data obtained in this study and with publicly available mRNA expression dataset, GSE89632. This dataset was selected as it contains brosis and NAS scores. The miRNAs demonstrating signi cant alterations in the expression patterns associated with total NAS score and brosis were selected. Using the same approach as a previous study in analyzing the GSE89632 dataset, we identi ed 385 and 359 genes with high and low expression, respectively. After considering that the expression of miRNA in the serum and liver tissue itself demonstrates a variable correlation, we matched two-by-two correlations between serum miRNA and mRNA expression patterns (higher miRNA-lower mRNA, lower miRNA-higher mRNA, lower miRNA-lower mRNA, and higher miRNA-higher mRNA). With respect to brosis, we observed 52 interactions corresponding to high miRNA-low mRNA expression and 30 interactions corresponding to low miRNA-high mRNA expression using the mirDIP program. One of the downregulated miRNAs, has-let-7b-5p, was related to a total of 30 mRNAs. (Fig. 4a) Five miRNAs were highly expressed in the high brosis samples, and 39 genes were downregulated. (Fig. 4b) With respect to the NAS score, 13 and 2 interactions were observed, corresponding to low miRNA-high mRNA (Fig. 4c) and high miRNA-low mRNA expression levels, respectively (Fig. 4d).

Discussion
As NAFLD is a broad-spectrum disease with pathologies ranging from simple steatosis to cirrhosis or HCC 17 , it is di cult to distinguish between the stages of NAFLD without liver biopsy, which is an invasive procedure. miRNAs are novel candidates for potential diagnostic and prognostic biomarkers of NAFLD because they are present in abundance in the liver and are associated with various kinds of liver diseases. 18 . We observed that several miRNAs showed signi cant differences inexpression in accordance with the severity of NAFLD and histological ndings.
In NAFLD patients diagnosed with NASH or advanced brosis, the brosis stage is determined to be 3 or 4, and the prognosis is typically poor compared to that of simple steatosis. 3 . In patients with NAFLD, NASH represents a severe form of hepatic damage due to the recruitment of pro-in ammatory immune cells, and it can eventually progress to cirrhosis and hepatocellular carcinoma. 19 Steatosis evaluation is important for the diagnosis of NAFLD, and change in steatosis is considered one of the indicators for improvement in clinical trials 20 . For the diagnosis of NASH with advanced brosis or evaluation of steatosis in patients, liver biopsy is the gold standard. However, it is impossible to perform a liver biopsy in all patients with NAFLD that account for 25-30% of the general population 2121212121 . Therefore, biomarkers for the noninvasive evaluation of NAFLD have been studied, such as serological and imaging biomarkers. 8 Some biomarkers are based on a single parameter and the others are developed in the form of a panel that combines several parameters. 7 Our study demonstrated that serum exosomal miRNAs showed signi cant differences according to the severity of NAFLD. NASH, an advanced form of NAFLD, is histologically characterized by in ammation, steatosis, and hepatocyte ballooning, and each histological feature is associated with different miRNAs. The diversity of miRNA expression pro les indicates that miRNAs are independently involved in the steatosis, in ammation, and ballooning process.
In addition, we identi ed correlations between laboratory data and expression of serum exosomal miRNAs. Alteration in the levels of liver enzymes, such as AST, ALT, ALP, and GGT, is important not only for the differential diagnosis of liver disease, but also for evaluating disease severity. 22 However, no study has reported on the relationship between liver enzymes and miRNAs. In this study, we identi ed four miRNAs that are positively correlated with liver enzymes: miR-133b, miR-4436a, miR-4709-3p, and miR-8079.
Although the NAFLD Activity Score (NAS) is de ned as the sum of steatosis, in ammation, and hepatocyte ballooning scores, when ANOVA analysis was performed, different miRNA expression patterns were found in NAS, steatosis, in ammation, and ballooning. However, in the Circos plot, the average miRNA expression level for each degree showed a similar pattern in each genomic region. Many miRNAs were observed corresponding to chromosomes 16, 17, and 19, which showed statistically signi cant differences based on the score, and fewer miRNAs were discovered corresponding to chromosomes 4 and 10. No miRNA correlated to the steatosis score was observed on chromosome 10. Interestingly, higher or lower correlation patterns between miRNA expression and two clinical parameters (ALT and PLT) were also detected with chromosomes 16, 17, and 19.
Functional network analysis between mRNA and miRNA is useful for evaluating the course of various diseases. 23,24 We identi ed relationships between the miRNA expression patterns observed in this study and the gene expression patterns demonstrated in the GSE89632 dataset. In the GSE89632 dataset, subjects were divided into healthy controls and patients with steatosis and NASH, and demonstrated characteristics of brosis (stage), steatosis (%), lobular in ammation (severity), ballooning score, and NAS 25 . The study methodology was similar to our strategy. Hence, we used this dataset for the comprehensive network analysis.
When it comes to hepatic brosis, we found that low expression of hsa-let-7b-5p correlated with high expression of genes associated with liver brosis. The interaction between pro-brosis genes and hsa-let-7b-5p has been identi ed in interstitial pulmonary brosis and renal brosis. 26,27 As for NAS score, low expression of miR-133b correlated with high expression of several genes. As miR-133b showed protective effects against allergic in ammation, 28 lower expression might increase hepatic in ammation. Taken together, the interaction between miRNA and mRNA is believed to play an important role in the progression of NASH.
The scarcity of publicly available gene expression data contributed to the limitation of our study. To the best of our knowledge, studies with similar approaches providing clinical parameters have not been previously performed, except for the GSE89632 study. The evidence demonstrated in our miRNA expression dataset needs to be validated in independent cohorts with larger sample sizes. The evidence and methods provided by us could be expanded to a comprehensive study of miRNAs and genes related to steatosis, brosis, ballooning, diabetes, and NAS. The small number of patients in each group constituting each variable is also a limitation. When ANOVA was performed with four variables, in ammation, steatosis, ballooning score, and NAS, it was di cult to ascertain a clear classi cation pattern by miRNAs describing each variable. Although the p-value indicating statistical signi cance was different for each variable, a minimum of 11 and a maximum of 25 miRNAs were identi ed that explained the classi cation of each variable. The data needs to be veri ed with larger patient cohorts in future studies.
We demonstrated differentially expressed miRNAs based on the groups of four variables with high, middle, and low values corresponding to the degree of severity (in ammation, steatosis, ballooning score, and NAS). Although NAS comprises three variables, different miRNAs have been identi ed, which are associated with each variable. Therefore, the pathways to explain each variable might be different.
Comparative visualization was attempted with the construction of a Circos plot. MiRNAs that were signi cantly different with respect to each variable showed similar expression patterns in similar genomic regions. The integration of miRNA and mRNA expression analyses and network analysis also enabled us to interpret the differentially regulated genes from the perspective of systems biology in NASH. Using our methodology, co-expression networks of miRNAs and mRNAs could help reveal the pathological pathways of NAFLD, as well as provide new insight into several biological pathways related to liver functions. These pathways may also be used in the diagnosis of many liver diseases.

Study population
A total of 41 patients who were diagnosed with NAFLD following biopsy, were enrolled. Patients with other chronic liver diseases or excessive alcohol consumption were excluded. Past medical history and clinical parameters were recorded, and biochemical analyses was performed during the admission period for liver biopsy. Biochemical analyses included the evaluation of hemoglobin, white blood cells (WBC), platelets, aspartate aminotransferase (AST), alanine aminotransferase (ALT), alkaline phosphatase (ALP), gamma-glutamyl transferase (GGT), total bilirubin, albumin, prothrombin time (PT), blood urea nitrogen (BUN), creatinine (Cr), and C-reactive protein (CRP). The approval for this study was obtained from the institutional review board of Korea University Guro Hospital (2016GR0302). All participants provided written informed consent and agreed to provide their sera and clinical data for this study. This study was conducted according to the Declaration of Helsinki, and all participating authors reviewed and approved the nal version of the manuscript.

Histopathologic evaluation
Two specimens of liver tissues were obtained via percutaneous liver biopsy. Following xation, the sections were embedded in para n blocks and were stained with hematoxylin and eosin and Masson's trichrome stain. Two experienced pathologists (HAK and BHK) evaluated each sample according to the NASH Clinical Research Network (NASH CRN) histological scoring system. 2 NAS was calculated as the sum of steatosis (0-3), in ammation (0-3), and hepatocytes ballooning scores (0-2). Fibrosis was classi ed into stages 0-4.
Serum collection, exosomal isolation, RNA extraction and microarray analysis Following exosomal isolation from the serum samples of patients with NAFLD, RNA extraction and microarray analysis were performed in accordance with a previously published protocol 29 . Microarray analysis was performed using an Affymetrix GeneChip™ miRNA 4.0 array, and image and signal data were extracted using GeneChip™ Scanner 3000DX, and Transcriptome Analysis Console 4.0. Normalization of each type of data was performed using Robust Multichip Average (RMA) and Detection Above Background (DABG) algorithms in Affymetrix Expression Console software. Normalized data were represented as a matrix, following which they were converted to the R dataframe; and statistical analysis was performed.
Identi cation of differentially expressed miRNAs in each variable and in a total of NAS scores We used analysis of variance (ANOVA) and ANOVA functions in R to identify differentially expressed miRNAs among each of the three groups based on the in ammation score (0/1, 2, and 3), steatosis score (0/1, 2, and 3), ballooning score (0, 1, and 2), and total NAS (1-3, 4/5, and 6-8).
The heatmap and hierarchical clustering plot were obtained using the 'pheatmap' package in R. The average expression of miRNAs in each group was calculated by 'rowMeans' function in R. The correlation between miRNA expression and clinical parameters of patients was evaluated using 'pingouin' package and 'corr' function in Python. The Circos plots were constructed using the 'circlize' package in R.

Assessing miRNA and mRNA interaction
To study the interaction between the miRNA expression evaluated in this study with the expression of the corresponding target mRNA or gene, we analyzed gene expression using the publicly available GSE89632 dataset. In a previous study, the GSE89632 dataset was selected as it depicts brosis and NAS score 25 . In this study, the results of 41 patients were compared to those of 24 healthy controls, 20 patients with simple steatosis, and 19 patients with NASH. Gene expression pro les were con rmed using the Illumina HumanHT-12 WG-DASL V4.0 R2 expression beadchip.

miRNA-mRNA interaction networks
For the prediction of targets modulate by miRNAs, we used the miRNA Data Integration Portal bioinformatics tool (http://ophid.utoronto.ca/mirDIP/). Cytoscape software was used to incorporate the identi ed miRNA and target genes to the interaction network (version 3.8.0).

Declarations
Grant Support: This study was supported by a National Research Foundation of Korea grant from the Korean government (the Ministry of Education, Science and Technology 2019M3E5D1A01068997 and 2018R1A2B2006183).