In a sample of 69 unrelated individuals, we identified a total of 773 polymorphisms (700 SNPs and 73 insertion-deletions) in pre-miRNA hairpins which passed strict quality control filters (Figure 1B). Of these 773 variants, 411 mutations occurred in pre-miRNA stem/loop regions, 242 in mature miRNA, and 120 in the seed sequences (Figure 1B and Additional file 1: Table S1). Among these polymorphisms, we identified 198 previously undescribed mutations that are currently not present in dbSNP v135. The number of alleles per base was calculated separately for each region of the miRNA (stem-loop, mature, and seed). In our dataset, allele frequencies were slightly higher in the stem-loop region compared to the mature miRNA and seed region (0.013, 0.011 and 0.011, respectively).

To control for somatic mutations, undescribed mutations found in a single individual (ie. singletons) were removed from analyses and 33 novel variants present in multiple individuals (non-singletons) were analyzed (Figure 1B). Among the non-singleton mutations identified in pre-miRNA hairpins, 5 novel mutations were in highly conserved miRNA seed sequences (Figure 1B and Additional file 1: Table S2). The first novel seed variant, a “C/T” SNP at chr3 128081086, was located in the 3’ strand of human miRNA hsa-mir-1280 and was present in one Hadza and one Sandawe, two hunter-gatherer populations from Tanzania. The second, a “C/T” SNP at chr1 62544469, was located in the 3’ strand of hsa-mir-942 and was present in one Hadza and two Sandawe individuals. Two novel seed sequence indels found at low frequency were an “ACA” deletion in miRNA hsa-mir-4483 at chr10 115537763-115537766, found in 2 Yoruban individuals, and a “T”-allele insertion in hsa-mir-3940 at chr19 6416444 in 2 individuals with northern European ancestry from Utah. Lastly, a “CT” deletion located on chromosome 6 in the 3’ strand of hsa-mir-4640 was found in 9 individuals from 7 global populations (namely, 1 Pygmy, 2 Sandawe, 1 Yoruba, 1 Maasai, 2 African-Americans, 1 Mexican-American and 1 Gujarati Indian individual). Using miRNA target prediction software, we determined that in the absence of the “CT” mutation in miRNA hsa-mir-4640, the 3’ strand was predicted to individually target 316 binding sites in 79 genes (Additional file 1: Table S3). With the “CT” deletion, however, the number of predicted targets dropped to 11 binding sites in 3 genes, where gene targets did not overlap between the original and modified miRNA.

We compared levels of nucleotide diversity in miRNA to other genomic sequence classes (for examples, exons and introns) in each population group using Watterson’s estimator of theta (θ_W) (Figure 2A). We found in our high-coverage sequence data that miRNAs were among the most conserved sequences in the genome, at the same level as exons (Figure 2A), consistent with results from prior lower-coverage sequencing studies [41]. Overall, however, African populations had the highest level of nucleotide diversity across all sequence classes (θ_W = 5.0 ± 0.7 × 10⁴) compared to European and Asian populations (θ_W = 3.6 ± 0.2 × 10⁴) (Figure 2A).

In addition, we identified 319 population-specific miRNA alleles (PSMAs), defined as variants present exclusively in 1 of the 14 globally sampled populations. About two-thirds (66.8%) of PSMAs were present in African populations and the proportion of population-specific alleles (ie. population-specific density) was highest in Africans relative to non-Africans, consistent with prior analyses of human genetic variation [41]. In particular, Pygmy hunter-gatherers had the highest population-specific density with 11 PSMAs per genome (Figure 2B). Among the Hadza hunter-gatherers, we found two high frequency PSMAs with allele frequency ≥ 50% in the Hadza only, a novel “A/C” SNP in the stem-loop of hsa-mir-1291 and an “A/G” SNP (rs111566161) in the 3’ mature sequence of miRNA hsa-mir-4711. This “A/G” SNP was found to be exclusively shared by two Hadza and a Southern Kalahari San individual, a population thought to share an ancient genetic ancestry with the Hadza and other African click speaking populations [43,44].

