Analysis of deep sequencing microRNA expression profile from human embryonic stem cells derived mesenchymal stem cells reveals possible role of let-7 microRNA family in downstream targeting of hepatic nuclear factor 4 alpha

Abstract

Background: Recent literature has revealed that genetic exchange of microRNA between cells can be essential for cell-cell communication, tissue-specificity and developmental processes. In stem cells, as in other cells, this can be accomplished through microvesicles or exosome mediated transfer. However, molecular profiles and functions of microRNAs within the cells and in their exosomes are poorly studied. Next generation sequencing technologies could provide a broad-spectrum of microRNAs and their expression and identify possible microRNA targets. In this work, we performed deep sequencing of microRNAs to understand the profile and expression of the microRNAs in microvesicles and intracellular environment of human embryonic stem cells derived mesenchymal stem cells (hES-MSC). We outline a workflow pertaining to visualizing, statistical analysis and interpreting deep sequencing data of known intracellular and extracellular microRNAs from hES-MSC). We utilized these results of which directed our attention towards establishing hepatic nuclear factor 4 alpha (HNF4A) as a downstream target of let-7 family of microRNAs.

Results: In our study, significant differences in expression profile of microRNAs were found in the intracellular and extracellular environment of hES-MSC. However, a high level of let-7 family of microRNAs is predominant in both intra- and extra- cellular samples of hES-MSC. Further results derived from visualization of our alignment data and network analysis showed that let-7 family microRNAs could affect the downstream target HNF4A, which is a known endodermal differentiation marker. The elevated presence of let-7 microRNA in both intracellular and extra cellular environment further suggests a possible intercellular signalling mechanism through microvesicles transfer. We suggest that let-7 family microRNAs might play a signalling role via such a mechanism amongst populations of stem cells in maintaining self renewal property by suppressing HNF4A expression. This is in line with recent paradigm where microRNAs regulate self-renewal and differentiation pathways of embryonic stem cells by forming an integral biological network with transcription factors.

Conclusion: In summary, our study using a combination of alignment, statistical and network analysis tools to examine deep sequencing data of microRNAs in hES-MSC has led to a result that (i) identifies intracellular and exosome microRNA expression profiles of hES-MSC with a possible mechanism of miRNA mediated intercellular regulation by these cells and (ii) placed HNF4A within the cross roads of regulation by the let-7 family of microRNAs.

Figure 1

Different stages of processing microRNA sequences to gain biological insight. Massive sequencing data can be overwhelming to examine, with each stage of processing choices are made to direct analysis towards biological significance and understanding. Initial processing stages emphases on visualization of datasets in different forms such as annotations and alignment in the genome. This is to allow for quick assessment of salient features. Following this are more advanced analysis exploiting existing tools for examining distributions and networks topology to conjure biological hypothesis for verification. The entire workflow can then be thought as a funneling process towards biologically interesting interactions for verification in our hES-MSC cell line.

Figure 2

High abundant expression signals of detected miRNAs in two independent replicates show similar frequency distributions. Distributions of the expression signal count (y-axis: representing the proportion of total transcript counted and the x-axis: representing the count of signals; data in a log scale) derived from biological replicates of intracellular samples show a high degree of similarity. The empirical distributions can be interpreted as a mixture of two essentially distinct distributions. The low-abundant count (presented on the left part of the distribution) provides noise-rich data. This data shows an exponential decline of the frequency distribution suggesting that a large portion of signal have low counts. This is followed by a fairly even distribution across the log scale of counts and ended with sporadic signal that have high counts at the far end of the log scale. The high-abundant count (presented on the right part of the distribution) has a long tail and shows log-normal-like frequency distribution. The last part of the distribution shows high level of similarity in the right-size trends, which is is exploited in statistical tests for determining significances/reliable signals from mostly noise signals.

Figure 3

