Computational Characterization of Exogenous MicroRNAs that Can Be Transferred into Human Circulation.

MicroRNAs have been long considered synthesized endogenously until very recent discoveries showing that human can absorb dietary microRNAs from animal and plant origins while the mechanism remains unknown. Compelling evidences of microRNAs from rice, milk, and honeysuckle transported to human blood and tissues have created a high volume of interests in the fundamental questions that which and how exogenous microRNAs can be transferred into human circulation and possibly exert functions in humans. Here we present an integrated genomics and computational analysis to study the potential deciding features of transportable microRNAs. Specifically, we analyzed all publicly available microRNAs, a total of 34,612 from 194 species, with 1,102 features derived from the microRNA sequence and structure. Through in-depth bioinformatics analysis, 8 groups of discriminative features have been used to characterize human circulating microRNAs and infer the likelihood that a microRNA will get transferred into human circulation. For example, 345 dietary microRNAs have been predicted as highly transportable candidates where 117 of them have identical sequences with their homologs in human and 73 are known to be associated with exosomes. Through a milk feeding experiment, we have validated 9 cow-milk microRNAs in human plasma using microRNA-sequencing analysis, including the top ranked microRNAs such as bta-miR-487b, miR-181b, and miR-421. The implications in health-related processes have been illustrated in the functional analysis. This work demonstrates the data-driven computational analysis is highly promising to study novel molecular characteristics of transportable microRNAs while bypassing the complex mechanistic details.

Introduction

Mature microRNAs (miRNAs) are a class of short non-coding RNAs, 21–25 nucleotides in length and endogenously transcribed in animals, plants, and viruses. These small molecules often regulate gene expression post-transcriptionally via base paring with complementary sites in target messenger RNAs (mRNAs) and either promote the degradation of mRNA or inhibit the translation of the mRNAs into proteins [1, 2]. In human, 2,588 known miRNAs (according to miRBase v21 [3]) have been estimated to target ~60% of human genes and regulate a vast array of fundamental cellular processes in different cell types [4].

Since miRNAs have been long considered to be synthesized endogenously, little has been studied on miRNA cross-species transportation during the past decade. It was very recently discovered that humans absorb a meaningful amount of certain exosomal miRNAs from cow’s milk, e.g., miR-29b and 200c; the endogenous miRNA synthesis does not compensate for dietary deficiency [5]; the biogenesis and function of such exogenous miRNAs are evidently health related [5-8]. While the evidence in support of milk-miRNA bioavailability is unambiguous, a recent report that mammals can absorb plant miRNAs (e.g. miR-168a) from rice [9], however, was met with widespread skepticism [10-13]. Based on these evidences, challenging questions may be raised regarding how human pick up miRNAs from dietary intake, why some exogenous miRNAs can be transferred into human circulation while others cannot, and what are the broader functional roles played by exogenous miRNAs in human disease processes.

A bioinformatics study is herein introduced to characterize the cross-species transportation of miRNA computationally where the following procedures have been employed. Firstly, through a comparative analysis across a large set of species, we systematically assessed the sequence conservation among all available miRNAs in the public databases. Current knowledge related to this issue is that miRNAs are well conserved in sharing common mature sequences, biosynthetic pathways and reaction mechanisms throughout evolution [14], while there is a large proportion newly evolved in each species and are considered to be species-specific [15]. Likewise, in this study, significantly different sequence profiles with some overlap are expected among species. Secondly, we applied a data mining strategy to identify discriminative molecular features that can classify miRNAs into different groups, e.g. different kingdom groups or circulating miRNAs versus the rest. Our initial list under evaluation covers the sequence features such as nucleotide composition, %G+C content and palindromic properties; the secondary structure of precursor miRNAs (pre-miRNAs); and the physicochemical properties, e.g., minimum free energy of the secondary structure. The rationale behind this collection is that functional study of miRNA has been largely depending on the target identification where sequences information is needed for identifying the complementary sites; and that miRNA gene recognition is mostly based on the prediction of pre-miRNA-like hairpin secondary structures that are conserved in closely related genomes. For example, current miRNA prediction methods have shown that sequential features, such as %G+C content and several normalized dinucleotide frequencies (%UA, %AA, %GC), are critical for detecting miRNAs from other types of non-coding RNAs [16-19]. In this study, all sequential and structural features that possibly capture the commonality and differentiation among miRNAs have been taken into account.

In addition, we know that extracellular miRNAs are found in circulation in two different forms: 1) associated with exosomes (also known as vesicles or microparticles) [20, 21], whose detailed molecular mechanism remains to be elucidated. Current studies show that microparticles exhibit highly distinct binding patterns with miRNAs, suggesting that there is a selection of miRNAs to be transported out of cells [22]. Hence the binding and transport mechanism may play a pivotal role in determining whether a miRNA will be excretory or not; 2) independent from exosomes/microvesicles, but instead bound to Argonaute (Ago) proteins as part of the RNAi silencing complex. Evidences suggest that the Ago-bound miRNAs may be the major form of miRNAs in blood circulation and their stability could be due to the binding with the Ago2 complex, which protects them from the RNAse degradation [23, 24], although the mechanism of miRNA-Ago2 complex secretion remains to be understood.