To measure population differentiation, we calculated pairwise F_ST at variants in miRNA genes (See Methods). F_ST ranges from 0 to 1 with an estimate of 0 indicating no population differentiation and an estimate of 1 indicating complete differentiation. Estimates of F_ST at miRNA variants were calculated hierarchically: (a) between individual populations (Figure 3A-B), (b) between major geographic regions (for example, Europe-Africa, Europe-Asia, Africa-Asia) (Additional file 2: Figure S1) and (c) between pooled African and pooled non-African populations (Figure 3C). We then compared the individual estimates of F_ST to empirical distributions of miRNA F_ST. The pairwise F_ST values that were outliers (ie. within the top 5% of the empirical distribution or above the 95^th percentile) were classified as population-differentiated (PD) and variants with these extreme values were inferred to be enriched for SNPs under recent selection [45]. Our data showed that among major geographic regions, African populations had the highest average F_ST values among populations, with 34 PD-miRNA alleles between Africa and Asia, 33 PD-miRNA alleles between Europe and Africa, and 18 PD-miRNA alleles between Europe and Asia (Figure 3B; Additional file 2: Figure S1).

In a comparison of pooled African and pooled non-African populations, we identified 31 PD-miRNAs containing variants with outlier F_ST values above the 95^th percentile of the distribution of miRNA F_ST values (F_ST ≥ 0.186 with p < 0.05) (Figure 3C and S2; Additional file 1: Table S4), with 4 PD-miRNA variants in seed sequences (Additional file 2: Figure S3). Furthermore, when we apply a more stringent criteria for population-differentiation among miRNA variants, we found 8 highly population-differentiated miRNAs (HPD-miRNAs) with variants above the 99^th percentile representing the top 1% of all pairwise F_ST estimates (F_ST ≥ 0.366; p < 0.05) (Figure 3C and Additional file 1: Table S4).

Experimentally-validated mRNA targets of the 8 HPD-miRNAs were identified by querying 45 publically available deep sequencing and microarray datasets (See Methods). We found that these 8 HPD-miRNAs experimentally down-regulated the expression of 2,139 unique human mRNA target genes in at least two datasets (Figure 4A). Target gene enrichment analysis revealed that 72 of the 2,139 mRNA targets were significantly over-regulated (p < 0.05) by these 8 HPD-miRNAs compared to a null or random set of 8 miRNA (See Methods) (Figure 4A-B and Additional file 1: Table S5). The genes of these over-regulated mRNA targets were involved in immune response, metabolism, developmental processes, cell communication, transport and response to stress (Additional file 2: Figure S4). Additionally, 14 of the 72 enriched gene targets were reported to be candidate loci in genome-wide association studies (GWAS) (Additional file 1: Table S6).

Functional enrichment was performed by two methods for the 2,139 mRNA gene targets of HPD-miRNA. First, genomic functional enrichment was used to determine if particular biological functions regulated by the 8 HPD-miRNAs were overrepresented compared to the null expectation across the genome (Figure 4B) [46]. Based on this analysis, we identified 319 statistically over-represented functions (p < 0.05; n = 319) (Figure 4). The top hits included sugar (mannose, fructose, and glucose) metabolism and the regulation of insulin (p < 10⁻⁷) (Additional file 1: Table S7). Secondly, we performed a bootstrapping analysis of miRNA functional enrichment to account for the non-random distribution of miRNA targets throughout the genome and identified biological pathways overrepresented by HPD-miRNAs as compared to other non population-differentiated miRNA sets (Figure 4B). In particular, the 2,139 experimentally-validated gene targets of the 8 HPD-miRNAs were found to function in 5,475 annotated biological processes based on the ENSEMBL gene ontology. Each of the 5,475 functions were assigned a probability of being regulated by a random null set of 8 miRNAs (See Methods) and this probability was compared to the observed regulation by the 8 HPD-miRNAs. We found 55 biological processes that were significantly enriched for regulation by our set of 8 HPD-miRNAs (p < 0.05) (Additional file 1: Table S8), and identified 120 HPD-miRNA gene targets involved in enriched biological processes; of particular interest, these targets include DICER1 and BRCA1 among other genes (Additional file 1: Table S8).

In addition, clustering analysis of the 55 significant biological processes, based on gene ontology similarity, revealed nine major pathways affected specifically by the 8 HPD-miRNAs: the mitosis pathway, axons and morphogenesis, response to toxins, signaling, acid regulation, transcription, ion sequestration, muscle movement and immune response (Additional file 2: Figure S5). The mitosis pathway contained the highest number of significantly over-represented processes, specifically mitotic transitions (G1 and metaphase/anaphase transition), cell cycle checkpoints, and cytokinesis, representing over 20% of the 55 biological processes. Interestingly, the altered expression of mitotic genes, which can occur through miRNA regulation, is one of the known characteristics of the onset and progression of cancer.

