Introduction

Cold preservation has been the mainstay of donor liver maintenance in transplantation for the past 50 years. However, a period of cold ischemia followed by reperfusion produces an ischemia reperfusion injury that is the main determinant of graft function and outcome1. To ameliorate ischemia reperfusion injury and the risk of primary non function (PNF) or dysfunction, cold ischemic times are ideally kept to a minimum. Relatively little is known about mitochondrial behaviour on reperfusion, but it is recognized that the integrity and resilience of mitochondria on reperfusion is a key determinant of the extent of injury and outcome.

Mitochondria are central to cellular energy generation and sit at the crossroads of a number of critical cell pathways e.g glycolysis/OXPHOS balance, senescence and apoptosis. Mitochondria have their own genomic, transcriptomic and proteomic machinery but are unable to be autonomous, needing both the nuclear and mitochondrial genomes to coordinate the expression of at least 2000 mass spectrometry identified mitochondrial proteins, of which 13 are encoded by the mitochondrial genome (mt genome). It is estimated that only 3% of the mitochondrial proteome (mt proteome) is involved in ATP production2,3.

Mitochondrial microRNA (MitomiRs) are a microRNA (miRNA) species that have been characterized to localize to mitochondria. The majority of MitomiRs originate from the nuclear genome and are imported into the organelle4,5. MitomiR function is presently being determined but it is thought to encompass metabolic control through their ability to adjust the composition of the mitochondrial proteome and potentially, provide a language for anterograde retrograde nuclear mitochondrial communication. Our previous work on miRNA expression in PNF demonstrated miR-22 to be of importance and with predicted influence on metabolic and apoptotic pathways, both of which centre on mitochondria6. This observation and the continuing paucity of data on PNF has motivated the present re-analysis. Using computational biology we have explored MitomiRs and miRNA interactions with the mitochondrial proteome in the clinical transplant model of PNF in the donor after cardiac death (DCD) liver.

Materials and Methods

miRNA Array source data and overview of computational biological strategy

Archival data, from our previous work that identified 16 miRNA to have differential expression in DCD PNF6, was re-analyzed using a number of publically available computational algorithms. Figure 1 summarizes computational algorithms, databases and strategy applied. The archival data analysed in this present manuscript originated from RNA that had been isolated from post perfusion formalin fixed paraffin embedded (FFPE) trucut liver biopsies using High Pure FFPE RNA micro kit (Roche Diagnostics Ltd, Hertfordshire, UK). The three DCD groups compared were 1. PNF retransplanted within a week (n = 7), 2. good functional outcome (n = 7) peak aspartate transaminase (AST) ≤ 1000 IU/L and 3. early graft dysfunction (n = 9) peak AST ≥ 2500 IU/L. High throughput expression analysis was performed using Affymetrix GeneChip miRNA 2.0 Arrays according to manufacturer recommendations.

Figure 1
figure 1

Computational work stream summarizing in silico analysis of mitochondrial microRNAs (MitomiRs) interactions with the mitochondrial and nuclear genomes. Archival miRNA array data from the donor after cardiac death (DCD) liver was filtered to identify MitomiRs. MitomiR seed interactions with the mitochondrial genome were explored by applying simple canonical seed rules to the both coding and non coding areas. Additional algorithms applied were rna22 and miRWalk 2.0 that are based on the principles of cystolic miRNA interactions. The nuclear encoded genome/MitomiR interaction was explored with miRDB. Gene targets were filtered with MitoMiner4.0 for mitochondrial proteins with an Integrated Mitochondrial Protein Index (IMPI) > 0.7. Identified mitochondrial protein was then characterized in Uniprot and Kyoto Encyclopedia of Genes and Genomes (KEGG) for biological function. If an open reading frame (ORF) protein was identified, putative biological role was explored in the Protein Data Base (PDB), Protein Atlas and HGNC Comparison of Orthology Predictions (HCOP). 2D and 3D modelling of ORF protein was performed using Raptor-X and I-Tasser respectively, to further characterize ORF protein function. URL shortcut to databases and algorithms used are listed.

