Identification of Translational microRNA Biomarker Candidates for Ketoconazole-Induced Liver Injury Using Next-Generation Sequencing

Abstract

Drug-induced liver injury (DILI) is a leading cause of acute liver failure. Reliable and translational biomarkers are needed for early detection of DILI. microRNAs (miRNAs) have received wide attention as a novel class of potential DILI biomarkers. However, it is unclear how DILI drugs other than acetaminophen may influence miRNA expression or which miRNAs could serve as useful biomarkers in humans. We selected ketoconazole (KCZ), a classic hepatotoxin, to study miRNA biomarkers for DILI as a proof of concept for a workflow that integrated in vivo, in vitro, and bioinformatics analyses. We examined hepatic miRNA expression in KCZ-treated rats at multiple doses and durations using miRNA-sequencing and correlated our results with conventional DILI biomarkers such as liver histology. Significant dysregulation of rno-miR-34a-5p, rno-miR-331-3p, rno-miR-15b-3p, and rno-miR-676 was associated with cytoplasmic vacuolization, a phenotype in rat livers with KCZ-induced injury, which preceded the elevation of serum liver transaminases (ALT and AST). Between rats and humans, miR-34a-5p, miR-331-3p, and miR-15b-3p were evolutionarily conserved with identical sequences, whereas miR-676 showed 73% sequence similarity. Using quantitative PCR, we found that the levels of hsa-miR-34a-5p, hsa-miR-331-3p, and hsa-miR-15b-3p were significantly elevated in the culture media of HepaRG cells treated with 100 µM KCZ (a concentration that induced cytotoxicity). Additionally, we computationally characterized the miRNA candidates for their gene targeting, target functions, and miRNA/target evolutionary conservation. In conclusion, we identified miR-34a-5p, miR-331-3p, and miR-15b-3p as translational biomarker candidates for early detection of KCZ-induced liver injury with a workflow applicable to computational toxicology studies.

Drug-induced liver injury (DILI) is the major reason for clinical trial termination in drug development and for drug withdrawal from the market (Kullak-Ublick et al., 2017). DILI is also a leading cause of acute liver failure (Lee, 2013; Weaver et al., 2020). Prediction of DILI using animal models has not always been successful; around 50% of drugs failed in clinical trials due to hepatotoxicity that was not detected in animal studies according to an analysis of 150 toxicity-inducing drug candidates (Olson et al., 2000). Serum levels of liver enzymes, such as alanine transaminase (ALT) and aspartate transaminase (AST), are considered standard diagnostic markers for DILI; however, elevation of serum enzyme levels is not always associated with significant liver injury (Baillie and Rettie, 2011). Thus, more reliable and translational biomarkers are needed for accurate prediction and early detection of DILI in preclinical investigations and clinical applications.

With the advancement of genomics technologies, increasing efforts have been directed toward identifying biomarkers using next-generation sequencing (NGS) beyond microarrays. Differentially expressed mRNA molecules are among the most common genomic biomarkers used for predicting drug safety, however, they have met challenges for DILI risk estimation (Kiyosawa et al., 2009). microRNAs (miRNAs) have received wide attention as a novel class of potential DILI biomarkers (Hornby et al., 2014; Krauskopf et al., 2015). miRNAs are a class of small non-coding RNA molecules (18–22 nt) that mostly function in post-transcriptional gene silencing and are involved in virtually every aspect of biology and pathology (Koturbash et al., 2012); they have been found in a number of bodily fluids including plasma, saliva, urine, bile, and breast milk, which potentially allows for non-invasive detection in clinical applications (El-Mogy et al., 2018; Weber et al., 2010). Circulating miRNAs can be produced via active secretion or passive leakage from tissues; they are often carried in membrane vesicles (eg, exosomes and apoptotic bodies) and thus protected from RNase degradation (Nik Mohamed Kamal and Shahidan, 2019). Many miRNAs also display specificity in disease subtypes and stages (Ludwig et al., 2016). For example, serum miR-122 and miR-192 were shown to predict hepatotoxicity induced by acetaminophen in a preclinical mouse model (Wang et al., 2009). Dysregulation of miRNAs has been associated with many diseases including liver diseases, which may show clinical symptoms similar to DILI (Schueller et al., 2018; Szabo and Bala, 2013). Moreover, altered miRNA expression appears to signal DILI development at an earlier stage than conventional DILI biomarkers (Zhang et al., 2010). Importantly, miRNAs are genetically conserved among species, which facilitates the extrapolation of discoveries in animal models to humans (van Rooij et al., 2012). These characteristics greatly encourage the exploration of miRNAs as important biomarkers for risk assessment in drug discovery and development. However, it is still unclear which miRNAs may serve as mechanism-driven translational biomarkers for DILI drugs other than acetaminophen.

