Distinct MicroRNAs Expression Profile in Primary Biliary Cirrhosis and Evaluation of miR 505-3p and miR197-3p as Novel Biomarkers

Background and Aims MicroRNAs are small endogenous RNA molecules with specific expression patterns that can serve as biomarkers for numerous diseases. However, little is known about the expression profile of serum miRNAs in PBC. Methods First, we employed Illumina deep sequencing for the initial screening to indicate the read numbers of miRNA expression in 10 PBC, 5 CH-C, 5 CH-B patients and 5 healthy controls. Comparing the differentially expressed miRNAs in the 4 groups, analysis of variance was performed on the number of sequence reads to evaluate the statistical significance. Hierarchical clustering was performed using an R platform and we have found candidates for specific miRNAs in the PBC patients. Second, a quantitative reverse transcription PCR validation study was conducted in 10 samples in each group. The expression levels of the selected miRNAs were presented as fold-changes (2−ΔΔCt). Finally, computer analysis was conducted to predict target genes and biological functions with MiRror 2.0 and DAVID v6.7. Results We obtained about 12 million 32-mer short RNA reads on average per sample and the mapping rates to miRBase were 16.60% and 81.66% to hg19. In the statistical significance testing, the expression levels of 81 miRNAs were found to be differentially expressed in the 4 groups. The heat map and hierarchical clustering demonstrated that the miRNA profiles from PBC clustered with those of CH-B, CH-C and healthy controls. Additionally, the circulating levels of hsa-miR-505-3p, 197-3p, and 500a-3p were significantly decreased in PBC compared with healthy controls and the expression levels of hsa-miR-505-3p, 139-5p and 197-3p were significantly reduced compared with the viral hepatitis group. Conclusions Our results indicate that sera from patients with PBC have a unique miRNA expression profile and that the down-regulated expression of hsa-miR-505-3p and miR-197-3p can serve as clinical biomarkers of PBC.


Introduction
MicroRNAs (miRNAs) are small endogenous RNA molecules of 19 to 24 nucleotides that control the translation and transcription of targeting mRNAs by base-pairing to the complementary sites [1] [2] [3] [4]. The expression of miRNAs in serum is reported to be stable, reproducible and consistent among individuals of the same species [5]. So far, specific expression patterns of serum miRNAs were identified as a fingerprint for numerous diseases and cancers [6] [5]. The serum miR-122 levels are elevated in patients with liver damage due to chronic hepatitis B (CH-B) and C infection (CH-C) [7] [8]. In addition, the miR-122 and miR-34a levels are positively correlated with the disease severity in CH-C and non-alcoholic fatty-liver disease [9]. However, there are some reports that miR-122 expression in healthy controls was significantly higher than that in patients with hepatitis C virus (HCV) infection [10]. Li et. al. described that 13 miRNAs were differentially expressed in hepatitis B virus (HBV) serum and that miR-25, miR-375 and let-7f could be used as biomarkers to separate a HBV-positive hepatocellular carcinoma (HCC) group from HBV-negative HCC [11]. However, little is known about the expression profile of miRNAs in autoimmune disease such as primary biliary cirrhosis (PBC).
PBC is female predominant, progressive autoimmune disease characterized by immune-mediated destruction of the intrahepatic bile ducts. The serological marker of PBC is the presence of antimitochondrial antibody (AMA) directed against the E2 subunit of the pyruvate dehydrogenase enzyme complexes located in the inner mitochondrial membrane [12] [13] [14]. The etiology of PBC is considered to be a combination of genetic predisposition and environmental triggers [15]. Particularly, concerning genetic predisposition, previous studies reported that common genetic variants at the HLA class II, IL12A, IL12RB2, STAT4, IRF5-TNPO3 and IKZF3 had significant associations [16] [17] [18] [19] [20] [21]. Recently, genome-wide association study in Japanese population showed TNFSF15 and POU2AF1 as susceptibility loci [22]. Several GWAS data suggested the important contributions of several immune pathways to the development of PBC. However, the results have differed among the study groups [21]. The diagnosis of PBC is established based on the following criteria: (1) biochemical evidence of cholestasis; (2) the presence of AMA; and (3) histopathologic evidence of nonsuppurative cholangitis and destruction of the interlobular bile ducts [23]. Though diagnostic criteria have been determined, the eventual progression to biochemically and clinically apparent disease is unpredictable. Many patients are recognized at an earlier stage of disease and respond well to medical therapy, while some patients will require liver transplantation [24] [25]. To revolutionize the diagnosis, treatment and prognosis of PBC, new biomarkers seem to be feasible, and miRNAs are emerging as highly tissue-specific biomarkers with potential clinical applicability [26] [6].
In this study, we employed a strategy of using Illumina small-RNA sequencing for the initial screening followed by quantitative reverse transcription PCR (qRT-PCR) validation to analyze serum samples, which were arranged in multiple trial and testing sets. Additionally, computer analysis was conducted to predict target genes and biological functions from the differentially expressed miRNAs in PBC. The results demonstrate that the unique expression pattern of serum miRNAs can serve as a noninvasive biomarker for the diagnosis of PBC.