Initial data analysis and mining was performed using Qlucore Omics Explorer (QOEv2.1) version 2.1. Statistical significance level was set at p = 0.05 with a false discovery rate q value of 0.956. The DCD PNF miRNA selected for further in silico analysis for this present manuscript are summarized in Table 1. Only miRNA with human homology has been considered. miRBase v21 published miRNA sequences and annotation have been adhered to7,8,9.

Table 1 Summary of Primary Non Function microRNA selected for in silico analysis. Mature microRNA (miRNA) identifiers have been used from miRBase. Table summarizes miRNA array fold change, accession number and sequence, total number of miRNA gene targets as identified by miRDB v21, total number of miRNA gene targets recorded in MitoMiner4.0 and number of MitoMiner4.0 identified genes with an Integrated Mitochondrial Protein Index (IMPI) score >0.7. miRNA that have been experimentally isolated from mitochondria is marked “yes”(Y) to identify as a MitomiR.

Liver PNF MitomiRs

To identify the miRNA that are MitomiRs i.e documented as being isolated from mitochondria, published next generation sequencing (NGS), miRNA array, RNA cross-linking immunoprecipitation (CLIP), RT-qPCR or functional data was cross referenced. The interactions of MitomiRs with the mt genome were visualized using MitoWheel10, a graphical summary of the circular human mt genome (16569 base pairs) that is based on the standard revised Cambridge reference sequence (GenBank NCBI database accession number NC_012920)11,12.

Initial exploration of MitomiR interactions with the mt genome used simple canonical binding rules of a 6 seed sequence, as identified from nucleotide 2 to 7 from the 5′ prime end of a given miRNA. An online sequence editor13 was used to interconvert between RNA and DNA to allow for consensus matching of miRNA seed sequence to the input sequence of the mt genome. See Fig. 2 for map summary of the human mt genome14. Further in silico analysis of MitomiR interactions with the mt genome were performed using rna22 and MiRWalk 2.0 both at default settings. rna22 uses the Teiresias algorithm to predict miRNA:mRNA heteroduplexes allowing for both bulges and G:U wobbles15,16. While MiRWalk 2.0 uses predicted and validated information on miRNA-target data to query a given sequence for identification of miRNA species and their binding points17,18. Both rna22 and MiRWalk are based on the known cystolic rules of miRNA seed behaviour.

Figure 2
figure 2

Map of the Mitochondrial Genome. The H (heavy) strand is the outer circle and L (light) strand is the inner circle. There are 22 transfer RNA (TRN) genes for the following amino acids: F, V, L1 (codon UUA/G), I, Q, M, W, A, N, C, Y, S1 (UCN), D, K, G, R, H, S2 (AGC/U), L2 (CUN), E, T and P, which are shown as white boxes. There are 2 ribosomal RNA (RRN) genes: S (small subunit, or 12S) and L (large subunit, or 16S) shown as blue boxes. There are 13 protein-coding genes: 7 for NADH dehydrogenase subunits (ND, yellow boxes), 3 for cytochrome c oxidase subunits (COX, orange boxes), 2 for ATPase subunits (ATP, red boxes), and one for cytochrome b (CYTB, coral box). Two gene overlaps are indicated (ATP8-ATP6, and ND4L-ND4, black boxes). The control region (CR) is the longest non coding sequence (grey box) and its three hypervariable regions are indicated (HV, green boxes). The most polymorphic area is HVRI (16024–16383) followed by HVRII (57–372) and HVRIII (438–574). Numbering is from the Cambridge Reference Sequence (van Oven 2009). Map of the human mitochondrial genome courtesy of Emmanuel Douzery (14).

Liver PNF miRNA and the mitochondrial proteome

In order to identify interactions of a given PNF miRNA with the mt proteome, the mature miRNA identifier was used for data mining of predicted gene targets using miRDB (MicroRNA Target Prediction And Functional Study Database)19. miRDB is an online database that uses a miRNA target prediction program based on a support vector machine (supervised machine learning). For target mining, threshold values were kept at the default producing generated target thresholds between 50–100, the higher the score the more confidence in a given prediction20,21.