Ketoconazole (KCZ) has been used for decades as an anti-fungal medication, and it has been shown to cause fatal hepatotoxicity and adrenal insufficiency (Garcia Rodriguez et al., 1999). Also, KCZ is a potent CYP3A4 inhibitor and may thereby cause serious drug-drug interactions (U.S. Food and Drug Administration, 2013). For these reasons, in 2013, the U.S. FDA and European Medicines Agency banned oral KCZ use for all but life-threatening fungal infections while recommending close monitoring for development of liver toxicity (European Medicines Agency, 2013; U.S. Food and Drug Administration, 2013). The effects of KCZ treatment on miRNA expression in the liver are unknown. Therefore, we chose KCZ as a classic hepatotoxin in our initial efforts to identify miRNA biomarkers for DILI. In this study we employed miRNA-sequencing (miR-seq) to analyze global miRNA expression in response to KCZ exposure in rat liver samples across multiple doses and treatment durations. We identified candidate miRNA biomarkers that were dysregulated by KCZ and specifically associated with a KCZ histopathological phenotype. In addition, we functionally characterized the candidate miRNA biomarkers for their gene targeting, target functions, and miRNA/target conservation among species using a set of comprehensive in silico approaches.

MATERIALS AND METHODS

Tissue sample collection, animal treatment, and ethics statement

Liver tissue samples were obtained by the U.S. FDA—National Center for Toxicological Research from the Japanese Toxicogenomics Project under a research collaboration agreement. Animal treatments for the original project were described previously (Uehara et al., 2010). Animal study protocols were approved by the National Institute of Health Sciences’ Ethics Review Committee for Animal Experimentation in Japan. All experiments were conducted according to the approval guidelines. Briefly, male Sprague Dawley rats (6-week old) were orally administered with KCZ at the low (10 mg/kg), middle (30 mg/kg), or high dose (100 mg/kg) per day for 3, 7, and 14 days, respectively. Rats treated with 0.5% methylcellulose (a vehicle for KCZ) for each of the three durations served as a corresponding control for that timepoint. Each treatment condition indicated a unique combination of a dose (control, low, middle, or high) and a timepoint (3, 7, or 14 days), totaling 12 conditions in this study. Rats were sacrificed and livers were collected 24 h after the final dosing. Liver tissue samples were immediately immersed in RNAlater RNA-stabilization solution (ThermoFisher, Waltham, Massachusetts) upon collection and were then stored at −80˚C with continuous temperature monitoring. Five rats were included at each treatment condition; liver samples from three out of the five animals were processed for miR-seq. Samples treated for 28 days in the original project were unavailable for miR-seq and thereby not considered in this study. The high dose was designed in the original study to induce minimal levels of toxicity over the 28-day treatment period to permit discovery of early biomarkers.

Histopathology and blood chemistry data analyses

Histopathology and blood chemistry data were downloaded from the data archive https://dbarchive.biosciencedbc.jp/en/open-tggates/download.html (last accessed May 2020). Pathological phenotypes and serum enzyme elevation in rats with matching miR-seq data were analyzed and graphed. Severity of liver injury phenotypes were color-coded.

miR-seq analysis

miR-seq experiments and analyses were conducted as described previously (Dweep et al., 2017). Briefly, total RNA was extracted from approximately 15 mg liver tissues using RNeasy Plus Universal Mini kits (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Quality of total RNA extracted from the liver tissue samples was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, California), and all samples yielded an RNA integrity number (RIN) greater than 8.6. Illumina TruSeq Small RNA Library Preparation kits (Illumina, San Diego, California) were used to construct miR-seq libraries from 1 μg extracted total RNA per sample following users’ manuals. Libraries were barcoded and pooled together for miR-seq on an Illumina HiSeq-2000 platform, resulting in about 10 million 50-bp single end reads per sample. Sequencing service was provided by the Genomics and Microarray Core Facility at the University of Texas Southwestern Medical Center (Dallas, Texas). Sequencing data were then processed for 3′ adaptor trimming and quality check, which were followed by mapping, annotation, and quantification using the rat miRNA dataset from miRBase (Release 21) http://www.mirbase.org (last accessed July 2020) and the miRDeep2 pipeline (Friedländer et al., 2012). Based on quality check, miR-seq data for one control sample on Day 14 were questionable and thereby excluded from subsequent analyses.

miRNA differential expression analysis

Analysis to discover differentially expressed miRNAs (DEMs) was conducted in RStudio (version 1.2.5019) using the DESeq2 (v1.24.0) package and raw reads without normalization or filtering following the DESeq2 vignettes (Love et al., 2014). First, miRNA expression at the low, middle, and high dose KCZ treatment at each timepoint was compared to that in the vehicle control from the same timepoint. miRNAs with fold change (FC) >1.5 or < 0.67 and adjusted p-value < .05 were considered differentially expressed. At each timepoint, a Venn diagram was generated by comparing DEMs at different doses; an area with DEMs associated with a DILI phenotype was filled with a specific color. Lastly, phenotype-specific DEMs were identified by overlapping colored pools from all three timepoints.