Disease associations were also identified for all 31 PD-miRNAs through the miRNA disease database MiRGator [47] (Figure 5A). Of the 31 PD-miRNAs, 7 miRNAs (hsa-mir-202, hsa-mir-196a-2, hsa-mir-423, hsa-mir-943, hsa-mir-520h, hsa-mir-1908 and hsa-mir-647) were identified in prior studies as being differentially expressed in human diseases representing a significant increase above expectation (1.1 out of 31) for disease-associated miRNAs (95% C.I. = 1.08 – 1.11; p = 2.2 × 10⁻¹⁶) (Figure 5B). The 7 PD-miRNAs were associated with 37 diseases, most extensively cancer, with all 7 PD-miRNAs being implicated in cancer risk or as biomarkers for one or multiple cancers, specifically, breast cancer, prostate cancer, testicular cancer, ovarian cancer, lung cancer, renal cell carcinoma, gastric cancer, pancreatic cancer, head and neck cancer, colon cancer, esophageal cancer, brain tumors, thyroid cancer, glioblastoma, glioma, and germ cell tumors (Figure 5C). The hsa-mir-196a-2 T-allele at SNP rs11614913 has been significantly associated with increased risk for esophageal cancer in non-smoking European males [48] but decreased risk for breast, lung and gastric cancers in Chinese populations [48]. In the present study, we observed a significantly lower frequency of the hsa-mir-196a-2 T-allele at SNP rs11614913 in Africans compared to non-African populations (F_ST = 0.41; p < 0.001) (Figure 5D).