The miRDB generated list of predicted gene targets for a given PNF miRNA were then exported into MitoMiner4.0v2016 APR(last accessed July 2017)22,23. MitoMiner4.0 uses data from the MitoCarta2.0 inventory of genes that encode for mitochondrial proteins originally derived from mitochondrial mass spectrometry data2,3. MitoMiner4.0 further integrates this data with data from Universal Protein Resource (UniProt), Gene Ontology, Online Mendelian Inheritance in Man, HomoloGene, Kyoto Encyclopaedia of Genes and Genomes (KEGG) and PubMed to calculate the Integrated Mitochondrial Protein Index (IMPI) score. The MitoMiner4.0 algorithm uses machine learning to evaluate these multiple data sources with a random forest learning classifier. Genes with an IMPI score of 0.7 or greater, being regarded as strong evidence for gene product mitochondrial localization22,23.

The mitochondrial proteins (IMPI > 0.7), computationally identified to be regulated by a given PNF miRNA were characterized in more detail using UniProt for protein sequence and annotation data24, and KEGG, to explore protein biological pathways25,26. If the PNF miRNA target was an open reading frame (ORF), the nucleotide sequence was translated in UniProt24 and the FASTA amino acid sequence downloaded for protein modelling and prediction. An ORF encodes for a hypothetical protein that presently does not have a characterized homologue in the protein database but has some data demonstrating gene product expression from microarray and mass spectrometry27,28. To explore ORF protein features and putative function, a search was performed via the Protein Data Bank (PDB) portal based on FASTA sequence using BLAST with the following limits of mask low complexity, e value cut off 10 and for retrieval of homologous proteins at 70% sequence identity29,30. For documented tissue and cell expression the Protein Atlas was searched31,32. For orthology exploration the HGNC Comparison of Orthology Predictions (HCOP) was used33. Protein secondary structure was explored with RaptorX34,35 and the dictionary of protein secondary structure nomenclature applied. I-Tasser (Iterative Threading ASSEmbly Refinement) that uses a hierarchial approach was used for further 3D protein structure modelling, ligand binding and enzymatic site predictions36,37.

Results

Liver PNF microRNA

Based on archival miRNA array data6, 16 miRNA were identified as potential distinguishers for DCD liver PNF with widespread interactions with the nuclear genome as predicted by miRDB. All of the identified DCD PNF miRNA were downregulated in expression on miRNA array analysis. Of these 16 miRNA, 10 have been identified to be MitomiRs based on published data that has demonstrated their presence in mitochondria. Whereas, MitoMiner4.0, characterized 6 PNF miRNA to be targeting transcripts that encoded for mitochondrial proteins IMPI > 0.7 (Table 1).

Considering PNF MitomiRs alone, seven had no direct regulatory action on the nuclear derived mt proteome as identified by MitoMiner4.0 (miR-378, miR-122-5p, miR-125b-5p, miR-24, let-7a, let-7d-5p, miR-22-3p). Leaving 3 MitomiRs (miR-107, miR-23b, miR-103a-3p) with a MitoMiner4.0 identified interaction with the nuclear encoded mt proteome. One PNF miRNA (miR-191-5p) was not recognized as either being a MitomiR or having any regulatory interaction with the nuclear encoded mt proteome (Table 1).

Predicted Interactions of MitomiRs with the Mitochondrial genome

Three MitomiRs miR-107, miR-103a-3p and miR-22-3p had no identified interactions with the mt genome based on canonical seed binding predictions (Table 2). However, miR-107 and miR-103a-3p had interactions with the mitochondrial encoded nuclear genome based on target mitochondrial protein having an IMPI > 0.7 (Table 1). The remaining 7 MitomiRs (miR-378, miR-23b, miR-122-5p, miR-125b-5p, miR-24, let-7a, let-7d-5p) had multiple potential sites of interaction with the mt genome that included both coding and non coding areas (Table 2). Non coding areas for interaction included the Hypervariable Region 1 (HVRI) for miR-23b and miR-122-5p; with miR-23b, let-7a and let-7d-5p additionally interacting with the control region. The majority of the MitomiR targets were on the light strand (L-strand) (Table 2).

Table 2 MitomiR seed predictions for interactions with transcripts from the mitochondrial genome. Table summarizes total number of seed binding sites based on canonical binding, area of mitochondrial genome affected, its functional area and the transcript affected. Small “r” after nucleotide number means reverse and is synonymous with the heavy strand. Non coding areas are Hypervariable Region 1 (HVR1) and the Control Region (CR).