miRNA conservation analysis

miRNA conservation analysis was performed based on sequence and annotation information from miRBase (Release 21) http://www.mirbase.org/ (last accessed July 2020). Sequences of phenotype-specific DEMs identified in KCZ-treated rats were compared to those of their human homologs. Comparison results were categorized as identical, when a miRNA sequence shared 100% similarity between rats and humans, or as mismatched, when a miRNA sequence had <80% similarity between rats and humans and contained at least one mismatching nucleotide in the seed region (Artzi et al., 2008). In the present study, DEMs with identical sequences were considered conserved between rats and humans with translatability. The “mismatched” category was not considered conserved or translatable because changes in the nucleotide sequence of a miRNA, especially within the seed region, may affect miRNA target recognition and binding and thus alter a miRNA’s function (Bartel, 2009; Lewis et al., 2005).

Cell culture and KCZ treatment

Terminally differentiated HepaRG cells were purchased from ThermoFisher. Cells were thawed and cultured in Williams’ E medium supplemented with the Thaw, Plate, & General Purpose Medium Supplement and 1 × GlutaMAX (ThermoFisher). Cells were seeded in collagen-coated 96-well plates (3.6 × 10⁴ cells/well) for cytotoxicity assays and collagen-coated 12-well plates (7 × 10⁵ cells/well) for intracellular and extracellular miRNA quantification. Four hours after seeding, cells were fully attached and then treated with 0.1% DMSO (a vehicle for KCZ) or KCZ (Sigma-Aldrich, St. Louis, Missouri) at indicated concentrations for 18 h. Cells were maintained at 37˚C in a humidified tissue culture incubator with 5% CO₂.

Cytotoxicity assays

MTS [3-(4,5-dimethylthiazol-2-yl)-5-(3-carboxymethoxyphenyl)-2-(4-sulfophenyl)-2H-tetrazolium] cell viability assays were performed using CellTiter 96 AQueous One Solution Reagent (Promega Corporation, Madison, Wisconsin). Briefly, after DMSO or KCZ treatment, culture media of HepaRG cells in 96-well plates were replaced with 10 µL of MTS reagent pre-mixed with 90 µL serum-free media per well. The same reagent solution was added to a set of cell-free wells to determine background signals. Following 1-h incubation in a tissue culture incubator, the absorbance at 490 nm was measured using a Synergy 2 Multi-Mode Microplate Reader (BioTek, Winooski, Vermont). Background signals were subtracted from sample signals. Cell viability was calculated by comparing the absorbance of KCZ-treated cells to that of DMSO controls. Three independent experiments were performed with four technical replicates for each treatment condition.

Cellular ATP levels were determined using a CellTiter-Glo Luminescent Cell Viability Assay kit (Promega Corporation) as described previously (Ren et al., 2017). The cellular ATP levels of KCZ-treated cells were calculated as the percentage compared to those of DMSO controls. Three independent experiments were performed with four technical replicates for each treatment condition.

Extracellular and intracellular miRNA quantification

After indicated treatment, HepaRG cells in 12-well plates and culture media from corresponding wells were collected for intracellular and extracellular miRNA quantification, respectively.

Cell culture media were centrifuged at 1000 × g for 5 min at room temperature to pellet dead cells. Supernatants were collected for extraction of total extracellular RNA using a MagMAX mirVana Total RNA Isolation kit (ThermoFisher). The extraction was carried out with a KingFisher Flex Magnetic Particle Processor 96DW (ThermoFisher) according to the manufacturer’s protocol. A synthetic miRNA, ath-miR159a (ThermoFisher) was added in lysis buffer during extraction as a spike-in control. For each sample, 100 µL of supernatants were used to yield 50 µL elute of total RNA. For each miRNA, miRNA reverse transcription was conducted using 5 µL of eluted RNA, TaqMan MicroRNA Reverse Transcription kits, and miRNA-specific TaqMan miRNA assays (ThermoFisher) according to the manufacturer’s instructions.

Intracellular miRNA levels were determined as described previously (Li et al., 2019a). Briefly, miRNeasy Mini kits (Qiagen, Valencia, California) were used to extract total RNA from HepaRG cells. The quality and quantity of extracted intracellular RNA were determined with a Nanodrop NanoDrop 2000 spectrometer (ThermoFisher). miRNA reverse transcription was conducted using 100 ng of the extracted intracellular RNA, TaqMan MicroRNA Reverse Transcription Kits, and TaqMan miRNA assays (ThermoFisher) according to the manufacturer’s protocol.