As there is a lack of prior knowledge of the secretory mechanism of miRNA to circulation, we plan to heavily rely on experimental data to identify features that can differentiate secreted miRNAs from the rest. Institutively, the secretory features should be highly associated with the intake and release mechanisms through transporting vesicles or the association with Ago proteins. In addition to the mature form of miRNA, we also include precursor sequences to possibly capture the editing associated features. Both structure- and sequence-based features are generated, including those related to the presence of branching and helical structure in pre-miRNAs and those describing the sequences with respect to their compositions of monomers and dimers, the existence of palindrome sequences, and the sequence length. While the precise effect of each feature on distinguishing secretory miRNAs from others is unclear, it is possible that these features could possibly contribute in recognizing whether the miRNAs are transportable by microvescicles, or measuring the strength of the miRNA-Argo2 complex formation. The binding strength between the miRNA and these proteins may inversely correspond to the likelihood of secretion. Based on the aforementioned features, we have conducted feature selection, followed by Manifold ranking analysis to infer the potential of exogenous miRNAs, particularly dietary miRNAs, being transported into human circulation. Experimental data was provided for validation.

Materials and Methods

A full description of the methods is provided in S1 Methods while a brief synopsis follows.

Data sets

The miRNA sequence and annotation data were downloaded from miRBase (Version 21) [3], which contains 34,612 mature miRNAs expressed from 28,421 stem-loop precursor sequence in 194 species. We first categorized the miRNAs into five kingdoms including Animalia, Plantae, Fungi, Protista and Viruses (detailed statistics is shown in Table 1). With the goal to find secretory miRNAs in human blood circulation, we adopted 360 human plasma miRNAs uncovered by Weber in 2010 [25].

Table 1

Detailed statistics of microRNA data, which includes a total of 34,612 mature sequences, 28,421 stem-loop precursor sequences, 194 species and 5 kingdoms.

Types	Animal	Plantae	Fungi	Viruses	Protista	Total
Mature miRNAs	26,705	7,645	84	152	26	34,612
Precursor miRNAs	21,257	6,990	53	91	30	28,421
Species	111	71	5	5	2	194

Various annotation information are collected from the following resources

1. miRBase [3], which complies the species/kingdom information of 34,612 mature miRNAs included in this study

2. DMD [26], which contains dietary species information of 5,217 miNRAs

3. Weber et al. [25] provides a list of 360 human circulating miRNAs

4. ExoCarta and EVpedia [27, 28] provide 370 exosomal miRNAs in human and mouse.

For assessment purpose, we have compiled a comprehensive collection of dietary miRNAs from literatures, a total of 5,217 miRNAs from 15 types of common food species such as cow’s milk, breast milk, tomato, grape, and apple fruit. All dietary miRNA information is accessible through our Dietary microRNA databases (DMD) [26]. In addition, annotation data also include exosome-associated information from ExoCarta and EVpedia [27, 28] for another dimension of assessment.

Feature collection

All features can be categorized into two classes: sequential features and secondary structural features. For each mature miRNA, a total of 1,102 features were generated including:

1,031 features calculated based on following sequences:

extend seed region sequence (first 8 nucleotides on 5’ end of mature miRNA sequence);

mature miRNA sequence;

corresponding precursor stem-loop sequence.

71 structural features identified based on the predicted secondary structure of precursor stem-loop sequence.

We note the key deciding factor of transportability might be related to the interaction between protein and miRNA. e.g. mature miRNAs may be associated with Ago proteins in cells [29], and the binding strength may inversely correspond to the likelihood of secretion. Hence, features that possibly associated with miRNA binding capabilities were examined, including the existence of palindromic sequences [30], sequence length and the compositions of monomers and dimers.