Repeating the above analysis using open source computational algorithms which are based on the principles of cystolic miRNA interactions. rna22, did not identify MitomiRs miR-107, miR-378, miR-103a-3p and 125b-5p to have an interaction with the mt genome. Otherwise, the remaining MitomiRs were characterized to interact with coding areas of the mt genome. Only miR-122-5p targeted non coding areas that included the HVR1 (see Supplementary Table S1). While MiRWalk did not recognize any of the PNF MitomiRs to be interacting with the mt genome (see Supplementary Table S2). For the remainder of this in silico analysis, the initial canonical MitomiR interactions with the mt genome based on conventional seed binding as documented in Table 2 has been applied as the rules of MitomiR binding in the matrix are not established.

Predicted PNF microRNA mitochondrial protein regulation

Considering PNF miRNA regulation of the mt proteome, 7 Open reading frame (ORF) proteins or conserved hypothetical proteins were identified with an IMPI > 0.7. The mitochondrial ORF proteins identified were small, ranging from 113 to 616 amino acids. Orthology predictions demonstrated them to be highly conserved. All appear to localize to mitochondria, with functions that range from cofactor for ATP Synthase assembly (C7orf55), REDOX balance (C5orf63, C2orf69, C8orf82), importin/exportin ankyrin motif protein (C19orf12) or remain as an uncharacterized conserved mitochondrial protein (C15orf40, C14orf159) (Table 3 and Fig. 3). The interaction of PNF miRNA with the mt proteome appears to be with little redundancy. In that one miRNA typically targets one ORF (Table 3). However, both miR-107 and miR-103 target the C7orf55 transcript, while miR-23b has three transcript targets in the form of C5orf63, C2orf69 and C15orf40, two of which are involved in REDOX (C5orf63, C2orf69).

Table 3 Primary Non Function (PNF) miRNA interactions with the mitochondrial proteome. The table summarizes PNF miRNA interactions with identified transcripts of open reading frames (ORF), Ensembl gene identification (ENSG), location of ORF, number of nucleotides (nt), homolog protein product of ORF and ontological biological process associated with protein as identified by the Protein Database. miRDB target score range 50–100, the higher the score the more likely that a given miRNA targets given gene. Further details of protein class, protein evidence and subcellular localization from Protein Atlas are also associated with identified ORF protein. If ORF protein is recognized by Mitocarta2.0 to be a mitochondrial protein (IMPI > 0.7) it is marked as “yes” (Y) and any ascribed function is summarized.
Figure 3
figure 3

Summary of nuclear encoded mitochondrial ORF proteins targeted by Primary Non Function (PNF) miRNA. Characteristics of the ORF protein include Ensembl gene identification (ENSG), amino acid (aa) length, weight in daltons (Da) and FASTA sequence of translated protein from ORF. In all cases isoform 1 has been selected, apart from C19orf12 where isoform 4 has been used for protein modelling. Predicted secondary structure is documented using the following nomenclature: alpha helix/4 turn helix (H), extended strand in parallel and/or anti-parallel β-sheet conformation (E) and coil (C). 3D model was generated using I-TASSER, estimated global accuracy of model is in −5 to 2, C score > −1.5 indicative of a model with correct global topology. Snap shot image of 3D model generated from I-TASSER has been included. Ligand binding and enzymatic activity protein prediction data have been generated using I-TASSER, the higher the C score, the greater the confidence of the prediction, range 0–1. Structural alignment predictions of function derived from the Protein Data Bank, TM-score is a metric for measuring the structural similarity of two protein models, TM-score has the value in 0 to 1, where 1 indicates a perfect match between two structures.

Based on homology, C5orf63 protein is a small glutaredoxin-like protein with some structural alignment to a chaperone protein38 with predicted ligand binding to glutathione, FAD and heme-s consistent with a role in REDOX (Fig. 3). The evidence for C2orf69 protein a hypothetical oxidoreductase is stronger, with supporting data at protein level and in MitoCarta2.0 (Table 3). Both of these REDOX ORF proteins are identified to have mitochondrial localization. Another mitochondrial ORF protein with a role in REDOX, is C8orf82 protein that is regulated by miR-296-5p, it has both protein and Mitocarta2.0 supporting evidence. From homolog functional data it has pleckstrin homology for nucleotide ligand binding with a proposed role for signal transduction during oxidative stress39.