Global analysis of miRNAs by deep sequencing
Circulating miRNAs were detected from human serum in 10 patients diagnosed with PBC and 15 non-PBC subjects with CH-B, CH-C and healthy controls, using Illumina GA IIx sequencing. Detailed clinical information is shown in table 1. We analyzed three samples by single-end deep sequencing on one lane and obtained about 12 million 32-mer short RNA reads on average per sample. After trimming the reads to exclude adaptor and tag sequences by using cutadapt, 9,245,752 high-quality reads per sample were subjected to analysis [27]. The mapping rates to miRBase were 16.60% on average and those to hg19 were 81.66% (Table 2). miRNA expression profile in serum affected of PBC We normalized the differential expression of miRNA count data using the trimmed mean of M values (TMM) normalization process and the number of individual miRNA reads was standardized by the total numbers of 1,000,000 reads in each sample [28]. Comparing the 4 groups (PBC, CH-C, CH-B and healthy control), the differential expression levels of miRNA were extracted using analysis of variance (ANOVA). 1594 miRNAs were detected by deep sequencing (Table S1). Due to the small number of miRNA detections, the statistical significance of 821 miRNAs could not be determined. The ANOVA was used to determine differentially expressed miRNAs and multiple comparisons procedure was applied to compare more than one pair of means at the same time. Therefore, in statistical significance testing in the remaining 773 miRNAs, the p-value by multiple comparisons was performed by calculating False discovery rate (FDR),0.1. The expression levels of 81 miRNAs were found to be differentially expressed in the 4 groups. Although, several types of clinical data were shown to be different from PBC, M2 negative for PBC-8, ANA positive for PBC-10 or past HBV infection for PBC-6, the histologies of all samples were characterized by PBC of chronic, nonsuppurative cholangitis affecting the interlobular and septal ducts. The heat map and hierarchical clustering demonstrated that the miRNA profiles from PBC clustered with those of CH-B, CH-C and healthy controls ( Figure 1). Of note, CH-B and healthy controls were not clearly distinguished.

miRNA validation study
We used qRT-PCR to verify the data obtained from the Illumina sequencing. The relative expression levels from 9 differentially expressed miRNAs were analyzed with the TaqMan MicroRNA assay (Applied Biosystems) or miRCURY LNA microRNA PCR system (Exiqon). In addition to the 25 samples used for Illumina deep sequencing, 10 cases of miRNA expression in PBC, CH-C, CH-B and healthy controls were quantified ( Table 4). The expression levels of selected miRNAs detected by qRT-PCR were normalized to spiked-in cel-miR-39 and presented as fold-changes (2 2DDCt ) above those of the CH-C-5. Among 9 miRNAs, the quantity of four miRNA expressions could be determined. The circulating levels of hsa-miR-505-3p, miR-197-3p and miR-500a-3p (p,0.01) were significantly decreased in PBC compared with healthy controls and the expression levels of hsa-miR-505-3p, miR-139-5p and miR-197-3p were significantly reduced compared with the viral hepatitis group (Figure 2). The quantities of miRNA expression by qRT-PCR were supported by the data of Illumina sequencing. Of note, for samples with only trace amounts, quantification of the reaction products after 30 cycles is uncertain. The quantity of 5 remaining miRNA expressions indicated the Ct values under 30 cycles even in the higher amount group or were undetected.

