An in silico argument for mitochondrial microRNA as a determinant of primary non function in liver transplantation

Mitochondria have their own genomic, transcriptomic and proteomic machinery but are unable to be autonomous, needing both nuclear and mitochondrial genomes. The aim of this work was to use computational biology to explore the involvement of Mitochondrial microRNAs (MitomiRs) and their interactions with the mitochondrial proteome in a clinical model of primary non function (PNF) of the donor after cardiac death (DCD) liver. Archival array data on the differential expression of miRNA in DCD PNF was re-analyzed using a number of publically available computational algorithms. 10 MitomiRs were identified of importance in DCD PNF, 7 with predicted interaction of their seed sequence with the mitochondrial transcriptome that included both coding, and non coding areas of the hypervariability region 1 (HVR1) and control region. Considering miRNA regulation of the nuclear encoded mitochondrial proteome, 7 hypothetical small proteins were identified with homolog function that ranged from co-factor for formation of ATP Synthase, REDOX balance and an importin/exportin protein. In silico, unconventional seed interactions, both non canonical and alternative seed sites, appear to be of greater importance in MitomiR regulation of the mitochondrial genome. Additionally, a number of novel small proteins of relevance in transplantation have been identified which need further characterization.

. 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.95 6 . 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 to 7-9 .
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 MitoWheel 10 , 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 editor 13 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 genome 14 . 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 wobbles 15, 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 points 17, 18 . Both rna22 and MiRWalk are based on the known cystolic rules of miRNA seed behaviour.
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 prediction 20,21 .
The miRDB generated list of predicted gene targets for a given PNF miRNA were then exported into MitoMiner4.0 v2016 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 data 2,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 localization 22, 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 data 24 , and KEGG, to explore protein biological pathways 25,26 . If the PNF miRNA target was an open reading frame (ORF), the nucleotide sequence was translated in UniProt 24 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 spectrometry 27,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 identity 29,30 . For documented tissue and cell expression the Protein Atlas was searched 31,32 . For orthology exploration the HGNC Comparison of Orthology Predictions (HCOP) was used 33 . Protein secondary structure was explored with RaptorX 34,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 predictions 36,37 .

Results
Liver PNF microRNA. Based on archival miRNA array data 6 , 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).
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).
Based on homology, C5orf63 protein is a small glutaredoxin-like protein with some structural alignment to a chaperone protein 38 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 stress 39 .
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 matrix 40 . 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 Ca 2+ ions to promote autophagy 41,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 analysis 32 . 3D modelling predictions for C19orf12 protein (isoform 4) give structural alignment predictions that identify it be a putative importin/exportin 43,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 outcome 6 . 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 49 . 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 structure 50 . 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). 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, TMscore has the value in 0 to 1, where 1 indicates a perfect match between two structures.  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 cytosol 51 . 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 respectively 52,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 transcription 56 . Whereas nuclear generated MitomiRs traffic between all compartments of the cell and human body to act as translational repressors 57,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 transplantation 59 .
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 cytosol 16,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′UTR 51,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 model 59 .
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 enzymes 1 . 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 stress 39 .
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 generation 40,61 . Graft recovery and viability is intimately involved with the ability of mitochondria to generate high energy phosphates 62 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 stress 41,42 . Structural based assignment identify the protein to be an importin/exportin chaperone 43,44,63 which could be involved in the control of nucleic acid traffic between the nuclear, cytosolic and mitochondrial compartments to define cell viability 64,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 experiments 6 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.