Secondary structural features were calculated based on the stem-loop structure of pre-miRNA. For example, RNAfold was employed to predict secondary structure and calculate Minimum Free Energy (MFE) [31]. Subsequently, 32 triplet features and 11 base-pairing features were calculated, such as A((( (frequency of 3 paired nucleotides leading by A) and %pairGC (length-normalized frequency of G-C pairing). NOBAI was utilized to compute Shannon Entropy (Q) and Frobenius Norm (F) [32]. The detailed descriptions and the references of each feature are given in Table A in S1 File.

Classification-based feature selection

Based on all aforementioned features, a support vector machine (SVM)-based feature elimination strategy was developed to identify features that can discriminate miRNAs of a certain class from others. The recursive feature elimination (RFE) based strategy has been employed to remove features irrelevant or negligible to the classification results in an iterative fashion [33-35]. Specifically, each iteration eliminates features with the lowest scores given by RFE. This process continues until a minimal subset of features is obtained while maintaining an acceptable level of classification performance.

We noted a major problem with our experimental dataset was its imbalance. For example, in the Plantae-against-Others case, the positive set that represents all Plantae miRNAs (7,645) was significantly outnumbered by the negative set (all miRNAs from other kingdoms, 26,967). To overcome the imbalance that presented challenges for SVM-based classification [36], synthetic minority over-sampling technique (SMOTE) [37] was utilized to produce a balanced dataset for each kingdom separation (Details in S1 Methods). We also grouped three minority kingdoms, namely, Fungi, Protista, and Viruses, into one virtual kingdom denoted as FPV.

Based on 5-fold cross validation, we evaluated the overall classification performance by calculating sensitivity, specificity, accuracy, and the Matthews correlation coefficient (MCC) [38]. It should be noted that, for each SVM-training and testing, we re-estimated the parameters by grid searching [39] and ensured optimized models were achieved for each classification. Last, the SVM-based feature elimination produced the minimal set of features that yields the best separation of one kingdom against others, and similarly, for the separation of circulating miRNAs against others miRNAs in human.

Manifold ranking to infer the miRNA transportability

Considering a large number of exogenous miRNAs might be transported into human circulation but have not been detected yet, which leads to a problem without well-defined negative sets, a different classification strategy, so-called ranking approach [40-42], can be alternatively employed. Here we built a model based on the identified discriminative features to rank miRNAs according to their potential of getting transported into circulation instead of predicting them to be transportable or not. The essence of such algorithms is as follows: the problem is defined on two datasets, a positive set, e.g. known secreted miRNAs, and a background set (an undetermined set which may include both positive and negative data); and the goal is to rank the individual members of the whole dataset according to their relevance to the positive data. A weighted graph is used to represent the whole dataset, with each data represented as a node, each pair of nodes as an edge and a weight defined as the similarity between the two nodes in the (to be identified) feature space. Then each positive data propagates its presence (as evidence) to its neighboring nodes to increase their relevance to the positive dataset, where this relevance is valued proportionally to the corresponding edge weight in the graph. An overall relevance score of each node is the sum over all the scores propagated to it from all the related positive data. One way to assess a ranking method is by checking the percentage of the positive training data that is ranked among the top X% of all the training data. Generally the higher the percentage is for each fixed X, the better the trained ranking algorithm is.

It has been well documented that Manifold Ranking algorithm (MR) helps in finding the most relevant samples from background to true positive datasets [43, 44]. In this study, we used all 360 human blood-detectable miRNAs as the positive set, and all other 34,252 miRNAs as background set in this experiment. The detailed description of MR can be found in the S1 Methods.

Functional inference through target analysis

The top-ranked miRNAs that are highly transportable were subject to further stratification according to their origins and if they are known exosomal miRNAs. As the functions of miRNA can be inferred based on its gene targets, we extracted the known human gene targets from CLASH dataset [45], miRTarBase [46] and DIANA-TarBase [47] if the dietary miRNA has identical sequences with human miRNA; otherwise, we predicted their targets in human using TargetScan [48] and miRDB [49]. Last, Gene ontology (GO) and pathway enrichment analysis [50] was carried out to infer the biological processes and functional pathways that the miRNA may get involved.

MiRNA-sequencing analysis on milk feeding study

A miRNA-sequencing analysis was conducted based on the archived human blood samples collected from a previous milk-feeding study [5]. These samples are from five health adult participants at four time points (0, 3, 6, 9 hours) after they consumed 1-liter bovine milk. In this study, both mRNA and microRNA were extracted from each blood samples at the BGI (Hong Kang, China) and the pooled miRNA was subject to small RNA sequencing analysis by using Illumnia-HiSeq2000. For bioinformatics analysis, the CAP-miSeq [51] was applied to identify both human and bovine microRNAs and calculate the expression. The miRBase (Version 21) [3] was used as reference library. We have carefully filtered out the low quality reads and strictly mapped the qualified reads to all known mature sequences, precursor sequences and the genomes of human and cow.

Data access

All the data and programs used in this analysis can be found at http://sbbi.unl.edu/publications/microrna.

Results

MiRNA sequence conservation across species

A total of 34,612 miRNA sequences from 194 species and five kingdoms are used for the initial comparative analysis. Although miRNA sequences have 21-25bps in length in general, skewed length distributions were shown with respect to the different kingdoms (Fig 1A). For example, compared to animal miRNA, the majority of viral miRNAs tend to have longer sequences.

Fig 1

Length distribution of mature miRNA sequences in 5 kingdoms (A) and schematic plot shows statistics of the cross-species sequence comparison (B).

(A) Length distribution of mature miRNA sequences in Animalia (red), Fungi (brown), Plantae (green), Protista (blue) and Viruses (purple) in both histogram (left) and Boxplots (right). (B) Schematic plot shows statistics of the cross-species sequence comparison. Within each species, light blue indicates the percentage of miRNAs that have homologues miRNA in Human, light purple represents the percentage of miRNAs that have homologues in other species within the same kingdom, and gray shows the percentage of miRNAs that have no homologues in any other species.

We doubt if the miRNA sequence conservation could be a feature contributing to the cross-species transportation. To test this, we compared all collected miRNA sequences across species using CD-HIT [52]. In total, 16,458 highly conserved clusters were derived (sequence identity higher than 0.98 with length variation no more than 1bp). We found most of species have miRNA homologs in other species within the same kingdom (Fig 1B, purple), e.g. 96 animal species share significant number of identical miRNA sequences with human (Fig 1B, blue). On the contrary, there are 18,154 (~52%) miRNAs that still lack of homologs in any other species (Fig 1B, gray), indicating each species gains specific miRNAs during evolution.

It seems to be quite rare that different kingdoms share identical mature sequences, which may partially explain why cross-kingdom transferring is challenging. For instance, among 7,645 plant miRNAs, none has identical or similar sequences in human, even using loose criteria allows up to 2 mismatches. In Fig 2, we illustrated the sequences conservation using a phylogenetic tree built on the precursor sequences of miR-190 and -171 families. It showed, among three miRNA gene clusters (miR-190a, miR-190b, miR-171), human miR-190a and -190b are close to many animal species, e.g. cow and mouse, within their respective clusters. However, a different gene cluster of plant miR-171 is closer to miR-190b, compared to miR-190a (Fig 2A). Specifically, human miRNA, hsa-miR-190b, show sequence identify of 79% and 77% with sly-miR-171a (tomato) and miR-190a (human), respectively (alignments shown in Fig 2B). It indicates while miRNA genes are often conserved among species or even across kingdom during evolution, the derived mature sequences, however, may vary from each other.

Fig 2

Phylogenetic tree of miR-190/171 family and sequence alignments of hsa-miR-190a/b and sly-MIR171a.

(A) Phylogenetic tree of miR-190/171 family based on the 162 precursor sequences from 89 species, where three major clusters are formed for miR-190a, miR-190b, and miR-171. (B) Alignments among precursor sequences of hsa-miR-190a/b and sly-MIR171a. Hsa-miR-190b show sequence identify of 79% (with blast similarity score = 28.3, E-value = 2E-05) and 77% (with blast similarity score = 76.1, E-value = 2E-18) with sly-miR-171a (tomato) and miR-190a (human), respectively.

A close look at the 2,588 human miRNAs shows that 930 of them share identical sequences with orthologs in other species. We suspect the exogenous miRNAs with identical sequences, if possibly getting into circulation, might be able to regulate the same gene targets in human; moreover, they might regulate the same homolog targets in their own species if other criteria are met, e.g. 3’ UTR of mRNAs are conserved across species.

MiRNA features related to cross-species transportation

Since sequence conservation alone cannot fully explain the miRNA cross-species bioavailability and molecular actions, we examined the aforementioned 1,102 features based on the sequence, structure and physicochemical properties to identify important features that can differentiate each kingdom group or distinguish human circulating miRNAs from the rest.

For each kingdom, we trained an SVM-based classifier wrapped by recursive feature elimination to select discriminative features associated with that kingdom. Based on 5-fold cross validation, we discovered a set of features that yields the best performance for each kingdom-against-others classification (Table 2). For example, in the Plants-against-other separation, we detected 147 features that produce a classifier with overall accuracy of 93.28% (Sensitivity = 89.71%, Specificity = 96.86%, MCC = 86.79%). Table 3 listed 21 features that contributed in two or more kingdom-wise classification. It is not surprising that the most top-ranked features were related with precursors, such as ensemble free energy, %pairGC and the %G+C content. Previous report shows that %G+C content may likely affects the stem-loop structure of pre-miRNA [53]. Moreover, several seed region features were included in this list, e.g. the frequency of “UUCC” in 5’ end strongly effected the Animalia- and FPV-against-others classification.

Table 2

Performance summary for kingdom-wise classification and human secreted miRNA prediction.

Classification	Selected Features	Accuracy	Sensitivity	Specificity	MCC
Animalia	166	93.29%	96.46%	90.11%	86.75%
Plantae	147	93.28%	89.71%	96.86%	86.79%
FPV a	126	89.79%	87.39%	92.19%	79.68%
Human secretion	96	90.03%	84.68%	95.37%	80.51%

The classification results are obtained through 5-fold cross validation with respect to different feature subsets selected.

a“FPV” denotes the virtual kingdom of Fungi, Protista and Viruses.

Table 3

Examples of overlapped discriminative features chosen by three kingdom-wise classifications and the human blood secretory prediction.

Features	Details	A	P	F	H	adj-P
Ensemble Free Energy (EFE)	Binding energy in kcal/mol	2	8	8		3.81E-02
Pairs	Number of pairing on stem-loop structure	1	2	19		9.43E-14
Minimum Free Energy (MFE)	Precursor sequence fold with least thermodynamic free energy in kcal/mol	7	9	7
%pairGC	Normalized frequency of G-C pairing on stem-loop structure.	10	24	5
Length_P	Length of precursor sequence	9	10	57
G(((	Frequency of 3 paired nucleotides leading with G	15	17	31
C(((	Frequency of 3 paired nucleotides leading with C	12	77	11		1.14E-03
freqUUCC_seed	Normalized frequency of UUCC in the seed region	29		33
freqCCA_seed	Normalized frequency of CCA in the seed region	27		35
freqCG_P	Normalized frequency of CG on precursor sequence	41		12
%G+C content_P	Normalized frequency of G and C on precursor sequence	16	12		57	5.14E-16
Max–stem-length	Longest stems block on stem-loop structure	39	20
A(((	Frequency of 3 paired nucleotides leading with A	26	1
freqGA_m	Normalized frequency of GA on mature miRNA	105	18
freqC_m	Normalized frequency of C on mature miRNA	126		94	30
freqCUG_P	Normalized frequency of CUG on precursor sequence	24	58	98	80
freqCUG_m	Normalized frequency of CUG on mature miRNA	122	120		92
freqCUGG_P	Normalized frequency of CUGG on precursor sequence	97	75		51
freqGU_P	Normalized frequency of GU on precursor sequence	137	101		52	5.78E-08
freqGU_m	Normalized frequency of GU on mature miRNA	51	103		5	2.27E-04
freqU_m	Normalized frequency of U on mature miRNA	114	53	127	48	6.89E-08
freqUGU_m	Normalized frequency of UGU on mature miRNA	101	98		93
Palindromes	Number of palindromes with length greater 3 occurring on precursor sequence	11	16	79	68	2.22E-16
freqC_seed	Frequency of C in the seed region	82		26	70
%G+C content_m	Frequency of CG that occurs on mature miRNA	93		70	37	3.90E-09

The numbers in the table represents the ranks of the selected feature in the corresponding classifiers (A: Animalia versus others, P: Plantae versus others, F: FPV versus others; H: human blood secretory miRNAs versus other human miRNAs) and the unselected ranks are not shown. The last column of “adj-P” shows that the adjusted p-value when analyzing the corresponding feature in the blood secretory prediction using Wilcoxon signed-rank test (insignificant p-values are not shown). Complete list was given in Table B in S1 File.

We also conducted the same feature selection on human circulating miRNA, where 96 features remained and the best performance for discriminating human blood miRNA from others can reach 90.03% accuracy (Table 2). We found most of these features are different from kingdom-wise features, except for 12 features such as number of palindromes of pre-miRNAs, %G+C content of mature miRNAs, and frequency of “C” in seed region (Table 3).

Taking into consideration all the features that are related to species and/or blood-secretion, we calculate a union of 221 features (categorized into 8 groups in Table 4) and believe the use of this hybrid feature set will render better prediction for transportable miRNAs in human circulation.

Table 4

221 features that have been selected for the final ranking of possible circulating miRNAs, which are categorized into eight groups according to the feature type and involved sequence type.

Feature groups	Number of selected features	Feature list
Frequency in seed region	28	AG, AGGU, C, CAGC, CAUC, CC, CCA, CCAG, CCAU, CCCA, CUUC, GA, GAG, GAGG, GCA, GCAG, GGU, GGUA, GU, GUA, GUAG, UA, UAG, UCC, UCCA, UGAG, UUC, UUCC
Frequency in mature miRNA	63	ACG, ACGG, AG, AGC, AGCU, AGG, AGGU, AGU, AGUA, AGUU, AUA, AUAG, AUCG, C, CAGU, CAUA, CC, CCA, CCAU, CCCA, CCG, CCGA, CG, CGA, CGAC, CGG, CGGA, CU, CUCC, CUG, GA, GAC, GACG, GAGC, GAGG, GCAC, GCUC, GGG, GGU, GGUA, GGUU, GU, GUA, GUAG, GUU, GUUG, U, UA, UAG, UAGG, UAGU, UC, UCC, UCCG, UCG, UG, UGA, UGU, UGUA, UGUG, UUCC, UUG, UUGU
Frequency in precursor sequence	80	ACCC, ACG, ACGA, ACGG, ACUA, AGCC, AGCG, AGG, AGGU, AGU, AGUA, AUA, AUAC, AUCG, C, CACG, CAG, CAGG, CAGU, CC, CCA, CCG, CCGA, CCGC, CCUG, CG, CGA, CGAC, CGAG, CGCC, CGG, CGGA, CGU, CGUG, CUA, CUAG, CUAU, CUCG, CUG, CUGG, CUGU, G, GA, GAAU, GACG, GCG, GCGA, GCUC, GCUG, GGCC, GGCG, GGU, GGUA, GGUU, GU, GUA, GUAG, GUUG, U, UA, UAC, UACG, UAG, UAGC, UAGG, UAGU, UAUA, UC, UCC, UCCG, UCG, UCGU, UCU, UG, UGCU, UGG, UGGC, UGGG, UUG, UUGU
Frequency of 3 nucleotides in stem loop structure	16	A(((, A((., A(.(, A.((, A.(., A. . ., C(((, C(.(, C.((, C. . ., G(((, G((., G(.(, G(., G. . ., U(((
Structure indicators	14	5 predicted shape type probability base on RNAshapes, MFE, Normalized_MFE, EFE, Normalized_EFE, freqMFEStructures, MFEI1, MFEI3, MFEI4, Shannon_Entropy
Stems/Pairs	12	%pairAU, %pairGC, %pairGU %A_unpaired, %C_unpaired, %G_unpaired max_stem_length, %G+C_stem, pairs, %pairs, %stems, Base-pairing-propensity
Percentage of nucleotides	4	%A+U_P, %A+U_m, %G+C content_P, %G+C content_m,
Length/Palindromes	4	Length_m, length_P, palindromes_P, palindromes_seed

Predicted transportable miRNAs

Since only 360 blood-detectable miRNAs (positive class) have been reported in previous study [25], we naturally assume that all other miRNAs may also possibly enter in human circulation. We performed a manifold ranking analysis on all 34,612 mature miRNAs based on the 221 selected features to rank miRNAs according to their transportable potential.

The final ranking list is given in Table C in S1 File. As expected, the query set of 360 known human plasma miRNAs were ranked among the top of the list. A close look at this list shows the top ranked entries are dominated by Animalia origin (Table 5). For example, 962 animal-borne miRNAs are ranked among top-1000 while 2812 are among the top-3000. Considering the percentages of miRNAs from Animalia, Plantae and Viruses in the original dataset are 77.16%, 22.09% and 0.44%, respectively, it indicates Animalia and Viruses miRNAs are highly enriched among the predictions of transportable miRNAs in blood circulation compared to others.

Table 5

Statistics of the top miRNA entries in the ranking list with respect to their origins.

	Animalia	Plantae	Viruses	Fungi	Protista	Dietary miRNAs
Original	26705 (77.16%)	7645 (22.09%)	152 (0.44%)	84 (0.24)	26 (0.08%)	5217 (15.07%)
Top-500	499 (99.8%)	1 (0.02%)	0	0	0	14 (2.8%)
Top-1000	962 (96.2%)	30 (3%)	8 (0.8%)	0	0	62 (6.2%)
Top-3000	2812 (93.7%)	163 (5.43%)	25 (0.87%)	0	0	273 (9.1%)
Top-5000	4678 (93.56%)	295 (5.9%)	27 (0.54%)	0	0	519 (10.38%)
Top-10000	9269 (92.69%)	670 (6.7%)	55 (0.55%)	4	2	1024 (10.24%)

There are 14 dietary miRNAs were ranked among top 500 and five of them have identical sequences in human including three bovine miRNAs (bta-miR-487b, -miR-421 and miR-216) and two chicken miRNAs (gga-miR-29a-3p and–miR-20b-5p). The identical sequence may indicate a higher chance that the exogenous miRNA will regulate human genes after transportation into circulation. As seen in Table 5, the number of dietary miRNAs scattered in the ranking list indicating the different likelihood of transportation. In particular, bta-miR-29b, a cow-milk miRNA, which we have previously validated in human blood circulation [5], is ranked as the 345th among all dietary miRNAs, which indicates there might be many other dietary targets to be explored in blood as a large screening is available. Among the top 345 dietary miRNAs including bta-miR-29b, there are 117 entries showing identical sequences with their homologs in human and 97 are exosome related. Intuitively, all exosomal miRNAs are highly likely to get into human blood circulation since exosomes are widely present in most of biological fluids.

In contrast, the brassica-specific miR-824 and miR-167a were ranked at the bottom of list, as the 31,502th and 29,669th, respectively, which is consistently with our previous discovery that they are the least detectable in circulation [5].

Validation of predicted transferrable miRNAs

From the prediction, the experimental data from cow milk study validated 9 transportable milk miRNAs in human blood, including bta-miR-487b, miR-181b, miR-421, miR-215, let-7c, miR-301a, miR-432, miR-127, and miR-184. The first three are highly-ranked in the dietary category and their functions are listed in Table 6.

Table 6

Gene targets and functional analysis of the three top predictions of the transportable miRNAs in cow’s milk, EBV, and rLCV.

	Dietary miRNAs	Human homologs with identical sequences	Number of targets		Related functional processes
			Experimentally validated	Predicated	Enriched pathway examples	Enriched GO term examples
Cow	bta-miR-487b	hsa-miR-487b-3p	46	468	Axon guidance (4.2E-2); Regulation of actin cytoskeleton (4.2E-2); Butanoate metabolism (5.2E-2); MAPK signaling pathway (3.4E-2).	Cell migration (8.0E-3); Positive regulation of locomotion (9.8E-3); Localization of cell (1.2E-2); Protein kinase activity (1.8E-2).
	bta-miR-181b	hsa-miR-181b-5p	1077	2084	Glioma (5.8E-3); Melanoma (8.1E-3); p53 Signaling Pathway (2.6E-2); Prostate cancer (1.0E-5).	Organelle lumen (7.3E-9); Nuclear envelope (9.5E-5); Intracellular non-membrane-bounded Organelle (1.4E-4).
	bta-miR-421	hsa-miR-421	793	808	Huntington's disease (3.5E-2); Hypoxia and p53 in the Cardiovascular system (4.6E-2); Apoptotic Signaling in Response to DNA Damage (3.6E-2).	Endomembrane system (4.6E-5); Nuclear lumen (7.8E-5); Establishment of protein localization (1.6E-4); Intracellular transport (2.8E-4).
EBV	ebv-mir-BART13-3p	-	-	208	Long-term potentiation (3.3E-2); Neurotrophin signaling pathway (3.9E-2); ErbB signaling pathway (6.1E-2); Melanogenesis (8.3E-2).	Phosphate metabolic process (1.2E-4); Wnt receptor signaling pathway (4.7E-2); hosphorylation (1.3E-3); Enzyme linked receptor protein signaling pathway (2.6E-3).
	ebv-mir-BART8-3p	-	-	549	Aldosterone-regulated sodium reabsorption (6.5E-3); Growth factors Survival factors Mitogens (2.3E-2); Prostate cancer (4.5E-2); Renal cell carcinoma (5.3E-2).	Nuclear lumen (5.6E-5); Transcriptional repressor complex (6.2E-5); Synapse (8.8E-5); Nucleoplasm (8.8E-5).
	ebv-mir-BART9-3p	-	-	568	Adherens junction (1.4E-2); Endocytosis (2.0E-2); Signaling Pathway from G-Protein Families (5.4E-2); Cell cycle (5.9E-2).	Membrane raft (1.5E-4); Regulation of transcription (1.8E-4); Regulation of transcription from RNA polymerase II promoter (3.4E-4); Positive regulation of macromolecule metabolic process (3.8E-4).
rLCV	rLCV-mir-rl1-16-3p	-	-	466	Fc gamma R-mediated phagocytosis (5.8E-4); Ubiquitin mediated proteolysis (6.0E-4); Endocytosis (6.6E-4); Melanoma (1.9E-3).	Regulation of endocytosis (1.0E-4); protein modification by small protein Conjugation or removal (1.4E-4); Extrinsic to membrane (5.7E-4).
	rLCV-mir-rl1-16-5p	-	-	403	Regulation of actin cytoskeleton (5.6E-4); Ca++/ Calmodulin-dependent Protein Kinase Activation (7.9E-3); Vascular smooth muscle contraction (9.3E-3); Focal adhesion (9.9E-3).	Positive regulation of cell migration (1.2E-4); Cytoskeleton (2.0E-4); Positive regulation of locomotion (2.6E-4).
	rLCV-mir-rl1-7-3p	-	-	228	Axon guidance (6.9E-3); Adherens junction (3.5E-2); Metabolism of Anandamide (5.7E-2).	Intrinsic to membrane (1.8E-4); Alkali metal ion binding (9.0E-4); Plasma membrane part (1.6E-3); Metal ion transport (1.7E-3).

The experimentally validated targets are collected from CLASH, MirTarBase and DIANA-TarBase; the complete list of the enriched pathways and GOs are listed in Table D in S1 File.

In addition, the top-ranked 9 Epstein–Barr virus (EBV) miRNAs (ebv-miR-BART9-5p, BART8-5p, BART9-3p, BART8-3p, BART14-3p, BART14-5p, BART15, BART13-5p, and BART13-3p) have been reported in [54]. These miRNAs show meaningful abundances in human B cells and they may cooperatively regulate several human genes in ebv-infected samples. Moreover, ebv-miR-BART13 and BART9 were proven to be involved in WNT signaling and cell cycle control in human [54], partially consistent with our analysis in Table 6.

Similarly, 14 miRNAs from Rhesus lymphocryptovirus (rLCV)(rLCV-miR-RL1-16-3p, RL1-16-5p, RL1-7-3p, RL1-7-5p, RL1-33-5p, RL1-33-3p, RL1-2-5p, RL1-24-3p, RL1-2-3p, RL1-24-5p, RL1-10-3p, RL1-10-5p, RL1-1-5p, and RL1-1-3p) that are highly transportable in our prediction have been reported in [55] where Raily et al. have found these rLCV miRNAs detectable in B cells of infected mammilla samples.

Based on all internal evaluation evidences, we provide a list of 368 exogenous miRNAs (23 viral miRNAs and 345 dietary miRNAs) as highly transportable miRNAs. The complete list can be found in Table D in S1 File. http://jvi.asm.org/content/84/10/5148.full.pdf

MiRNA-mediated gene regulations in human

For each miRNA that is potentially transferred into human circulation, 208 to 4,000 targets were collected through database search and computational prediction. The function and pathway enrichment analysis indicated that the 368 exogenous miRNAs may regulate human genes participating in immune development, metabolism and cancer. The detailed information for 9 exogenous miRNAs is provided in Table 6 while the full list is given in Table D in S1 File.

Theoretically, when human absorb meaningful amount of exogenous miRNAs from food, these confounders must successfully bind to human genes in order to make subsequent regulatory impacts on certain biological processes in human. To further assess this binding potential, we examined the sequence conservation among the targets in human and other species. Specifically, we collected the 3’UTR sequence of the target genes from different organisms and performed multiple sequence alignment based on the binding sites reported in TargetScan [48] and DIANA-TarBase [47]. For example, the top ranked cow-milk miRNA, bta-miR-487b, was confirmed in our validation and it shows identical sequence with hsa-miR-487b in human circulation. We compared the sequences of 15 predicted bovine target genes of bta-miR-487b and 46 experimentally validated targets of hsa-miR-487b in human. As shown in Fig 3, three conserved alignment blocks were observed among miRNA-mRNA binding regions in human and bovine. The consistency may provide more confidence if such exogenous miRNAs enter into human circulation, they may be able to play regulatory roles in human pathways by interacting with human genes. Based on our analysis, hsa-miR-487b targets 464 human genes targets and may be able to regulate human pathways related to MAPK signaling, actin cytoskeleton regulation, axon guidance, and Butanoate metabolism (Fig 4).

Fig 3

The multiple sequences alignment among the binding regions of 15 miR-487b targets in cow (_bta) and 46 targets in human (_hsa).

Fig 4

Regulatory network of bta-miR-487b in human.

Blue octagon nodes indicate genes that are involved in MAPK signaling pathway (adjusted fisher test p-value = 0.034); purple circle nodes indicate genes that are involved in regulation of actin cytoskeleton (p-value = 0.042); green triangle nodes represent genes that are involved in axon guidance (p-value = 0.042); pink square nodes denote genes that are involved in butanoate metabolism (p-value = 0.052). All light blue small circle nodes represent other predicted targets of miR-487b.

Another example is bta-miR-29b, which has also been experimentally validated in human blood [5]. Based on the 301 predicted mRNA targets, miR-29b is found to be involved in leukocyte transendothelial migration, cancer, and bone development. Overall, the transportable exogenous miRNAs predicted in this study are involved in many major biological processes including development, differentiation, cell proliferation, and metabolism [56], e.g. miR-27b, miR-34a, miR-106b, and miR-130 that are related to immune or development [6-8].Discussion

While our knowledge of miRNAs secretion and circulation is still limited, compelling evidences has indicated there is an selective intake and release mechanism involved in these processes. Our study has followed this line to explore the mechanistic features that may contribute in miRNA cross-species transfer and gene regulation in human using an integrative approach. Through sequence comparison, miRNAs from different species show moderate conservations among mature sequences throughout phylogeny. Subsequently, various sets of features related to sequence, structure and physicochemical properties are found to be discriminative for miRNAs in different kingdom groups and blood secretory group. The selected feature contributing to blood secretion may reflect molecular mechanism related to selective package and exportation [57], carrier-mediated transport realized by its encapsulation in exosomes and microvesicles or Ago2-bound complexes, and the microparticles exhibit highly distinct binding patterns with miRNAs [22] in which, intuitively, involved certain molecular sequence, structure, or physicochemical properties.

Selected features may bring new insights of transposable miRNAs. For example, the length of pre-miRNAs and %G+C content of mature miRNAs show different patterns between human circulating miRNAs and the rest of human miRNAs (shown in S1 Fig), suggesting human blood miRNAs are produced by longer pre-miRNAs and often show higher percentage of C, G nucleotides. In the kingdom-wise classifications, several selected features were related to the frequency of nucleotide G in the first segment of miRNAs, i.e., the 6–7 nucleotides of 5’ end of miRNAs. This could result from the following. For target recognition by two groups of miRNAs, each recognizes its mRNA targets by 5’ or 3’ end complementary pairing. The first 6 or 7 nucleotides on the 5’ end are known to be used for target recognition with little or no support from the 3’ miRNA end [58]. This suggests that 5’ end and its nucleotide composition are important factor in determining the fate of miRNAs. A recent study showed that strand bias selection exists for miRNAs in incorporation into the RISC complex; and highly expressed strands tend to have nucleotide G-bias and U-bias at 5’ end [59]. All these clues suggest that miRNAs enriched with G and U nucleotides at 5’ end are more likely to bind to the Ago2 protein, forming a RISC complex.

Within the top-1000 ranked prediction, 96.1% miRNAs are from animal origin and only 3% are from plant, which is consistent with our intuition that animal-borne miRNAs are subject to more significant absorption in human compared to plant miRNAs. However, it should be noted the bioavailability of milk miRNAs has not been investigated at a large scale, and the uptake mechanism is still ambiguous regarding which and how miRNAs enter blood circulation. In contrast, it was shown that rice miR-168a (osa-miR-168a) is also detectable in human and animal sera, and it decreases the expression of low-density lipoprotein receptor adapter protein 1 (LDLRAP1) mRNA [60]. Nonetheless, the low concentration reported by multiple follow-up studies seems to exclude any impact of these miRNAs on gene expression. For example, the levels of osa-miR-168a in human plasma were only about 3% of the bta-miR-29b levels observed in our preliminary studies. It is possible that the miRNAs from plant have sequential or structural features that prevent their secretion into blood, or that the methylation of the 3’-terminal ribose in position C2 in plant miRNAs by the methyltransferase HEN1 [61], impairs the intestinal transport of miRNAs, but this hypothesis is currently untested. We also expect the interaction between exsome and host intestinal cells may influence the transport. An in-depth investigation of transport mechanisms and kinetics of milk-borne miRNAs was beyond the scope of this study, but is currently pursued in the investigator’s lab.

Another critical challenge for uncovering the diverse biological roles of miRNAs lies in the efficient identification of targeting genes where current computational methods are still at a very early stage of focusing on static miRNA target prediction [62], while new observations have revealed the dynamic nature of miRNA-mRNA interactions that may vary in different phenotypic conditions [63-66]. Our on-going efforts are focused on the integration of gene expression information into target prediction toward identifying the real regulatory events under a pathway context. Empowered by the next-generation sequencing technology, we can study miRNA existence and expression in different specifies. However, sequencing based analysis on cross-species transportation study still encounter challenges in terms of the sensitivity of detecting exogenous miRNAs with low abundance and differentiation of the sources when identical sequences are involved. With that has been said, such computational study is important to provide an efficient tool that can facilitate a targeted search for exogenous miRNAs in human circulation rather than profiling in the old fashion.

Conclusion

Here we presented an integrative study where comparative analysis and computational prediction have been applied to assess the cross-species transportation of miRNAs, particularly focusing on inferring the likelihood of exogenous miRNA in human circulation. Given the limited understanding about miRNA circulation, this study will contribute substantially in overcoming the aforementioned scientific limitations and dramatically reducing the extensive lab-load in miRNA biology research by using a revolutionary systems-driven strategy to study this complex problem. Specifically, this bioinformatics-driven study enables bypass the following key issues: (1) Lack of supporting information to discern between endogenous miRNA synthesis or dietary miRNA absorption in the miRNA expression change in human blood test subjects; (2) Inference from endogenous miRNA synthesis [67] that might compensate for dietary miRNA deficiency; (3) potential distinct metabolism of dietary miRNAs in the intestinal mucosa. Substantial follow-up studies will be conducted to extend the analysis and clarify in greater detail the information generated by this study in revealing information on miRNA exchange and functional regulation in human disease prevention. We anticipate the novel computational tools developed for characterizing miRNA circulation and targeting will be useful for other miRNA and nutrigenomics research areas.

Distributions of three example features.

(A) Number of palindromic sequences in precursor sequence in each kingdom. The x-axis contains five kingdoms and the y-axis represents the number of palindromic sequences. (B) Length of pre-miRNA (left) and %G+C content on mature miRNAs (right) over the miRNAs found in blood and the other human miRNAs. The x-axis contains the value of the corresponding feature, and the y-axis represents the frequency.

(PPTX)

Click here for additional data file.

Supplementary Tables.

Table A. Detailed descriptions of all 1102 miRNA features; Table B. All overlapped discriminative features chosen by three kingdom-wise classifications and the human blood secretory prediction; Table C. Final manifold ranking list of all 34,612 miRNAs; Table D. Gene targets and functional analysis of 368 predicted transferrable exogenous miRNAs.

(XLS)

Click here for additional data file.

Supplementary Methods.

(DOCX)

Click here for additional data file.