For both extracellular and intracellular miRNAs, miRNA levels were determined by real-time PCR using TaqMan Universal Master Mix II, no UNG (ThermoFisher) according to the manufacturer’s instructions. Real-time PCR reactions were carried out by a BioRad CFX96 Touch Real-Time PCR Detection System (BioRad, Hercules, CA). The levels of extracellular or intracellular miRNAs were normalized against that of ath-miR159a (a spike-in control) or RNU6 (an endogenous small RNA control) and analyzed using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). Three independent experiments with five replicates in each experiment were performed.

miRNA target gene prediction and conservation

Putative target genes of indicated miRNAs were predicted using miRWalk 3 at http://mirwalk.umm.uni-heidelberg.de/ (last accessed July 2020) via identification of miRNA binding sites within the 5′ untranslated region (UTR), coding sequence (CDS), and 3′ UTR of rat and human mRNA transcripts (Sticht et al., 2018). In rats, potential target genes were selected by meeting these criteria: (1) they possessed miRNA binding sites with a binding probability > 0.95; and (2) they were concurrently predicted by miRWalk 3 and the miRDB algorithm [a filter provided by miRWalk 3 (Sticht et al., 2018)]. These target genes were then visualized in a miRNA/mRNA network using Cytoscape 3.7.2 (https://cytoscape.org/, last accessed July 2020) and subjected to subsequent enrichment analyses.

For a stringent target conservation analysis, potential target genes of indicated miRNAs in rats and humans were selected by meeting these criteria: (1) they possessed miRNA binding sites with a binding probability > 0.99; (2) the predicted miRNA/target binding energy was less than −20 kcal/mol [indicating high binding strengths (Yu et al., 2015)]; and (3) they were concurrently predicted by more than one algorithm (miRWalk 3, miRDB, and TargetScan in humans, and miRWalk 3 and miRDB in rats). Binding probability and energy were calculated by miRWalk 3 for each miRNA target prediction, and miRDB and TargetScan algorithm filters were provided by miRWalk 3 for result comparison and sorting (Sticht et al., 2018). The TargetScan filter in miRWalk 3 showed no target prediction for any rat miRNAs, therefore, was not applied to our target prediction analysis in rats.

Pathway and gene ontology-term enrichment analyses

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and gene ontology (GO)-term enrichment analysis for predicted target genes were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 at https://david.ncifcrf.gov/home.jsp (Huang da et al., 2009a,b; last accessed July 2020). Pathways or terms with a p-value < 0.05 were considered significantly enriched. Graphs were generated in RStudio (version 1.2.5019).

Correlation analysis

Normalized miR-seq and RNA-seq data from the Liver Hepatocellular Carcinoma (LIHC) dataset of the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/, last accessed July 2020) were downloaded from http://firebrowse.org/ (Broad Institute, Boston, Massachusetts, last accessed July 2020). Expression levels of indicated miRNAs and mRNAs in the 50 nontumor human liver samples of LIHC were extracted and used for Pearson’s correlation analysis in GraphPad Prism 5 (GraphPad Software, San Diego, California). Correlations with a p-value < 0.05 were considered statistically significant. Expression levels of the miRNA and mRNA were presented as RPM (reads per million miRNA mapped) and RSEM (RNA-Seq by Expectation-Maximization), respectively.

miRNA:mRNA duplex formation analysis

Formation of miRNA:mRNA duplexes was analyzed using RNAhybrid https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/, last accessed July 2020). Minimum free energy (ΔG) < −20 kcal/mol was set as the threshold. miRNA:mRNA duplexes with at least 6-nt base paring at the seed sequence of miRNAs were considered.

Statistical analyses

Statistical analyses were conducted using GraphPad Prism 5. Data were analyzed using one-way analysis of variance (ANOVA) followed by the Dunnett’s test. The results were displayed as mean with standard deviation from at least three independent experiments. A p-value less than 0.05 was considered statistically significant.

Data availability

miR-seq data generated and analyzed in this study were deposited at the NCBI’s Gene Expression Omnibus (GEO, accession number GSE159200).

RESULTS

The overall study design and workflow are depicted in Figure 1 and each step was detailed in each section below.

Identification of Hepatic DEMs in KCZ-Treated Rats

To identify reliable miRNA biomarkers for KCZ-induced liver injury, we defined two critical criteria for candidate selection: (1) the expression of a miRNA significantly changes in response to KCZ treatment (FC > 1.5 or < 0.67, adjusted p-value < .05); and (2) its differential expression is associated with a DILI phenotype. We first conducted miR-seq to examine hepatic miRNA expression in KCZ-treated rats. Male Sprague Dawley rats were orally administered with vehicle control or KCZ at the low, middle, or high dose per day for 3, 7, and 14 days. Using DESeq2, we compared miRNA expression in KCZ-treated rats to that in the corresponding control at each timepoint (Figure 2A). Among 763 rat mature miRNAs, 0, 79, and 87 miRNAs were differentially expressed (FC > 1.5 or < 0.67, adjusted p-value < .05) at the low, middle, and high dose of KCZ treatment, respectively, on Day 3 (Figure 2A, left). On Day 7, 115, 135, and 121 DEMs were identified at the low, middle, and high dose, respectively (Figure 2A, middle). On Day 14, 0, 66, and 84 miRNAs exhibited significantly altered expression at the low, middle, and high dose, respectively (Figure 2A, right). The number of DEMs increased along with the dose at each timepoint, except for the middle dose on Day 7. Also, under all conditions with DEM identification, we observed more upregulated than downregulated miRNAs, suggesting that KCZ tended to induce, more than suppress, hepatic miRNA expression in rats.