C7orf55 protein, also known as Formation of Mitochondrial Complex 5 Assembly Factor 1 homolog (FMC1), has a documented role in mitochondrial biogenesis according to MitoCarta2.0 and supporting protein expression evidence. Both miR-107 and miR-103a-3p appear to be involved in FMC1 translational control. Putatively, FMC1 is a matrix localizing protein that is needed for Complex 5 (ATP Synthase) stability and organization at high temperatures. When there is no FMC1 available at times of cell stress, the subunits of ATP Synthase are unable to make a functional oligomer and aggregate in the matrix40. FMC1 has been identified to belonging to the family of proteins carrying the LYR motif (pfam13233: Complex1_LYR_2) that contributes to Fe-S cluster biogenesis (Table 3 and Fig. 3).

Out of the 7 PNF ORF proteins identified, C19orf12 protein appears to be the most well studied but within the clinical context of neurological disease. It has 9 splice versions generating a protein of variable size (range 68–152 amino acids). The protein has been identified to be a component of the endoplasmic reticulum and the mitochondrial associated membrane. The C19orf12 protein co-localizes in mitochondria under oxidative stress, where it supports the transfer of lipids and Ca2+ ions to promote autophagy41,42. However, it is not recognized to be a mitochondrial protein by MitoCarta2.0 but there is some reported evidence for mitochondrial protein expression based on tissue analysis32. 3D modelling predictions for C19orf12 protein (isoform 4) give structural alignment predictions that identify it be a putative importin/exportin43,44. This importin/expotin has predicted ability to bind basic residues such as nucleic acids (Fig. 3).

Discussion

There is a lack of data on predictive markers in DCD liver transplantation and this has led to the re-analysis of archival miRNA array data, despite its limitations and prior demonstration that miRNA is a poor discriminator of functional outcome6. However, by focusing on mitochondrial biology, this in silico work presents an argument for how PNF miRNA interacts with the mt genome and nuclear encoded mt proteome to form the mitochondrial events that culminate in PNF. Using this strategy a number of novel/hypothetical ORF proteins have been identified, as well as giving insight into MitomiR behaviour in the context of PNF. Out of the 16 PNF miRNAs identified, 10 can be classified as MitomiRs as they have previously isolated from mitochondria (Table 1). Four of these MitomiRs have been functionally characterized (miR-378, miR-24, miR-23b and let-7a) with roles in control of hypoxia induced apoptosis45, systemic energy homeostasis46, oxidative capacity47, insulin signalling48 and ROS generation49. All of these processes are of importance for good graft function and survival in liver transplantation.

The mt genome is compact, with little or no intergenic non coding sequences. Additionally, mitochondrial messenger RNA lack a non translated leader and trailing sequence. Other unique features of the mt genome are the control region and the hypervariability region (HVR). The control region contains the point where transition from RNA to DNA occurs, thereby being the point of control for both replication of the mt genome and mitochondrial transcription (Fig. 1). While HVR, despite its name, encodes for a RNA transcript that has a conserved secondary structure50. These features of the mt genome, combined with canonical seed predictions, identify that the majority of the PNF MitomiRs are interacting with coding areas of the mt genome. While non coding area interactions are limited to miR-23b and miR-122-5p for the HVR1, and miR-23b, let-7a and let-7d-5p interacting with the control region. Additionally, the majority of the MitomiR interactions in the cellular context of PNF appear to be directed against the L-strand, suggesting that in PNF there is a loss of mitochondrial polycistronic balance (see Table 2).

In silico, MitomiRs miR-107, miR-103a-3p and miR-22-3p do not have any identified interactions with the transcribed mt genome. Suggesting that unconventional seed interactions, both non canonical and alternative miRNA seed sites are of greater importance in the mitochondrial matrix than the accepted miRNA regulatory events of the cytosol51. Some of these atypical seed rules have already been described for miR-107 and 103a-3p in the cytosol, which target mRNA ORFs and the 5′ UTR respectively52,53.