Notably, one of our PD-miRNAs, hsa-mir-202, contains a T-allele with known effect on miRNA expression and a protective effect on breast cancer mortality [42]. In our dataset, we found that the frequency of the T-allele at hsa-mir-202 was lower in African populations (26%) compared to non-African populations (65%), on average (Figure 6A and Additional file 1: Table S9). We observed similar frequencies of the T-allele in African and non-African populations in the 1000 Genomes Project data, that sampled only one African population (YRI) (16%) and 2 non-African populations (CEU, CHB + JPT) (83.3% and 92%, respectively) (See Methods). In addition, we also observed variability in the frequency of the T-allele within Africa; specifically, the T-allele occurred at lower frequency in the non-hunter-gatherer Yoruba, Maasai, and Luhya populations ((0.08%) compared to hunter-gather Pygmy, Hadza, and Sandawe populations (43%), on average (Figure 6A and Additional file 1: Table S9). We also found the T-allele at moderate frequency in African-American populations (Figure 6B).

Based on analysis of high-coverage whole genome sequence data in 69 individuals from 14 worldwide populations, we identified 33 novel polymorphisms in miRNA genes and 31 variants with high levels of genetic population differentiation between Africans and non-Africans (PD-variants). Five novel variants (2 SNPs and 3 indels) were located in miRNA seed sequences where they could have a large effect on the number, strength, and specificity of each miRNA/mRNA target interaction. Specifically, a novel “CT” deletion found in 7 of the 14 globally diverse populations altered the predicted mRNA targets of hsa-mir-4640 to include 3 additional targets and removed all of its 79 original predicted targets from regulation, indicative of the disruptive potential of seed sequencing indels on miRNA function. In addition, among the 31 PD-variants we found 8 HPD-miRNAs with variants that lie in the top 1% of pairwise estimates of F_ST, indicating extensive population differentiation at these loci between African and non-African populations consistent with a model of local adaptation (Figure 4) [5,45,49].

From functional enrichment analysis we observe that HPD-miRNAs are significantly enriched for regulation of genes involved in the glucose/insulin metabolism pathway (p < 10⁻⁷), and cellular division (p < 0.001), notably the regulation of genes involved in mitosis (specifically, mitotic checkpoints, transitions, and cytokinesis) (Figures 4 and 5; Additional file 1: Tables S7 and S8). Aberrant cell division during mitosis or aberrant gene expression levels during chromosome segregation often result in chromosomal instability, which is a key diagnostic feature of most cancers [50]. The disruption of cell division by altering gene expression levels in response to mitotic instability has been strongly correlated with tumor development and progression; both in vitro and in vivo evidence have demonstrated that in the absence of other cell cycle and DNA repair defects, mitotic disruption can transform cells and predispose them toward cancer [50]. Given that mis-regulation of cellular division is the hallmark of cancer, it is striking that miRNAs with highly population-differentiated alleles are observed to be significantly enriched for regulation of mitotic pathway genes including genes such as S100A8 and P2RX3 whose expression profiles are currently used as biomarkers for multiple cancers.

In addition, we identified 7 population-differentiated miRNAs where expression level differences of these PD-miRNAs have been correlated with cancer and other disease phenotypes (Figure 5B-C). Among the identified cancers, higher mortality rates have been reported for breast, ovarian, gastric, prostate and testicular-germ cell cancers in individuals of recent African ancestry compared to individuals of either European or Asian descent (p < 0.001) [51-54]. Of particular interest is hsa-mir-202 which contained one of the most highly population-differentiated variants in our dataset and is one of two miRNAs currently under investigation as a circulating blood-based marker for the detection of non-Hodgkin lymphoma and early stage breast cancer [42,55]. Recent in vitro functional data demonstrated that the T-allele was protective against breast cancer mortality by first increasing mature hsa-mir-202 expression levels, leading to subsequent down-regulation of its gene targets, including cancer related genes CRYBB2, DICER1, SART1, S100A8, P2RX3, and BRCA1[42]. Diminished expression of mature hsa-mir-202 in individuals harboring at least one non T-allele resulted in a significantly elevated risk of non-Hodgkin lymphoma (OR = 1.83, 95% CI: 1.17–2.85; P = 0.008) [42]. Our data showed that African and African American populations had a lower frequency of the T-allele compared to European and Asian populations, suggesting decreased baseline expression levels of mature hsa-mir-202 in African populations. In the context of tissue specific expression, hsa-mir-202 expression varies considerably by tissue type and is more highly expressed in prostate cells compared to breast or ovarian tissues (Additional file 2: Figure S7). Therefore, lower frequencies of the T-allele in African populations have the potential to reduce gene expression levels down below critical thresholds in the breast and ovarian tissues. This reduction may contribute to population differences seen in the onset of breast and ovarian cancer in women of African and European descent (Additional file 2: Figure S7).

Similarly, the other 6 PD-miRNAs identified (hsa-mir-423, hsa-mir-196a-2, hsa-mir-520h, hsa-mir-1908, hsa-mir-647 and hsa-mir-943) have also been implicated in cancer susceptibility in human populations [56]. First, the hsa-mir-423 SNP rs6505162 has been shown to confer reduced risk for breast cancer in women of European decent in GWAS [57]. Second, the hsa-mir-196a-2 SNP rs11614913 CC genotype has been significantly associated with increased risk for breast, lung and gastric cancers in Chinese populations; conversely, the homozygous TT genotype has been significantly associated with esophageal cancer in non-smoking European males [48]. Based on these results, recent studies have called for more detailed analysis of the frequency of this allele in different ethnic groups [58]. We observed a significantly higher frequency of the hsa-mir-196a-2 C-allele at SNP rs11614913 in African compared to non-African populations (F_ST = 0.41; p < 0.001) (Figure 5D). Third, hsa-mir-520h expression was determined to be significantly associated with E1A-mediated tumor suppression and cell migration during cancer metastasis and inhibition of hsa-mir-520h significantly decreased the downstream ability of cancer cells to migrate and invade other areas of the body [59]. This pattern has also been observed consistently in different types of cancer including pancreatic, breast and ovarian cancer [59-61]. Also, up-regulation of hsa-mir-520h was shown to increase the effects of the anticancer drug resveratrol in slowing lung cancer tumor mobility [62]. Finally, multiple studies have linked miRNAs hsa-mir-1908, hsa-mir-647 and hsa-mir-943 expression to various cancers known to have ethnic specific disparities [63-65]. Overall, these studies demonstrate that genetic variability within miRNA has the potential to vary miRNA expression and/or mRNA target binding which can be strongly correlated with the onset of multiple cancers, the progression of cancer metastasis and the response to drug therapies. We demonstrate that the frequency of clinically important genomic miRNA variants varies significantly among ethnic populations, particularly between African and non-African groups. Thus, we suggest that population-differentiated variation in miRNA may contribute to ethnic disparities seen in certain forms of cancer.

High quality whole genome sequencing (~60× coverage) was obtained for 69 globally diverse individuals from publically-available datasets. Fifteen African hunter-gathers were obtained from Lachance et al. 2012 [39], including 5 Pygmy (Py) (three Baka, one Bakola, and one Bedzan), 5 Hadza (Hz) (plus two technical replicates), and 5 Sandawe (Sw) using the Complete Genomics sequencing platform [68,69].

Additionally, 54 unrelated individuals were obtained directly from Complete Genomics including 9 individuals of Northern European ancestry (combined as CEPH - including 4 CEPH and 5 CEU individuals), 9 individuals of Yoruban ancestry (YRI), 5 individuals of Mexican ancestry (MEX), 5 African-Americans living in Dallas (ASW), and 2 individuals of Puerto Rican ancestry (PUR), and 4 individuals each of Toscani Italians (TSI), Japanese ancestry (JPT), Han Chinese (CHB), Gujarati Indian (GIH), Maasai Kenyan ancestry (MKK), Luhya Kenyan ancestry (LWK) [70]. Variants were defined as autosomal alleles that differ from the human reference genome build (GRCh37/hg19), and novel variants are defined as variants that are absent from dbSNP (db135). Complete Genomic Data is available through their public site: http://www.completegenomics.com/public-data/69-Genomes/.

Pre-miRNA locations and mature miRNA sequences were downloaded from the database mirbase (v18) and included 1524 known pre-miRNAs. We filtered known pre-miRNAs into high and low confidence miRNA, based on experimental validation of their function as target suppressors [71]. High confidence miRNAs were classified as those with at least one experimentally-validated mRNA target, using 45 public datasets (See Messenger RNA Target Enrichment and Analysis). Among novel variants, 165 mutations were identified within high confidence miRNA, 83 were found in miRNA stem-loops, 82 in mature miRNA, and 28 in seed sequences. Seed sequence locations were computed in R software and defined as the 2^nd – 8^th base pair of the mature miRNA [72]. Novel miRNA target prediction was done in DIANA software [73] for hsa-mir-4640 “CT” deletion. Predicted targets were also filtered for those with experimental validation, based on 45 public datasets (See Messenger RNA Target Enrichment and Analysis). Nucleotide diversity for each sequence class was measured using Watterson’s estimator of θ (θ_W) [74].

Whole genome sequence data were annotated using the UCSC genome browser for sequence class determination. Pairwise population F_ST values were calculated using Weir and Cockerham’s weighted equations adjusting for small sample size using the R statistical software [75,76]. MiRNAs were determined to be population differentiated (PD-miRNA) if they contained an allele with a F_ST value in the top 95^th percentile of the empirical distribution of F_ST values (F_ST ≥ top 5% with p < 0.05). Highly population differentiated miRNA (HPD-miRNA) were defined as miRNA containing variants with F_ST values above the 99^th percentile between African and non-Africans (F_ST ≥ top 1% with p < 0.05). P-values for F_ST’s were determined by testing the allele frequency difference at each allele between African and non-African samples using a Welsh two-sided t-test (Additional file 2: Figure S2).

For HPD-miRNA, messenger RNA (mRNA) targets were identified from the consensus of 45 paired miRNA/mRNA experimental datasets, and subject to the following quality control filters: (1) the mRNA targets have to be correlated with miRNA expression with a correlation coefficient of r < −0.5 and (2) that this level of correlation was observed in at least 2 publicly available experimental deep sequencing or microarray datasets. Experimental datasets consisted of: 6 ENCODE (Encyclopedia of DNA elements) cell line RNA sequencing datasets (GM1287- a lymphoblastoid cell line produced from the blood of a female donor with northern and western European ancestry, H1_hESC- a human embryonic stem cell line, Hela_S3- an immortalized cell line from an African-American female patient with cervical cancer, K562- an immortalized cell line from a female patient with chronic myelogenous leukemia (CML), HepG2- a cell line produced from a male patient with liver carcinoma, and NHEK- a epidermal keratinocyte cell line), 13 cancer deep sequencing datasets from the Cancer Genome Atlas (TCGA) (BLCA - Bladder Urothelial Carcinoma, BRCA - Breast invasive carcinoma, COAD - Colon adenocarcinoma, HNSC - Head and Neck squamous cell carcinoma, KIRC - Kidney renal clear cell carcinoma, KIRP - Kidney renal papillary cell carcinoma, LAML - Acute Myeloid Leukemia, LIHC - Liver hepatocellular carcinoma, LUAD - Lung adenocarcinoma, LUSC - Lung squamous cell carcinoma, READ - Rectum adenocarcinoma, STAD - Stomach adenocarcinoma, and UCEC - Uterine Corpus Endometrioid Carcinoma), 2 deep sequencing Gene Expression Omnibus (GEO) datasets (GSE31999 and GSE37765), and 24 microarray datasets from GEO (GSE2564, GSE9234, GSE11255, GSE12250, GSE14224, GSE14473, GSE14794, GSE14834, GSE15387, GSE15745, GSE16558, GSE16654, GSE16759, GSE17306, GSE17491, GSE17498, GSE18155, GSE18693, GSE18899, GSE19350, GSE20692, GSE21032_1, GSE21032_2, GSE21321). Dataset correlation was done through MiRGator v3 software [77].

Two HPD-miRNA (hsa-mir-5007 and hsa-mir-4634) had no targets that passed quality control filters and were excluded from downstream analyses. The number of targets that passed quality control filters for the remaining HPD-miRNA were hsa-mir-202 (1145 targets), hsa-mir-1304 (710 targets), hsa-mir-1269a (488 targets), hsa-mir-4482-1 (144 targets), hsa-mir-449c (69 targets), and hsa-mir-4707 (1 target). Each unique mRNA target was individually analyzed using bootstrapping analysis to test for significant enrichment in HPD-miRNA. Specifically, targets were identified for (n = 2000) randomly chosen sets of 8 miRNA. The distribution of random targets was then compared to the 2,139 targets regulated by the 8 HPD-miRNA, and 72 targets were significantly enriched for regulation by HPD-miRNA (p < 0.05) (Figure 4 and Additional file 1: Table S5). Enriched targets were then annotated with the PANTHER pathway software (Additional file 2: Figure S4).

MiRNA target prediction was done with and without the novel “CT” deletion observed in hsa-mir-4640 using DIANA prediction software [73]. Canonical hsa-mir-4640 targets were compared among 6 different target prediction algorithms (TargetScan, miRNAorg, Microcosm Targets, PITA, PICTAR and miRDB) using MiRGator v3 software [77].

Experimentally validated gene targets were analyzed for reported genome wide significance in a GWAS study using the NHGRI GWAS catalog (Additional file 1: Table S6) [78]. We determine the extent of linkage disequilibrium (LD) in the regions surrounding the 8 HPD-miRNA SNPs, by using common tagging SNPs from the Illumina 1 M-Duo array in African samples (Additional file 2: Figure S6). The program PHASE v.2.1, which implements a Bayesian statistical method [79], was used to reconstruct multi-site haplotypes from genotype data for 15 SNPs from the Illumina 1 M-Duo SNP array flanking the miRNAs of interest on chromosomes 4, 5, 10,11,13 and 14 in 697 African individuals (Additional file 2: Figure S6). Haploview [80] was then used to calculate pairwise measures of LD among SNP loci, creating a graphical representation of the LD relationships among these loci (Additional file 2: Figure S6).

Genomic functional enrichment was analyzed with GOEAST software to determine statistically overrepresented GO terms within our set of 2,139 mRNA gene targets [46]. GO terms include annotated biological processes, molecular functions, and cellular components. The GOEAST algorithm assumes genes should be evenly distributed across the genome. MiRNA, as a class, may favor targeting a specific range of biological processes. To address this fact we perform miRNA functional enrichment.

MiRNA functional enrichment identifies biological pathways overrepresented by HPD-miRNA as compared to a set of 8 randomly sampled miRNA. We found the 2,139 mRNA gene targets of HPD-miRNA to be involved in 5,475 annotated biological processes using the ENSEMBL Gene Ontology database (Figure 4) [81]. For miRNA functional enrichment, the null distribution was created by, (1) randomly resampling 8 miRNA from the mirbase database, (2) identifying all mRNA targets (meeting the quality control filters described above in Messenger RNA Target Enrichment and Analysis), (3) identifying all biological processes associated with mRNA targets (using the ENSEMBL Gene Ontology database). Steps 1–3 are replicated n = 2000 times. The frequency of obtaining each of the 5,475 biological processes seen in HPD-miRNA is calculated. Finally, the actual frequency of each of the 5,475 biological processes associated with HPD-miRNAs is then compared to the null distribution above to ascertain enrichment of biological processes in our sample. Biological processes overrepresented in HPD-miRNA, having p-values < 0.05, were considered significantly enriched in HPD-miRNA and further clustered by gene ontology similarities using Revigo software [82].

Disease associations were assessed through the MiRGator v2 disease database [47]. Disease associations are defined as differentially expressed miRNA in a published gene expression profile [47]. The significance of disease association results was tested by bootstrapping 31 miRNAs from mirbase and testing for disease association, n = 10,000 times. Network analysis was done using R with disease network visualization performed with Gephi software [83].