Venn Diagrams depict the difference of miRNA transcript counts identified in intracellular and extracellular (secreted) miRNA sets. Due to statistical results presented on Figure 3, the transcript count domains are broadly separated into 3 main levels: Low (below 32 transcripts), Mid-range (32-10000) and High (>10000). Within each sub-interval of the miRNA expression value, Venn diagrams are used to analysis the degree of overlapping between experimental replicates. Kappa correlation coefficient, reflecting the degree of agreement between two sets (Cytel Studio@7) was calculated for each pair of replicates. In general, for intracellular compartment miRNAs the degree of overlapping (agreement between experiments) for the Low category is much lower compared to the mid-range and high data subset. This result suggests a good experimental responsibility of occurrence of miRNAs expressed at the moderate- and high-abundant expression levels. The low-abundant set of miRNA shows poor reproducibility, perhaps due to experimental and biological noise. For the extracellular miRNA samples, there were no transcripts at the high expression category. Moreover, a low probability of co-occurrence of the same miRNAs (non-significant kappa correlation) in low expression category or in moderate expression category is found.

Figure 4

Additional Venn diagram analysis. A) Considerable overlap in top ranked genes from comparing of highly-expressed microarray and sequencing data for intracellular miRNAs. B) Overlap in top ranked genes from comparing of highly-expressed microarray and sequencing data for extracellular miRNAs is low. C) Top-level sequence data shows strong agreement between intra and extra cellular samples. Expression microarray for miRNAs were performed and compared with our deep sequencing data. Ranking of the top miRNAs by their abundance from both techniques was performed. Top hundred miRNAs for the intracellular samples were compared using kappa correlation analysis and a large degree of overlap between studied sets was observed. Similarly for the extra cellular samples, the top 50 miRNAs were compared. The degree of overlap between studied sets was less than that of the intracellular samples. A final comparison involves that of the common genes amongst replicates of the intracellular samples with the extra cellular samples. In this case we found that most of the genes of the extracellular microRNA form a common subset of the genes of the intracellular miRNA.

Figure 5

Different hES-MSC replicates show similar alignment peaks. Visualization of Seqmap mapping via UCSC genome browser reveals peaks along specific genomic regions that have large numbers of miRNA binding to these regions. In the top figure A Blue and Purple trend lines represent biological replicates of samples derived from intracellular environment. High degree of correlation is observed between the replicates as can be observed from the similarity in locations where peaks were found. The height of each such peak corresponds to the number of transcripts detected from deep sequencing. Each peak now represents genomic locations where a large number of specific transcripts bind to. The bottom figure B depicts the extra cellular sample transcripts that are aligned to the human genome using Seqmap. Peaks occur in similar region after mapping across the replicates. Each replicate is visualized with a different color and the salient feature reveals peaks from different sample aligning in similar locations.

Figure 6

Histograms of peak magnitudes show a similar distribution after an initial uneven trend. Distributions of peak magnitudes [Y-axis shows the frequency of a particular transcript associated with the peak occurring; x-axis is a measure of the transcript counts in log scale] after mapping shows the same trend after an initial uneven distribution. Green and red trend lines are least square fits of gamma distribution function reflecting that both samples showed great similarity after the initial uneven distribution possibly suggesting that the initial phase to be noise.

Figure 7

32 transcripts is an optimum threshold for achieving similar cumulative frequency distribution amongst replicates. KS statistics was applied iteratively to our biological replicates (red and black lines), each graph depicts a gradual change in the threshold value. The KS test statistics can be thought of as a cost function that we seek to minimize to ensure that the distributions between the two replicates are similar. Applying such a strategy iteratively as we change the threshold value gradually, we arrive at a point where both line converge together indicating similar distributions. The point where this first occurs is designated as the minimum threshold count value of biological significance.

Figure 8

HNF4A is a common hub for networks derived from alignment data and TargetScan predictions. Gene interaction network on the left (A) is derived from the dataset of genes with overlapping regions corresponding to peaks from previous mapping. The other gene interaction networks (B) is derived from computationally predicted gene targets from TargetScan. Comparing both gene interaction network, similar topology was observed with HNF4A as a node amongst the interactions suggesting HNF4A as a possible downstream target for let-7 family miRNAs.

Figure 9

Expression of genes in the HNF4A alignment network follows a similar pattern as HNF4A in hES-MSC and HEPG2 cells. This figure shows the relative gene expression of genes from gene interaction network (A) in Figure 8 in HEPG2 cells compared to hES-MSC. It is interesting to note that most genes in this network are suppressed when let-7 family miRNAs are over-expressed (hES-MSC) and up-regulated when let-7 family gene expression goes down (HEPG2 cells).