miRNA target genes and Gene ontology (GO) analysis
To determine possible target genes for the differentially expressed miRNAs, we searched the miRror 2.0 database for predicted miRNA targets in human. The ranking of miRror was according to the miRror Internal Score (miRIS). This score is a balance between (i) the proportions of predicting miRNA-target prediction databases (MDBs) out of all tested MDBs and (ii) the fraction of the potentially regulated genes from the entire input genes [29]. A total of 75 genes were predicted as targets for 9 differentially expressed miRNAs by Illumina sequencing of hsa-miR-1273g-5p, miR-33a-5p, miR-3960, miR-766-5p, miR-505, miR-30b, miR-139-5p, miR-197 and miR-500a in PBC serum (Table 5). Then, these 75 genes were submitted to Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7, which was used for the GO biological process categorization, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and BIOCARTA pathway. The predicted target genes involved biological processes such as blood circulation (5.4%), circulatory system process (5.4%), positive regulation of I-kappaB kinase/NF-kappaB cascade (4.1%), regulation of phosphoprotein phosphatase activity (2.7%), cell volume homeostasis (2.7%), regulation of phosphatase activity (2.7%) and cellular amino acid derivative catabolic process (2.7%), some cellular components such as cytosol (13.5%), cell fraction (12.2%), membrane fraction (9.5%) and protein serine/threonine phosphatase complex (4.1%), and some molecular functions such as ion binding (33.8%), Cation binding (33.8%), metal ion binding (33.8%) and transition metal ion binding (25.7%) ( Table S2). The functional annotation analysis of BIOCARTA showed that the genes of catenin (cadherinassociated protein), alpha 1 and similar to breast cancer antiestrogen resistance 1 predicted target genes of the listed miRNAs, and played a role in cell-to-cell adhesion signaling pathway ( Figure  S1). The KEGG pathway indicated that the genes of baculoviral IAP repeat-containing 2, protein phosphatase 3 catalytic subunit beta isoform and tumor necrosis factor ligand superfamily member 10 were related to apoptosis ( Figure S2).

Discussion
MiRNA changes in the liver have been reported in diseases such as HCC or chronic viral hepatitis. However, there is only limited information about their detection in blood and their correlations in PBC patients. The current study provides the first evidence that PBC is associated with altered miRNA expression. We have demonstrated that a number of miRNAs, especially hsa-miR-505-3p and miR-197-3p, were significantly differentially expressed in patients with PBC, leading to a unique miRNA expression profile in the diseased liver. Recently, many studies have examined several PBC associations with genes and there have been significant differences in the genetic risk loci reported [21] [30]. Therefore more carefully constructed studies will be needed to clarify, the pathogenesis of PBC, and the study of these differentially expressed miRNAs could serve in identifying biomarkers or lead to a better understanding of the underlying molecular mechanism that perpetuates PBC.
In our study, miRNA Illumina deep sequencing was first used to screen 10 PBC patients' sera. We were able to match the sample's sex because of particular importance for X-linked miRNA [31]. Then, qRT-PCR was used to confirm the result of deep sequencing.
Quantitative differential expression analysis identified a 81-miRNA signature distinguishing PBC, CH-C, CH-B and healthy controls. A hierarchical clustering analysis was performed utilizing the 81-miRNAs and their patterns separated the PBC from viral hepatitis and healthy controls. In addition, there were three subgroups (PBC-1,2,4,5, PBC-3,6,10 and PBC-7, 8,9) in the PBC cluster. When comparing the subgroups, the PBC-3, 6,10 group showed an expression pattern that differed from those of the other two subgroups. As for the clinical background, PBC-6 patient that had been infected with hepatitis B virus were HBsAg negative and anti-HBc and anti-HBs positive, and PBC-10 patient had positive antinuclear antibodies (ANA) titers of 1:2560 or greater, while other patients had no serological evidence of HBV infection and no positive ANA titers of 1:160 or greater. However, there was no clear difference between the three subgroups in terms of the clinical stage. At present, there is no conclusive proof whether  Figure 1. Heat map and hierarchical clustering. Individual miRNA expression were calculated by R platform and heat map was computed and described using a function of heatmap.2 in gplots. It uses hierarchical clustering with Euclidean distance; Pearson Linear Correlation and Ward's method to generate the hierarchical tree [56]. ANOVA was applied to extract differentially expressed miRNAs and adjustment of the p-value by clinical or pathological differences can be found in clustering subgroups. Nine miRNAs were confirmed to be significantly differentially expressed between the PBC group and viral hepatitis group or healthy control by Illumina deep sequencing. Among these 9 miRNAs, the serum levels of hsa-miR-505-3p and miR-197-3p were significantly lower in patients with PBC than in those with viral hepatitis and healthy controls, hsa-miR-139-5p was lower in patients with PBC than in those with viral hepatitis and miR-500a-3p were lower in patients with PBC than in healthy controls. Of note, we conducted qRT-PCR on the sera of some PBC who had already been treated with ursodeoxycholic acid. The serum levels of miRNA showed improvements in some samples (data not shown). Accordingly, quantifying these miRNAs may yield reliable diagnostic information. However, one problem is that the quantification of miRNA in this study used standardization by the total numbers of 1,000,000 reads in the deep sequencing or a comparative method in qRT-PCR. In other words, it was assumed that the same amount of miRNA was contained in each serum sample. Therefore, if one miRNA is quantified in a single specimen, we will not be able to accurately assess the result.
Specific circulating miRNA profiles have been reported for various diseases [32] [5] [33] [34] [35]. These circulating miRNA profiles have been described as correlating with differentially expressed miRNA in diseased tissue, such as liver injured by drugs or stomach afflicted with gastric cancer [36] [37]. Moreover, some disease-specific profiles can inform both the diagnosis and prognosis [38] [39]. Therefore, to determine if any of these diferentially expressed miRNAs could lead to better understanding of the molecular mechanism that perpetuates PBC, we examined for gene targets that may be reflected by this particular miRNA expression signature. Of the several target prediction algorithms prepared, we selected mirror 2.0. There has been evidence that a seed region of miRNA positioned within a limited range in the 39 UTR of a target gene degrades the mRNA function [40]. We predicted 75 genes as targets for 9 differentially expressed miRNAs and conducted a functional analysis of DAVID. This analysis revealed that the genes of catenin (cadherin-associated protein), alpha 1 and similar to breast cancer anti-estrogen resistance 1 predicted target genes of the listed miRNA and played a role in cell-to-cell adhesion signaling, and the genes of baculoviral IAP repeat-containing 2, protein phosphatase 3 catalytic subunit beta isoform and tumor necrosis factor ligand superfamily member 10 were related to apoptosis. The onset of autoimmune disorders with PBC can be linked to apoptosis. A previous report described that the expression of TRAIL receptors is up-regulated by an increased bile acid level and that the serum level of soluble TRAIL is elevated, which may be involved in the development and progression of PBC [41,42]. However, further work will be required so that these miRNAs can serve not only as biomarkers but also for the elucidation of the pathogenesis of PBC.
GO analysis provides representations of biological annotations using precisely defined terms [43]. A previous report has described a number of genes involved in the signaling, regulation of I-kappaB kinase/NF-kappaB cascade and homeostasis that are associated with PBC [44] [45,46]. Additionally, our study indicated the biological processes, cellular component and molecular functions affected by the target genes included those associated with cell or membrane fraction, various kinds of ion binding and protein serine/threonine phosphatase complex, all of which are potentially related to PBC. Further studies will need to examine the relationship between differentially expressed miR-NAs, both in serum and liver tissue, and target genes, which may provide more insights into the role of miRNAs in the pathology of PBC.
In conclusion, our results indicate that sera from patients with PBC have a unique miRNA expression profile compared to viral hepatitis and healthy controls and down-regulated expression of hsa-miR-505-3p and 197-3p may represent new clinical biomarkers in PBC. This study suggests that the amounts of miRNAs in serum have potential as diagnostic and prognostic biomarkers for PBC.