Identification of DEMs Associated With KCZ-Induced Liver Injury in Rats

The DEMs identified in Figure 2A may indicate exposure to KCZ but do not necessarily indicate liver toxicity induced by KCZ. Thus, we next examined traditional liver toxicity indicators in the liver and blood of matching rats used for miR-seq analysis. In Figure 2B, each pie chart represented a phenotype at a treatment condition, with 1/3 of the area representing each rat and colors indicating severity levels. Two histological phenotypes were detected, whereas elevation of serum ALT and AST levels was not detected in the KCZ-treated rats. Slight cytoplasmic vacuolization was observed in the interlobular bile duct of all three rats treated with the high dose KCZ on Days 3 and 7 (Figure 2B, 3/3 green). On Day 14, cytoplasmic vacuolization was slight in all rats at the middle dose (3/3 green) but worsened to moderate in two of the three rats at the high dose (2/3 brown, 1/3 green). Additionally, the two rats with moderate vacuolization on Day 14—high dose exhibited slight centrilobular hypertrophy, a more severe liver injury phenotype than cytoplasmic vacuolization. Under all conditions, no significant ALT and AST elevation occurred. Due to the small sample size (n = 2) for hypertrophy, we focused on vacuolization for the subsequent DEM-phenotype association.

To discover miRNAs that may be associated with KCZ-induced liver injury, we attempted to identify DEMs associated with vacuolization by selecting DEMs in only samples exhibiting the vacuolization phenotype. We first compared DEMs from the different doses at each timepoint. On Day 3, 30 DEMs were identified specifically at the high dose (Figure 2C, blue area), where vacuolization was observed (Figure 2B, left). On Day 7, 18 miRNAs exhibited differential expression only at the high dose (Figure 2C, red area), in association with vacuolization (Figure 2B, middle), but not at the low or middle dose. On Day 14, 57 common DEMs were found at the middle and high dose concurrently (Figure 2C, green area), distinguishing these two conditions with vacuolization from the low dose (Figure 2B, right). By overlapping the three vacuolization-specific pools (blue, red, and green), we identified four common DEMs, namely rno-miR-331-3p, rno-miR15b-3p, rno-miR-34a-5p, and rno-miR-676 (Figure 2D, yellow area). As shown in Figure 3, rno-miR-331-3p, rno-miR-34a-5p, and rno-miR15b-3p were significantly upregulated (FC > 1.5, adjusted p < .05), whereas rno-miR-676 was significantly downregulated under all four treatment conditions with vacuolization (Day 3—high dose, Day 7—high dose, and Day 14—middle and high doses). Differential expression of these four miRNAs was not detected at any conditions where vacuolization was undetectable. Our findings suggest that rno-miR-331-3p, rno-miR15b-3p, rno-miR-34a-5p, and rno-miR-676 may serve as potential biomarkers of KCZ-induced liver toxicity in rats specifically indicating cytoplasmic vacuolization in the bile duct.

Characterization of Phenotype-Specific DEMs for Translatability

To investigate whether these putative miRNA biomarkers can be translated from rats to humans, we analyzed the conservation of the candidate miRNAs based on sequence annotations from miRBase (Release 21). Among the four phenotype-specific DEMs, miR-331-3p, miR15b-3p, and miR-34a-5p displayed 100% sequence identity between rats and humans (Table 1). Sequences of miR-676 in the two species showed 73% similarity, with one of the six mismatching nucleotides (underlined in Table 1) in the seed region (2nd–7th nucleotide from 5′ end, in bold). Altered nucleotide sequences (<80% sequence similarity), especially in the seed region, may affect a miRNA’s target recognition and binding and thereby its function (Artzi et al., 2008; Bartel, 2009). Therefore, we considered miR-331-3p, miR-34a-5p, and miR-15b-3p as conserved candidate miRNA biomarkers.

Extracellular and Intracellular Levels of Newly Identified miRNA Candidates in HepaRG Cell Culture Treated With KCZ