The majority of MitomiRs are derived from the nuclear genome with a few from the mt genome (miR-4485-3p, miR-1973, miR-4284)4,54,55. There is also some evidence that small noncoding RNA (20-40 nucleotides) derived from the mt genome stay within mitochondria and promote transcription56. Whereas nuclear generated MitomiRs traffic between all compartments of the cell and human body to act as translational repressors57,58. In silico, the PNF MitomiRs (miR-107, miR-103a-3p, miR-22-3p) that do not have any identified interactions with the mt genome based on canonical seed binding rules, do however interact with the nuclear encoded mitochondrial proteome. Potentially, this distribution and differential activity of MitomiRs in the nucleus and mitochondria could provide a mechanism of anterograde retrograde nuclear mitochondrial communications that contribute to mitochondrial dysfunction and failure in this clinical model of liver transplantation59.

Most miRNA computational algorithms are a combination of validated knowledge bases and predicted canonical base pairing of the miRNA seed sequence with mRNA in the context of the cytosol16,18,51. Nevertheless, canonical predictions limited to the 3′UTR can be an over simplification of events as RNA binding protein immunoprecipitation studies have demonstrated that miRNA can bind both non canonically and to sites outside the 3′UTR. In plants, it is established that miRNA binding to coding regions is more common, with binding to the 5′UTR (promoter and enhancer area of a gene) leading to up regulation. Whereas in mammals, miRNA has mainly been characterized to bind to the 3′UTR51,60. Additionally, the mt genome is highly compact. Therefore, MitomiR seed binding rules and behaviour in the matrix are likely to be different to that of the cytosol. MitomiR binding events and regulatory effect remain unknown and need further characterization, making it presently difficult to categorically state in silico, whether a given MitomiR will be a translational repressor or enhancer in the present model59.

Mitochondria are a major source of ROS and maintenance of REDOX balance is essential, in order to protect from oxidative damage, and its cellular sequelae. While previous work has associated poor outcome in liver transplantation with low levels of oxidant enzymes1. The mitochondrial ORF proteins, whose translation in silico that have been identified to be regulated by PNF miRNA (Table 3), have predicted functional roles that appear to be biologically plausible and relevant to the clinical model of PNF in liver transplantation. C5orf63 and C2orf69 proteins are both characterized to be oxidoreductases. While C8orf82 protein act as a signal transducer during oxidative stress39.

C7orf55 protein, or FMC1 identified homolog is involved in the formation and stability of ATP synthase and thereby bioenergetic mitochondrial efficiency. ATP synthase is the concluding protein complex of OXPHOS and is needed to drive the proton gradient for ATP generation40,61. Graft recovery and viability is intimately involved with the ability of mitochondria to generate high energy phosphates62 and an inadequate amount of FMC1 protein would be predicted to lead to ATP synthase failure. In contrast, C19orf12 protein, with an ankyrin motif, appears to act as a shuttle between the cytosol and mitochondria during oxidative stress41,42. Structural based assignment identify the protein to be an importin/exportin chaperone43,44,63 which could be involved in the control of nucleic acid traffic between the nuclear, cytosolic and mitochondrial compartments to define cell viability64,65.

There is presently no biological data on the in silico identified protein targets and the MitomiRs computationally characterized are yet to be associated with their gene targets. The original experiments6 were based on the DCD liver, which is subjected to a warm ischemic period that is one of the determinants of this graft’s unpredictable behavior on reperfusion. However, in order to validate and understand whether the identified PNF MitomiRs and ORF proteins are biologically true and relevant to liver transplantation, our future work will focus on their expression profile in different donor liver types (neonatal, living related, DCD and Donor after Brainstem Death), combined with in vitro experiments to characterize their interactions with the mitochondrial proteome.

In conclusion, the behaviour of MitomiRs within the mitochondrial matrix are ill understood and more study is needed. However, the emerging data on their compartmentalized functional roles could explain, how a cell sets and adjusts its mitochondrial proteome according to cellular energy needs and stresses. Additionally, a number of small mitochondrial ORF proteins that could be of relevance in liver transplantation haven been identified. Despite the limitations of this computational work, validation of these in silico observations on MitomiRs and ORF proteins could provide a way to understand the limits of a DCD liver in transplantation.