Patients and sample processing
We included sera of 10 patients with PBC who were treatmentnaïve, sera of 5 patients with CH-B, sera of 5 patients with CH-C and sera of 5 healthy controls in this study. Initially these serum samples were enrolled to be analyzed by the Illumina miRNA deep sequencing (Illumina). The diagnosis of all cases was based on internationally established criteria [23].

Library preparation and Illumina sequencing
A ten ml venous sample was collected from each participant. The whole blood was separated into serum and cellular fractions by centrifugation at 2,500 r.p.m. for 10 min, followed by 10 min centrifugation at 10,000 r.p.m. to completely remove cell debris. The supernatant serum was stored at 220uC until analysis. Total RNA was extracted from 800 ml of serum using Trizol LS (Invitrogen, Carlsbad, CA). The libraries were constructed from total RNA using the TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, CA) following the manufacturer's protocol. Briefly, RNA 39 and 59 adapters were ligated to target microRNAs in two separate steps. Reverse transcription reaction was conducted to the ligation products to create single stranded cDNA. The cDNA was amplified by PCR using a common primer and a primer containing the index sequence. One ml of each library was loaded on an Agilent Bioanalyzer (Agilent, Santa Clara, CA) to check the size, purity, and concentration. Libraries were sequenced on an Illumina GA IIx (SCS 2.8 software; Illumina, SanDiego, CA), with a 32-mer single end sequence. Image analysis and base calling were performed using RTA 1.8 software.