To evaluate the responses of these miRNA candidates to KCZ treatment in human liver cells, we used metabolically competent HepaRG cells for in vitro assays. First, we examined the cytotoxicity in HepaRG cells treated with 50, 75, 100, 125 µM of KCZ compared to DMSO control. Our results showed that both cell viability and cellular ATP levels decreased in a concentration-dependent manner in response to KCZ treatment, with 50 µM and 100 µM causing little and moderate toxicity, respectively (Figures. 4A and 4B). Next, we analyzed the extracellular levels of hsa-miR-34a-5p, hsa-miR-331-3p, hsa-miR-15b-3p, and hsa-miR-676-3p in the culture media of HepaRG cells treated with 50 or 100 µM KCZ and their intracellular counterparts in the corresponding cells. Our data showed that 100 µM KCZ treatment for HepaRG cells caused a significant increase (p < 0.05) in the extracellular levels of hsa-miR-34a-5p, hsa-miR-331-3p, and hsa-miR-15b-3p but not hsa-miR-676-3p, whereas 50 µM KCZ had no significant impact on any miRNA candidates (Figure 4C). Notably, in the corresponding KCZ-treated HepaRG cells, the intracellular level of hsa-miR-15b-3p was significantly decreased by about 50% at both 50 and 100 µM, contrasted to the extracellular level of hsa-miR-15b-3p that was significantly increased by 50% at 100 µM KCZ treatment (Figure 4D). There appeared to be an increase trend in the intracellular levels of hsa-miR-34a-5p, hsa-miR-331-3p, and hsa-miR-676-3p at 100 µM of KCZ treatment, however, the alteration was not statistically significant (Figure 4D).

Prediction and Functional Enrichment of Target Genes for Phenotype-Specific DEMs in Rats

miRNAs exert their biological functions by regulating genes involved in related cellular processes. To further understand the functions of the four phenotype-specific DEMs, we predicted their target genes using miRWalk 3 to identify miRNA binding sites within the 5′ UTR, CDS, and 3′ UTR sequences of rat mRNAs. To reduce false-positive prediction, we selected targeting genes with a miRNA-binding probability greater than 0.95 and concurrently identified by miRWalk 3 and the miRDB algorithm (a filter provided by miRWalk 3). Our analysis identified 185 unique target genes for rno-miR-34a-5p, 89 for rno-miR-331-3p, 16 for rno-miR15b-3p, and 11 for rno-miR-676. A miRNA/mRNA network was visualized using Cytoscape (3.7.2), as shown in Figure 5. Among the four miRNAs, rno-miR-34a-5p was predicted to target the highest number of genes, suggesting that it may be highly involved in gene regulation and potentially diverse biological processes. Among all 298 unique target genes identified, only three were targeted by more than one miRNA, indicating that these four miRNAs may operate via independent gene targeting networks.

To further understand the functions of these target genes, we conducted KEGG pathway and GO-term enrichment analyses using the DAVID database (v6.8). Figure 6 shows the top 10 most significantly affected pathways or GO terms (p < .05), involving signaling transduction, transcription regulation, and protein and enzymatic activities that were essential for cell survival. Our results suggest that the 298 target genes may contribute to enzyme-dependent gene and protein regulation that are essential to cell homeostasis. Such observation reflects a combination of regulatory effects from these four phenotype-specific DEMs.

Conservation Analysis for miRNA/Target Relationships Between Rats and Humans

To further examine the translatability of the three conserved miRNA candidates, we analyzed the conservation of their miRNA/target relationships between humans and rats. We applied stringent cutoffs to target gene prediction using miRWalk 3 for this conservation analysis. Potential target genes of miR-34a-5p, miR-331-3p, and miR-15b-3p in rats and humans were selected by meeting these criteria: (1) the probability for miRNA/target binding was greater than 0.99; (2) the energy for the miRNA/target binding was less than −20 kcal/mol; and (3) target genes were concurrently predicted by more than one algorithms (miRWalk 3, miRDB, and TargetScan in humans, and miRWalk 3 and miRDB in rats). The TargetScan filter in miRWalk 3 showed no target prediction for any rat miRNAs and was therefore not applied to rat target gene prediction. Our analyses showed that 13 and 4 genes targeted by miR-34a-5p and miR-331-3p, respectively, were common between rats and humans (Figure 7A). No target genes of miR-15b-3p in rats met the stringent cutoffs; thus, no comparison was performed for miR-15b-3p target genes between rats and humans.

To analyze whether miR-34a-5p and miR-331-3p may regulate their putative targets concurrently identified in rats and humans, we examined the correlation between the expression levels of miRNAs and those of their corresponding target genes in humans. Expression levels of the two miRNAs and 17 target genes (13 for miR-34a-5p and 4 for miR-331-3p) in 50 nontumor human liver samples were extracted from the miR-seq and RNA-seq datasets of the TCGA-LIHC project. We performed Pearson’s correlation analyses for all 17 pairs of miRNA/mRNA (Table 2, in order of increasing r values for each miRNA). The scatter plot in Figure 7B shows that the expression of hsa-mir-34a correlated inversely and significantly with that of TBL1XR1 [Transducin (β)-like 1 X-linked receptor 1, r = −0.304, p < .05]. This result suggested that hsa-miR-34a-5p may negatively regulate the expression of TBL1XR1.

