|Home | About | Journals | Submit | Contact Us | Français|
MicroRNAs are short, non-coding RNA sequences that regulate genes at the post-transcriptional level and have been shown to be important in development, tissue differentiation, and disease. Limited attention has been given to the natural variation in miRNA expression across genetically diverse populations even though it is well established that genetic polymorphisms can have a profound effect on mRNA levels. Expression level of 577 miRNAs in the livers of 70 strains of inbred mice was assessed, and we found that miRNA expression is highly stable across different strains. Globally, the expression of miRNA target transcripts does not correlate with miRNA expression, primarily due to the low variance of miRNA but high variance of mRNA expression across strains. Our results show that there is little genetic effect on the baseline miRNA levels in murine liver. The stability of mouse liver miRNA expression in a genetically diverse population suggests that treatment-induced disruptions in liver miRNA expression, a phenomenon established for a large number of toxicants, may indicate an important mechanism for the disturbance of normal liver function, and may prove to be a useful genetic background-independent biomarker of toxicant effect.
Microribonucleic acids (miRNA) are short non-coding RNAs that bind to complementary sequences on target transcripts (usually 3′ UTRs), thereby regulating translation. Translation is inhibited either directly or indirectly by nucleolytic degradation of the mRNA transcript . While the mechanistic details underlying miRNA-mediated translational silencing are not fully understood, several key details are known. In plants, miRNAs bind with complimentary target sites which induce transcriptional degradation . In mammals, miRNAs bind to target sites with imperfect complementarity and inhibit elongation of the polypeptide chain during translation . Transcript degradation in response to miRNA delivery has also been reported in mammals . In Drosophila, both translational repression and transcript degradation have been observed, with some genes showing concordant decreases in both protein and mRNA levels and some genes showing decreases in protein levels without a corresponding decrease in mRNA levels . These different modes of action by miRNAs upon their targets have made computational prediction of miRNA target genes difficult, with accuracy rates (correctly predicted/total predicted) of only 50% .
There are over 10,000 catalogued miRNAs in various species, with 579 reported in the mouse (miRBase, v.14 ). It is estimated that between 20 and 60% of human genes have conserved miRNA target sites in the mouse [8, 9]. Several human and mouse studies have surveyed the tissue-specific expression of miRNAs using microarrays. Tissue-specific miRNAs, such as miR-122a in the liver, are expressed in both humans [10, 11] and mice [12, 13]. miRNA microarrays have also been combined with transcript expression arrays to search for correlations between miRNA and transcript expression, often in an effort to identify miRNA target genes. These studies have shown that miRNA expression and the expression of their target transcripts are often negatively correlated [14, 15], and that miRNAs physically located within or near mRNA transcripts are positively correlated [16, 17]. Microarray data can also be used to infer miRNA function by performing gene set analysis on transcripts that are correlated with miRNA expression, revealing functional coherence among the correlated transcripts.
Few investigators have studied the effect of genetic variation among individuals on miRNA expression. In humans, two groups have investigated miRNA variation in human cells. In lymphoblastoid cell lines, it was observed that miRNAs are both negatively and positively correlated with mRNA expression, and that miRNAs that are located near or within transcripts tend to be positively correlated . Another study used human cell lines to show that there are many miRNA that are positively correlated with each other, that miRNAs which are located within ~106 base pairs tend to be positively correlated with each other and that the transcripts correlated with miRNAs show functional coherence . In the mouse, quantitative trait locus mapping of miRNA expression showed that very few loci exist that regulate miRNA levels . In the present study, we profiled the expression of 577 miRNAs using microarrays in two panels of mice - a BXD recombinant inbred cross and a panel of commonly used laboratory inbred mice. We combined our data on miRNA levels with gene expression data from the same strains. We report that, in contrast to studies in multiple cancer tissues, liver miRNA expression is highly conserved across genotype. This observation, combined with the well known disease- and toxicant-induced changes in miRNA levels in liver, suggests that disruption in liver miRNAs expression may be a useful marker of liver injury, independent of genotype.
Archived frozen liver tissue samples collected from 70 strains of mice (males only, one mouse per strain) were used in these studies. The details of breeding, housing and tissue collection were described previously [21, 22]. Liver tissue collection (synchronized with respect to the time of day) from male mice of 36 C57BL/6JxDBA/2J (BXD) recombinant inbred strains, C57BL/6J and DBA/2J parentals, and B6D2F1 was conducted at the University of Tennessee Health Science Center in Memphis and approved by the local Institutional Animal Care and Use Committee . Additional liver tissues were collected from male mice of 34 inbred strains from the Mouse Diversity Panel (MDP) as detailed elsewhere  following a protocol approved by the Institutional Animal Care and Use Committee at the University of North Carolina at Chapel Hill.
Total RNA was isolated using the Ambion (Austin, TX) miRVANA kit according to the manufacturer’s instructions. RNA quality was assessed using the Agilent (Stanta Clara, CA) Bioanalyzer RNA Nano chip. Small RNAs (<40 nucleotides) were concentrated using the Ambion FlashPAGE system, applying 500 ug to each gel.
The Agilent G4472A mouse miRNA microarray was used to measure miRNA expression. Labeling and hybridization were carried out according to the manufacturer’s instructions using 10 ng of concentrated small RNA on each array. The Agilent miRNA microarray is a one color array and we performed quantile normalization on all arrays because a comparison of miRNA normalization techniques demonstrated the quantile normalization performed best by minimizing the mean square error in their experiment .
The BXD liver data set (GSE17522) and the MDP liver data set (GSE14563) were both run on the Agilent G4121A microarray platform. The data was downloaded from the University of North Carolina microarray database using the log2 Lowess normalized red/green ratio (mean background subtracted spot intensity). A spot was excluded if it had a Lowess normalized net signal < 10 in either channel. Transcripts detected on >90% of the arrays were retained and arrays with non-null values for >70% of transcripts were kept. A total of 16,738 (BXD) and 14,214 (MDP) transcripts were retained. The arrays were run in separate batches and batch effects were removed by subtracting the batch mean of each gene from the gene expression values. Final values were normalized by fitting a linear model with strain, sex, and the strain by sex interaction term to the gene expression data.
Transcripts with putative miRNA target binding sites were downloaded from the Targetscan , Miranda  and Pictar . For each miRNA target prediction list and for each of the 144 selected miRNAs, we calculated the correlation of the predicted target transcripts with each miRNA, as well as the correlation of the non-target transcripts with the same miRNA. The p-values were calculated by converting the Pearson correlations to t statistics and then calculating two-sided p values from the Student’s t distribution. We then performed a one-sided Student’s t test against the alternative hypothesis that the target-to-miRNA correlations are lower than the non-target-to-miRNA correlations. To determine an empirical p value for this test, we found the number of target transcripts per miRNA for each prediction algorithm, randomly selected 10,000 transcripts sets of the same size, and performed the t test. We then used the p value distribution from the random sampling to obtain an empirical p value that indicated whether the target transcripts were more negatively correlated with the miRNA than the non-target transcripts. Finally we calculate Bonferroni corrected empirical p values to account for multiple testing across 144 miRNAs and select miRNAs with significant negative correlation to their target transcripts (corrected p ≤ 0.05). This procedure was carried out separately in BXD and MDP.
For each of the 20,868 transcripts on the Agilent G4121A array, we downloaded the Miranda  predicted miRNA binding sites. We then intersected these sites with predicted SNPs, insertions, and deletions between C57BL/6J and DBA/2J from the mouse resequencing track in Ensembl to select 1,160 transcripts which may be influenced by differential miRNA binding. Of these, 943 were detected on the array in the BXD data set. For each transcript, we performed a two-sided Student’s t test of the transcript expression versus the BXD genotype in the 3′ UTR. Gene expression quantitative trait locus (eQTL) mapping data for the BXD strains was previously reported . A set of 2,379 significant local eQTL was selected using a q value  of 10% or less (p≤0.05). Distant or trans eQTL were not considered because there is no direct mechanism by which a SNP in a transcript miRNA binding site could cause a distant eQTL. Fisher’s exact test was conducted on a 2×2 contingency table with local/distant eQTL on one axis and SNP/no SNP in a miRNA binding site on the other axis.
BXD genotype data for 3,795 informative markers was downloaded from GeneNetwork and single marker mapping using FastMap  was performed on each of the 144 detected miRNAs. MDP genotype data, consisting of 156,000 SNPs, was used to perform haplotype association mapping in a 3–SNP sliding window . Multiple testing across SNPs was addressed by performing 1,000 permutations of the SNP data for each miRNA and calculating a permutation p value from the empirical distribution of maximum test statistics under permutation. Multiple testing across miRNAs was addressed by applying the q value algorithm  to the permutation p value from each miRNA and retaining those with qFDR<0.1. A local or cis eQTL was defined as a maximum QTL located within 5 Mb of the miRNA genomic location.
We profiled the expression of 577 miRNAs using the Agilent G4472 miRNA microarray in two panels of isogenic mice: a BXD panel consisting of 36 strains and an inbred laboratory mouse panel (Mouse Diversity Panel, MDP) consisting of 34 strains. The difference in expression among different miRNAs in liver covers a 10-fold range (Figure 1a). We calculated average intensity for each miRNA and selected the detected miRNAs as those with average intensity above the 75th percentile of the average intensity distribution of all miRNAs. The number of miRNA above various percentiles is shown in the legend of Figure 1a. Importantly, we found that liver miRNA expression was remarkably consistent among strains (Figure 1b), with the well-known liver-specific miR-122 showing the highest expression in all mouse strains.
Next, we compared the 29 miRNAs within the top 95th-percentile of liver miRNA expression in this study with that of other mouse and human studies (Figure 2). The most highly expressed miRNA in both the human and mouse liver is miR-122 . miR-21, miR-22, and miR-192 are also consistently reported as highly expressed in many studies. miR-720, the fourth most highly expressed miRNA in our study, was not reported by the other studies possibly because it was first reported in 2006 . The miR-30 family, of which miR-30c and miR-30a are the 9th and 15th most highly expressed in this study, has been reported as essential for hepatic development .
Next, we examined the correlation between miRNAs in the BXD and MDP data by hierarchically clustering the miRNA correlations in each panel (Figures 3a & b). Several clusters appear in both panels and are marked with colored bars between the two panels. We converted the Pearson correlation values to t statistics and then calculated p-values from the Student’s t distribution. Surprisingly, mRNA expression was not significantly correlated with the miRNAs in any cluster (p≤5×10−6, q≤0.1).
It has been reported that the expression of miRNAs that are closely located in the genome (approximately <0.5 Mb) tend to be positively correlated in human cell lines . We examined the correlation between miRNAs on the same chromosome and plotted this versus the genomic distance between the miRNAs (Figure 4). Indeed, miRNAs that are co-located have more highly correlated expression than more distantly linked miRNAs. We performed a parallel analysis of mRNA expression in each mouse panel and found no significant relationship between genomic distance and mRNA co-expression (data not shown).
In order to determine whether there are genomic loci that modulate expression of miRNAs in the liver, we performed miRNA QTL mapping in both the BXD and MDP data sets (Figure 5). At a threshold of q≤0.1 (p≤0.05), there were 33 and 18 significant QTL in the BXD and MDP data sets, respectively. Mir-345, located on Chr 12 at 109.2 Mb, has a local QTL at 108 Mb and is the only QTL that replicated between the two data sets. miR-31 and miR-34a have local QTLs in the BXD data set and this is consistent with a previous report of local QTLs for these two miRNAs . All the remaining QTLs are distant.
In order to understand how miRNAs regulate gene expression and vice versa, we correlated miRNA expression with gene expression in both panels of mouse strains. Calculating the Pearson correlation of all 144 expressed miRNAs with all expressed transcripts in both panels yielded 234 significant miRNA-to-transcript correlations in the BXD panel and 1,176 in the MDP at a q-value ≤0.1. We noted that the correlation p value distribution for some miRNAs was skewed toward 1.0 (Figure 6a). This diluted the significance of miRNAs which appeared to have significant correlation with transcripts (Figure 6b). We hypothesized that the large number of statistical tests along with the possibility that transcript degradation is not widespread in mammals diminished our ability to detect miRNA to target transcript correlations.
We examined the correlations in an miRNA-specific manner to see if there were certain miRNAs that were significantly positively or negatively correlated with transcripts. For each detected miRNA in each data set, we calculated the Pearson correlation of the miRNA with all transcripts and selected significantly correlated transcripts (q≤0.1). This is in contract to the analysis above in which all transcripts were correlated with all miRNAs together. There were 30 miRNAs with at least one significantly correlated transcript in the BXD data set and 53 in the MDP, with an overlap of 12 miRNAs. However, there was no overlap in the transcripts that were correlated with these 12 miRNAs between the two panels of mice. In the BXD panel, only miR-202-3p showed any GO or KEGG pathway enrichment for sensory perception of chemical stimulus (GO:0007606, p=3.275×10−10), a category with many G-coupled receptors. In the MDP, the transcripts correlated with miR-23a showed enrichment for ATP-dependent DNA helicase activity (GO:0004003, p=6.673×10−5) and the transcripts correlated with miR-27a showed enrichment for endonuclease activity (GO:0004519, p=1.107×10−4). Overall, the above approaches did not generate miRNA-transcript correlations that were reproducible in the two data sets or that provided a clear indication of the biological function of the transcripts that were correlated with each miRNA.
We also attempted a more targeted approach by downloading predicted miRNA gene targets from Targetscan , Miranda , and Pictar . For each target prediction algorithm and each miRNA, we calculated the correlation of all predicted targets and non-targets with the miRNA and performed a Student’s t test to see if the predicted target transcripts were more likely to be negatively correlated with miRNA expression than non-target transcripts. This procedure was performed in the BXD and MDP separately. MiR-101a had a significant negative correlation with its target transcripts in both mouse panels using all three prediction algorithms. MiR-101a has been associated with cyclooxygenase 2 (Cox2) expression and cell proliferation in mammary gland development , but we were unable to replicate this finding because Cox2 is not queried on the array platform that was used. The target transcripts of miR-101a showed an enrichment for the “positive regulation of gene expression” GO Biological Process (p = 4.108×10−10). We then looked for miRNAs that had consistent positive correlations with transcripts in both mouse panels. Among the Targetscan predictions, miR-142-3p was positively correlated with a set of 61 transcripts in both panels and miR-30e with 38 transcripts, but neither set exhibited any GO or KEGG enrichment.
It has been reported that miRNAs that are located within or near transcripts are co-expressed [16, 17]. We downloaded the genomic locations of all miRNAs from the Sanger miRNA data base and intersected these locations with the transcript locations of all transcripts on the array. Of the 144 detected miRNAs, 76 (52.8%) were located within transcripts. For each data set, we calculated the Pearson correlation of each miRNA with the transcript nearest to it and found no enrichment of significant positive correlations (p < 0.05) (Figure 7).
It is possible for SNPs in miRNAs to modulate target transcript levels by shifting binding affinity to target sites on transcript 3′ UTRs. We downloaded the sequence and location of all 600 mouse miRNA from the Sanger miRNA database and searched for SNPs in those regions using the mouse resequencing data on Ensembl. There were 21 miRNAs with SNPs that differed between C57BL/6J and DBA/2J haplotypes. Of these, miR-24–2 and miR-680 were detected in the BXD liver data set. miR-24–2 had a mean log2 expression of 9.4 in the BXD strains and contains a single SNP (rs49245174) in the pre-miRNA that does not affect the mature miRNA sequence. We investigated whether this SNP might nevertheless affect miRNA expression or processing by performing a Student’s t test of the miR-24–2 expression in strains containing the B allele versus the D allele of miR-24–2 and found no significant difference (p = 0.648). miR-680 has much lower expression (mean log2 5.95) in the BXD liver, but has 6 SNPs in two of its three genomic locations. miR-680–1 (Chr 6: 129.6415 Mb) contains 4 SNPs (rs51107210, rs45902352, rs1356740, rs50887828) in the pre-miRNA, none of which affects the mature sequence whereas miR-680–3 (Chr 12: 35.8795 Mb) contains 2 SNPs (rs49621977, rs48662351) in the pre-miRNA, one of which is within the mature sequence. However, despite these SNPs there were no detectible differences in miR-680 expression based on the allelic status at the miR-680 locus on Chr 12 (p = 0.549) and none of its target transcripts had a significant association to this locus (p ≤0.05).
miRNAs are thought to regulate translation by binding to partially complementary sites in 3′ UTRs. For each of the 20,868 transcripts on the array, we downloaded the Miranda  predicted miRNA binding sites. We then intersected these sites with predicted SNPs, insertions and deletions between C57BL/6J and DBA/2J from the mouse resequencing track in Ensembl, leaving 1,160 transcripts, of which 943 were in our BXD data set. Using the BXD genotype nearest the 3′ UTR of each gene and the expression of these transcripts, we performed a two-tailed Student’s t test against the null hypothesis that there was no difference in expression between the two parental alleles. We found 196 transcripts with B-vs-D allelic differences at a p values ≤0.01, consistent with the hypothesis that SNPs in miRNA binding sites could modulate expression. It should be noted that other local sequence variants could also contribute to these differences. We hypothesized that cis eQTL might be caused by SNPs in miRNA binding sites of these transcripts. We downloaded eQTL for the expressed transcripts in the BXD data set and selected 2,379 significant eQTL with q values ≤ 0.1 (p≤0.05). We then tested if the proportion of transcripts associated with local eQTL and with SNPs in potential miRNA binding sites was greater than the proportion of transcripts with distant eQTL (>5 MB from the transcript) with SNPs in miRNA binding sites using Fisher’s exact test. This test was highly significant (p=2.357×10−12) corroborating that genes with SNPs in their miRNA binding sites and more likely to have cis than trans eQTLs.
Most tissue-specific studies of miRNA expression have used a few samples of each tissue and did not vary the genotype; however, three studies have examined the effect of genetic variation on transcript and miRNA expression. Genome-wide transcriptional and miRNA profiling in human lymphoblastoid cell lines showed that more transcripts are positively correlated with miRNA expression than are negatively correlated . The transcripts that were correlated with miRNA expression also exhibited some functional enrichment. A similar analysis in 16 human cells lines found that some miRNAs were negatively correlated with their known target transcripts and that correlated transcripts showed functional coherence . Both studies also showed that miRNAs that are closely located in the genome (<0.5 Mb) tend to be positively correlated with each other. The only mouse study examined liver miRNA expression in wild type and ob/ob mutants derived from C57BL/6J and BTBR-T+tf/J strains and found approximately 150 consistently expressed miRNAs .
In this study, we evaluated the expression of 577 miRNAs in the livers of two panels of genetically diverse isogenic strains and observed that 144 miRNAs are expressed consistently in all samples. The first and most striking outcome of this study is the great degree of similarity in miRNA levels across all samples, reflected by the low variance of miRNAs in both mouse panels. This is consistent with the hypothesis that miRNA expression is important for liver homeostasis and the observation that Dicer1 is required for the overall mouse development . Paradoxically, liver development and function are largely preserved when miRNA processing is abolished through the selective knockout of Dicer1 in the liver [35, 36]. However, as the selective Dicer1 knockouts age, progressive liver damage is observed including increased apoptosis, cell proliferation, and inflammation , as well as development of hepatocellular carcinomas . These recent observations, combined with our data suggest that while Dicer and miRNAs have critical roles in hepatocyte survival, metabolism, developmental gene regulation, and tumor suppression in the liver, their effect on transcript abundance is subtle.
Several sets of miRNAs are co-expressed in both mouse panels and, consistent with previous observations , many of these are physically located near each other in the genome. Surprisingly, we were unable to detect significant positive correlation between miRNAs and the mRNAs within which they are in very close proximity, sometimes even within introns [15, 17]. Further, while several other approaches have successfully elucidated miRNA target genes [14, 37], we were unable to detect a significant number of miRNA targets using our approach. We attribute this to the low variance of miRNA expression which limited our ability to detect both positive correlations with host transcripts and negative correlations with putative target genes. These results suggest that miRNA target discovery is better conducted using tissue surveys , case/control studies  or targeted knock-down of specific miRNAs . The lack of GO category or KEGG pathway significance among transcripts that were correlated with miRNA expression may be due to the still poor understanding and annotation of genes in these databases .
miR-345 was the only miRNA with a QTL that replicated in the two panels of mice. While there is currently no known function for miR-345 in the liver, over-expression has been associated with progression toward oral squamous cell carcinoma  and miR-345 has been shown to target multidrug resistance-associated protein 1 (MRP1) in human MCF-7 cell lines . In the mouse, Abcc1, the ortholog of human MRP1, plays an important role in regulating cellular concentrations of endogenous and exogenous compounds. While there was no significant correlation between the expression of miR-345 and the Abcc1 transcript in our dataset, and while Abcc1 has no significant eQTL in either panel [21, 22], there may be an undetected association with protein levels. Mir-34a and miR-31 have both been shown to have local QTLs in the liver  and the replication of this finding in a separate inbred mouse cross between C57BL/6J and BTBR T+tf/J leads to increased confidence in these results. miR-34a was previously shown to be up-regulated by p53 and to be involved in the promotion of apoptosis [43, 44]. We searched for a SNP in the p53 binding site identified by Chang and colleagues  and found no polymorphisms between C57BL/6J and DBA/2J to indicate that the QTL is caused by genetic variation at this site. Thus, the replication of the QTL for miR-34a in two separate crosses indicates that more complex regulatory mechanisms may exist. miR-31 is part of a group of miRNAs that, while transcribed, are not always processed to the mature form and exported from the nucleus . There are no polymorphisms between C57BL/6J and DBA/2J haplotypes within 100 bases of miR-31. However, it is possible that differential processing of miR-31 leads to different levels of the mature miRNA in different strains.
This study shows that liver miRNA expression is highly conserved, in spite of significant genetic variation among the mouse strains used in this study. However, alterations in miRNA expression are associated with disruptions in liver homeostasis including non-alcoholic steatohepatitis , hepatocarcinogenesis , and as a result of a toxic insult [48–50]. Minor genetic perturbations have strong effects on mRNA expression [21, 22]; thus, even though toxic insults more easily perturb mRNA expression than miRNA expression, mRNA changes across a population are difficult to attribute to the effects of treatment alone and use as biomarkers in a genetically diverse population . In contrast, basal liver miRNA expression shows low variation in a genetically diverse population and alterations in miRNA expression in response to toxicants may prove more to be more informative as population-wide biomarkers.
In conclusion, we have surveyed the expression of liver miRNAs in 76 strains of inbred mice and found miRNA expression to be highly stable and conserved across strains. MiRNAs that are co-located in the genome tend to have expression profiles that are positively correlated, which is consistent with the hypothesis that these miRNAs are transcribed simultaneously from the same pre-miRNA transcript. Surprisingly, we found the expression of miRNAs that lie within the introns of genes did not have significant positive correlation with their host transcripts. We also found little negative correlation with the target transcripts of miRNAs and attribute this to the low variance of miRNA expression that occurs across genotypic variation. We identified several miRNA QTLs, some of which have been reported before, and note that the only one that replicates is a local QTL.
Gatti et al. highlights:
Financial support for these studies was provided, in part, by grants from the National Institutes of Health: R01 AA016285, P42 ES005948 and R01 ES015241.
Conflict of Interest Statement:
The authors declare that there are no conflicts of interest
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.