Sequence and statistical analysis
Raw miRNA sequence reads were conducted as a quality check and the 39 and 59 adapter sequences were removed by cutadapt while discarding reads shorter than 20 nucleotides [27]. The sequence reads were mapped with miRBase (Release 18) and UCSC (hg19) by use of bwa (0.5.9-r16), allowing one nucleotide base mismatch [47] [48].
Digital expression levels were normalized by taking into account the length of miRNAs and the total number of miRNA reads generated in each library using TMM normalization [28]. Read counts of each identified miRNA was normalized to the total number of miRNA reads, and then the ratio was multiplied by a constant set to 1610 6 in this study. ANOVA was applied to extract differentially expressed miRNAs among the four groups. Adjustment of the p-value by multiple comparisons was performed by calculating FDR [49]. Those miRNAs with FDR,0.1 were extracted as differentially expressed and used in the following analysis. Hierarchical clustering was performed using an R platform and a heat map described as using a function of heatmap.2 in gplots [50].

qRT-PCR validation study
In addition to 25 samples analyzed by Illumina sequencing, five more serum samples of CH-B, CH-C and healthy controls (a total of 10 samples in each group) were used in qRT-PCR validation study. We followed the protocol previously reported by Mitchell et al. to determine the endogenous miRNA levels with spiked-in miRNA. Spiked-in miRNA was designed against Caenorhabditis elegans microRNA-39 (cel-miR-39) (59-UCA CCG GGU GUA AAU CAG CUU -39) and was synthesized by Sigma Aldrich Japan [32]. After total RNA isolation from 300 ml serum, reverse transcription was conducted using a TaqMan miRNA RT kit for identification of the cel-miR-39 expression (Applied Biosystems) with 5 fmol/ml for the internal control. qRT-PCR were conducted for detection of hsa-miR-1273g-5p, miR-505-3p and miR-139-5p in 20 ml PCR reactions using TaqMan MicroRNA assay with StepOne Plus detection system at 50uC for 2 min and 95uC for 10 min, followed by 40 cycles of 95uC for 15 s and 60uC for 1 min (Applied Biosystems). For detection of hsa-miR-33a-5p, miR-3960, miR-766-5p, miR-30b-3p, miR-197-3p and miR-500a-3p expression, we used the Exiqon system. Total RNA was reverse transcribed using the miRCURY LNA TM Universal RT miRNA PCR, Polyadenylation and cDNA synthesis kit (Exiqon). cDNA diluted 506 was assayed in 10 ml PCR reactions according to the protocol for miRCURY LNA TM Universal RT miRNA PCR with StepOne Plus detection system at 95uC for 10 min, followed by 40 cycles of 95uC for 10 s and 60uC for 1 min (Exiqon). The data were analyzed by the 2 2DDCt method.

Statistical methods
Expression levels of the selected miRNAs detected by qRT-PCR were normalized to cel-miR-39 and presented as the foldchange (2 2DDCt ) above the control (CH-C-5): DDCt = (Ct miRNA -Ct cel-miR-39 ) patients -(Ct miRNA -Ct cel-miR-39 ) CH-C-5 . Results for normally distributed continuous variables are given as means (6standard errors of the mean) and compared between groups by Student's ttest. Results for non-normally distributed continuous variables are summarized as medians (interquartile ranges) and were compared by Mann-Whitney U test.
In silico analysis of miRNA target gene The number of candidate genes and the number of miRNAs are indicated for each of the major MDBs. We selected as the search mode: miR2Gene; and as the search parameters of organism: human; and of selected tissue: all. Advanced parameters were inputted, cutoff: 0.01; database hits: 2; and target counts: 3. We created a list of common target genes for miRNAs. Then, these common targets were annotated by an annotation tool at the DAVID v6.7 (January 2010 release, http://david.abcc.ncifcrf. Table 3. Differentially expressed miRNAs in serum from PBC patients compared with the second group (CH-C, CH-B, Healthy).    gov/) [52]. DAVID can detect functional enrichment of a gene list based on the GO terms, KEGG pathway and BIOCARTA pathway. Differences were considered significant when the P value was less than 0.05.

Ethics statement
This study was approved by the Ethics Committee of the Tohoku University School of Medicine (2010-404) and written informed consent was obtained from each individual. Figure S1 The pathway of cell-to-cell adhesion signaling. The functional annotation analysis of BIOCARTA showed that the genes of catenin (cadherion-associated protein), alpha 1 and similar to breast cancer anti-estrogen resistance 1 played roles in this pathway. The stars indicate the related genes. (TIFF) Figure S2 The pathway of apoptosis. The functional annotation analysis of BIOCARTA showed that the genes of baculoviral IAP repeat-containing 2, protein phosphatase 3 (formerly 2B), catalytic subunit, beta isoform and tumor necrosis factor (ligand) superfamily, member 10 was related to apoptosis. The gene is indicated with the stars. (TIFF)

Supporting Information
Table S1 The list of differential expression levels of miRNA in each sample. (XLS)