To assess interaction between hsa-miR-34a-5p and TBL1XR1, we analyzed miRNA:mRNA duplex formation using RNAhybrid. We found that hsa-miR-34a-5p interacted with TBL1XR1 mRNA at multiple sites. Figure 7C shows the targeting site with the highest number of continuous base pairing at the seed region (9 mer) and lowest binding energy (ΔG = −27.6 kcal/mol), suggesting high binding strength in miRNA:mRNA duplex formation.

In summary, our analyses suggested that miR-34a-5p may serve as a reliable and translatable candidate biomarker for KCZ-induced liver injury.

DISCUSSION

In the present study, we established a workflow to identify and characterize candidate miRNA biomarkers (Figure 1) that may provide valuable insights for predictive and computational toxicology studies. We utilized conventional toxicity indicators to direct our identification of candidate miRNA biomarkers for DILI. It depicts a multidimensional perspective of DILI development to associate miRNA expression with DILI phenotypes at various treatment doses and timepoints; this parallel analysis approach can be adopted by other toxicogenomic studies. We also used an in vitro model, drug-treated HepaRG cell culture, to evaluate the translatability of miRNA candidates from rats to humans. For candidate miRNA biomarker characterization, the types of analyses may be adapted depending on the results from biomarker identification. For example, miRNA functional enrichment may be performed when a large number of miRNA candidates have been identified. Association of miRNA and mRNA expression profiles from the same samples may be analyzed to provide mechanistic insights when a sufficiently large sample size is available. Involvement of transcription factors that regulate miRNA expression may also be predicted computationally based on their binding to the promoter regions of miRNA genes. Our computational analyses were guided by toxicological phenotypes and biological mechanisms; therefore, they should improve high-confidence prediction of miRNA candidates for DILI.

Translatability from preclinical models to clinical applications is critical in miRNA biomarker discovery and development. It would be ideal to measure and compare the levels of candidate miRNAs in human and rat samples, preferably serum samples, in response to drug treatment. However, limited sample availability may hinder this straightforward approach. In addition, utilization of in vitro systems is in line with the 3 Rs principle for the welfare of animals in research (Sneddon et al., 2017). Human hepatoma HepaRG cells are metabolically competent with the expression of drug metabolizing enzymes comparable to primary human hepatocytes; thus, they serve as a good surrogate for results extrapolation (Bernasconi et al., 2019; Paini et al., 2017). Furthermore, this in vitro system allows for paired examination of extracellular and intracellular miRNA levels (Figure 4) and is compatible with high-throughput analyses (Leite et al., 2012; Seo et al., 2019).

Intracellular and extracellular miRNAs may respond differently to drug treatment and toxicity. Our results showed that extracellular hsa-miR-15b-3p was upregulated while intracellular hsa-miR-15b-3p was downregulated in HepaRG cell culture exposed to 100 µM KCZ, in association with KCZ-induced cytotoxicity (Figure 4). Previously, similar observations have been reported previously in a mouse model (Wang et al., 2009). In mice with liver injury caused by acetaminophen treatment, several miRNAs, including miR-122, miR-15a, miR-29a, and miR-192, showed higher levels in plasma and lower in liver tissues; conversely, miR-710 and miR-711 exhibited an elevation in the liver but a decrease in plasma (Wang et al., 2009). Importantly, the downregulation of intracellular miR-15b-3p was unlikely caused by reduced RNA extraction yield due to cell death, because the same amount of RNA (100 ng) were used in miRNA reverse transcription for all the samples. Meanwhile, dead cells were removed from cell culture media before extracellular RNA extraction to prevent elevation in extracellular miRNA levels caused by lysis of dead cells. Finally, the discrepancy in miR-15b-3p dysregulation between HepaRG cells (Figure 4D) and rat liver tissues (Figure 3) may be attributed to heterogenous components in liver tissues, compared to a single HepaRG cell type. For example, bile ducts and blood vessels in liver tissues may also carry bile and blood, respectively, in which extracellular miR-15b-3p may be present and elevated upon KCZ treatment.

We have successfully demonstrated the feasibility of our integrated strategy in experimentally characterizing miRNA functions in several previous studies; the in silico analyses and cutoffs used for miRNA:mRNA targeting applied in this study were shown to result in high confidence validation in human liver cells (Li et al., 2019a, 2020; Yu et al., 2018). Meanwhile, others have reported additional functions in the liver for our candidate miRNA biomarkers. For example, miR-34a-5p was upregulated in the liver of rats exposed to aflatoxin B1 (AFB1) and might serve as a serum biomarker for AFB1-induced liver injury (Liu et al., 2015). miR-34a was also overexpressed in steatosis-induced hepatocytes and influenced lipid metabolism by suppressing peroxisome proliferator-activated receptor alpha in mice in nonalcoholic fatty liver disease (Ding et al., 2015). Also, the level of serum miR-331-3p was elevated in patients with early stage hepatocellular carcinoma (Sun et al., 2020). We focused on designing a computational strategy for miRNA biomarker screening for DILI here, and functional characterization for candidate miRNA biomarkers will be conducted in future studies.

In addition to the discovery of useful miRNA biomarkers for DILI, we endeavored to identify miRNA/mRNA target pairs that are conserved between preclinical animals and humans. Interaction between miRNAs and their putative mRNA target may be more biologically relevant than miRNAs alone, as miRNAs exert cellular functions though their target genes. However, it is worth noting that the expression of protein coding genes is modulated by a multiplex network including different types of cellular factors. First, one gene may be regulated by multiple miRNAs and one miRNA may regulate multiple genes (Shivdasani, 2006; Shu et al., 2017; Xu et al., 2020). Second, gene expression is often simultaneously regulated by other factors besides miRNAs, such as transcriptional regulation by transcription factors and genetic variations, multilevel modulation by long noncoding RNAs, and epigenetic regulation including DNA methylation and histone modification (Barter et al., 2012; Bianchi et al., 2017; Humphries et al., 2019; Li et al., 2019b). Third, miRNAs tend to exert moderate regulatory effects on gene expression, which may not be detected under all conditions due to technical sensitivity or biological override by other factors (Bartel, 2009; Schira-Heinen et al., 2020; Specjalski and Jassem, 2019). Fourth, gene regulation is dependent on distribution and local miRNA concentration in the cells, which is difficult to determine experimentally (Ning et al., 2019). The candidate miRNA biomarkers and the miRNA/mRNA pair identified in this study require experimental validation and mechanistic elucidation. Furthermore, extracellular levels of candidate miRNAs may be measured in the blood of treated laboratory animals to indicate their suitability for noninvasive detection.

Both sensitivity and specificity are critical factors in assessing the suitability of biomarkers. In this study, we were limited by a small sample size, which is typical for current toxicogenomics studies, but insufficient for generating precise estimations with narrow confidence intervals that are informative in clinical use (Hajian-Tilaki, 2014). Samples size required for statistically meaningful estimation of specificity and sensitivity depends on many factors related to study design and predetermined values in clinical applications; for example, a sample size of 1254 will be required for a sensitivity of 80% with the prevalence of a disease being 10%, a marginal error less than 7%, and a 95% confidence level (Hajian-Tilaki, 2014).

U.S. FDA employs stringent rules and a multistep submission process to qualify biomarkers for regulatory use https://www.fda.gov/drugs/cder-biomarker-qualification-program/about-biomarkers-and-qualification (last accessed September 2020). Development of miRNA biomarkers for DILI will require extensive investigation. In this study, our goal was to establish a strategy for accurate prediction of miRNA biomarker candidates and to provide preliminary insights on a few candidates for KCZ-induced liver injury.

CONCLUSIONS

In summary, we profiled global miRNA expression in the livers of rats treated with KCZ at multiple doses and durations. Using a series of computational analyses, we identified four potential candidate miRNA biomarkers: miR-34a-5p, miR-331-3p, miR-15b-3p, and miR-676. Significant upregulation of rno-miR-34a-5p, rno-miR-331-3p, and rno-miR-15b-3p and downregulation of rno-miR-676 were associated with cytoplasmic vacuolization, a DILI phenotype in rat livers exposed to KCZ. Such dysregulation preceded AST and ALT elevation in serum, indicating the potential of using these miRNA candidates for early detection of KCZ-induced liver injury. Evolutionarily conserved between rats and humans, miR-34a-5p, miR-331-3p, and miR-15b-3p showed an increase in extracellular levels in HepaRG cell culture exposed to KCZ, at a concentration that induced cytotoxicity. miR-676 showed 73% sequence similarity between the two species and nonsignificant changes in HepaRG cell culture in response to KCZ treatment. Additionally, a miRNA/target pair, miR-34a-5p/TBL1XR1, was predicted in both rats and humans with high binding strengths and exhibited inverse correlations of expression in nontumorous human liver samples. The toxicogenomics workflow we created will provide valuable insights for predictive and computational toxicology studies.

AUTHOR CONTRIBUTIONS

Conceptualization: D.L., W.T., and B.N.; project management: D.L. and B.N.; data acquisition: D.L., B.K., S.C., and L.G.; data analysis: D.L., B.G., and Z.L.; writing—original draft preparation: D.L.; writing—review and editing: W.T. and B.N. All authors have read and agreed to the published version of the manuscript.

FUNDING

This study was funded by the National Center for Toxicological Research/U.S. Food and Drug Administration, Project E0753201.

DECLARATION OF CONFLICTING INTERESTS

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Disclaimer: This article is not an official guidance or policy statement of the U.S. Food and Drug Administration (FDA). No official support or endorsement by the U.S. FDA is intended or should be inferred.

REFERENCES