Genomic features and computational identification of human microRNAs under long-range developmental regulation

Defining self-transcribed miRNAs and the extent of their pri-mRNAs

Based on the assumption that some miRNAs have co-evolved with their cis-acting regulatory elements, in this study, we focused on conserved human self-transcribed miRNAs (ST miRNAs), which are assumed to be transcribed from their own promoters and have an ortholog in at least one other vertebrate species (Additional file 1, Table S1). These requirements are consistent with the work describing the original definition of a GRB [17] as well as the subsequent analysis of the features of GRB target genes [25]. In general, ST miRNAs include all conserved human intergenic miRNAs except those transcribed by RNA polymerase III (Pol III) [30], as well as known cases of intronic miRNAs whose fate after whole-genome duplication in fish indicates that they have their own promoters, independent of their host genes (see Methods).

Since proximal miRNAs can appear in genomic clusters that share the same promoter and are transcribed as a single transcript, we clustered ST miRNAs according to their genomic vicinity (see Methods) and analyzed only one promoter per cluster,

The annotation of the pri-miRNAs we used in this study was the result of a combinatorial approach by Saini et al. [14] that employed predicted and experimental evidence to identify both transcription start and end sites. One of the most important reasons for choosing this dataset was the fact that they considered a larger range of possible miRNA sizes, thereby avoiding biases introduced by excluding very long transcripts. We estimated the general scale of pri-miRNAs based on these data, where the TSS is within 50 kb upstream of pre-miRNAs and the transcription end is within 20 kb downstream of the pre-miRNAs. This approach was sensitive enough to include around 84% of the pri-miRNAs in the initial dataset (see Figure 1 and Methods).

Enhancer density

HCNEs

HCNEs represent putative enhancers of gene expression, whose frequency rises strongly around the conserved developmental genes they usually target [17, 22]. We therefore first investigated whether this applied to conserved ST miRNAs as well, by comparing the HCNE density in a 300 kb window centered on the ST miRNA with random coding and non-coding genomic regions with similar gene density (see Methods). In a few cases (~8%), the location of the ST miRNA search window overlapped with a previously defined GRB making it difficult to identify the actual target gene, which could either be the miRNA, the annotated GRB target gene, or both. For this reason, we compared the HCNE densities both including and excluding cases in which the ST miRNA search window overlapped with a previously defined GRB. As shown in Table 1, the HCNE density was significantly higher around ST miRNAs than the random coding and non-coding regions in both comparisons, a trend that is independent of the lineages compared (p-value ≤ 0.05, two-sided bootstrapped version of the Kolmogorov-Smirnov test: see Table 1, Figure 2 as well as Additional files 2 and 3). Thus, in support of our original hypothesis, we can conclude that HCNEs are highly over-represented in the genomic vicinity of conserved human ST miRNAs. The HCNE density was also compared between the conserved ST miRNAs and the GRB target genes annotated in Akalin et al. [25]. We found that the HCNE density was significantly lower for the ST miRNAs (p-value ≤ 0.05, two-sided bootstrapped version of the Kolmogorov-Smirnov test: see Table 1, Figure 2 as well as Additional files 2 and 3). This is due to the fact that the GRB target genes were annotated based on their high HCNE density [25], which was not a prerequisite for constructing our dataset of ST miRNAs.

	Comparison between human conserved ST miRNAs and human random coding regions		Comparison between human conserved ST miRNAs and human random non-coding regions		Comparison between human conserved ST miRNAs and human GRB target genes
Lineage comparison	p -value (all regions)	p -value (excluding regions overlapping GRBs)	p -value (all regions)	p -value (excluding regions overlapping GRBs)	p -value (all regions)	p -value (excluding conserved human ST miRNAs overlapping GRBs)
human: mouse	0	0	0	0	0	0
human: dog	0	0	0	0	0	0
human: opossum	0	0	1.0e-04	0	0	0
human: platypus	0	0	0	2.0e-04	0	0
human: chicken	0	0	0	0	0	0
human: frog	0	0	0	0	0	0
human: zebrafish	4.0e-03	8.0e-03	1.8e-03	4.0e-04	0	0

The comparisons represented by each column were performed after selecting HCNEs and ST miRNAs conserved between the lineages shown in the left column (see Methods). The p-values were computed using the two-sided bootstrapped version of the Kolmogorov-Smirnov test. All p-values lower than 1.0e-20 were set to 0.

Experimentally verified enhancers

In addition to HCNEs, we utilized experimental data describing the location of the transcriptional co-activator p300 to analyze the enhancer density around ST miRNAs. This transcription factor is a ubiquitous component of enhancer-associated protein assemblies. It co-localizes with active enhancers and plays a critical role during embryonic development [31–36]. Visel et al. [36] mapped the genome-wide binding of p300 in forebrain, midbrain and limb tissue of developing mouse embryos using ChIP-seq technology. This generated ~5000 p300 binding regions which are associated with active enhancers. Since conservation criteria did not play a role in determining these p300 binding regions, they represent a dataset of experimentally verified enhancers that is completely independent of the HCNEs, allowing us to test our hypothesis in an unbiased manner. We then compared the p300 binding density around ST miRNAs that are conserved between mouse and human (mouse:human) with that of mouse coding and non-coding regions using the same method as the previous HCNE density analysis and found that the density of p300 binding sites around mouse:human ST miRNAs was indeed significantly higher than in the control sets (p- value < 1e-20 for both comparisons, two-sided bootstrapped version of the Kolmogorov-Smirnov test) (Figure 3).

These results are independent of HCNE density, since the same trend was still present when we performed the analysis excluding p300 binding sites that overlap with mouse orthologs of our HCNEs (p-value < 1e-20 for both comparisons, two-sided bootstrapped version of the Kolmogorov-Smirnov test) (Figure 3, Additional file 4 and Table S2).

Comparison between enhancer enriched and enhancer poor ST miRNAs

As the density of HCNEs/p300 binding around conserved ST miRNAs is significantly higher than in random non-coding and protein coding regions and there are conserved ST miRNAs that are known developmental regulators similar to GRB target genes, we wanted to see whether there was an association between the individual features of GRB target genes and HCNE/p300 densities. We again limited the analysis to ST miRNAs conserved between human and mouse in order to be able to incorporate the p300 binding data.

We then compared the features of GRB target genes between the HCNE/p300 enriched and HCNE/p300 poor miRNAs excluding intragenic ST miRNAs due to the uncertainty in associating genomic features with the ST miRNA or the host gene.

GRB target gene features can identify ST miRNAs subject to long-range regulation

The miR-9 family of miRNAs is a known, experimentally verified GRB target [17] and therefore a prime example for illustrating how the genomic features we analyzed could serve to annotate miRNAs under long-range regulation. This family is specifically expressed in brain and affects the fate of ES cell-derived neural precursor cells differentiating along the glial or neuronal pathways [41]. The expression of the miR-9 family of genes is dynamically regulated [42] during differentiation and development, and the human miR-9 family has three members: hsa-mir-9-1, -2, and -3. In the human genome, hsa-mir-9-2 and hsa-mir-9-3 are intergenic miRNAs, and hsa-mir-9-1 is found in the second intron of C1orf61. Since the zebrafish ortholog of hsa-mir-9-1 is located in an intergenic region, it qualifies as a ST miRNA based on our definition. Moreover, hsa-mir-9-1 and hsa-mir-9-2 are paralogs that were most likely separated by the whole-genome duplication at the root of jawed vertebrates [43–45], which suggests that their common ancestor was intronic, but that the host gene is not required for their transcription [17]. Therefore, all miR-9 family members can be classified as ST miRNAs and are likely to have their own promoters. A detailed examination of their genomic environments showed that all of them share the features of GRB target genes (Figure 4), since they map to genomic regions with high HCNE densities, the putative promoters of their mouse orthologs map to bivalent domains and they are all associated with several proximal CGIs.

We further examined the annotation of each gene within the investigated region around the miRNAs. Since there is no other putative GRB target gene within the region of hsa-mir-9-3, we conclude that hsa-mir-9-3 is most likely the only target of long-range enhancers in that region. In the neighborhood of hsa-mir-9-1 and hsa-mir-9-2, we find myocyte enhancer factors MEF2D and MEF2C, respectively. Both of these genes are transcription factors that play a role in myogenesis and are therefore plausible GRB target gene candidates. However, it was shown that two regions with enhancer activity located ~10 kb downstream of dre-mir-9-1, the zebrafish ortholog of hsa-mir-9-1, and ~100 kb downstream of dre-mir-9-5, the zebrafish ortholog of hsa-mir-9-2, gave the reporter gene an expression pattern similar to that of zebrafish miR-9, but not the zebrafish myocyte enhancer factors [17]. Thus, it may very well be that both miRNAs and the myocyte enhancer factors are regulated by distinct elements in the regions, and that there is a boundary between the two, similar to that of neighboring GRBs targeting PAX6 and WT1 genes [46]. In summary, the genomic features of regions around members of the miR-9 family display characteristics equivalent to those of protein-coding GRB target genes, lending further support to the use of these features for predicting novel miRNA targets of long-range regulation.

Genome-wide identification of putative GRB target miRNAs

Taken together, the results of our analyses show that HCNE enriched ST miRNA genes are more likely to be associated with the features of GRB target genes than HCNE poor ST miRNA genes, which leads us to believe that we can utilize them for the prediction of GRB target miRNAs.

Assemblies and annotation

The following genome assemblies were used: human NCBI 36, mouse NCBI 37 and mouse NCBI 36, chicken v2.1, Xenopus tropicalis v4.1 (US DoE JOINT genome Institute), zebrafish Zv7 (The Wellcome Trust Sanger Institute). Genomic coordinates and sequences of known miRNAs were obtained from miRBase (Version 12)[1] and the genomic coordinates of HCNEs were made available from the Ancora database ([29]; http://ancora.genereg.net; Additional file 6, Table S4). CGI annotations, Broad Institute Chromatin State Mapping using ChiP-Seq and genomic sequences were extracted from the UCSC Genome Browser Database ([47]; http://genome.ucsc.edu) and gene annotations for human, mouse and zebrafish were obtained from Ensembl ([48]; http://www.ensembl.org; version 52) using Biomart ([49]; http://www.biomart.org). In addition, because of the limited annotation of the zebrafish genome, we also included zebrafish orthologs of human genes predicted based on the similarity of exonic sequences and gene structures between zebrafish and human genomes [50]. The list of putative GRB target genes was constructed by inspecting HCNE density peaks, tracing the teleost gene/HCNE synteny, and analyzing functional gene annotation [25]. GRBs are human-zebrafish synteny blocks overlapping with GRB target genes. Synteny blocks were defined through the joined zebrafish to human high scoring net alignments downloaded from the UCSC Genome Browser [47].

miRNA classification

Human and zebrafish miRNAs were classified into intergenic and intragenic miRNAs based on the overlap of the pre-miRNA with known protein coding genes (at least 1 bp overlapping). Our ST miRNAs included both human intergenic miRNAs and human intragenic miRNAs that were intergenic in zebrafish genome, since these are likely to be self-transcribed. We removed the Pol III-transcribed miRNAs mentioned by Borchert et al [30] and subsequently retrieved orthologs of human ST miRNAs conserved at different evolutionary distances. For each lineage comparison we clustered human ST miRNAs into genomic clusters if they were less than 10 kb apart.

Identification of orthologous miRNAs

Orthologs of human ST pre-miRNAs were identified in seven other species (mouse, dog, opossum, platypus, chicken, frog and zebrafish). We used BLAST with default parameters, with the exception of word size, which was set to six, to compare all known pre-miRNAs and mature miRNAs between all species. The best reciprocal BLAST hits highlighted miRNA pairs that were chosen as putative orthologs. In cases where we could not find orthologs at the pre-miRNA level, we used the mature miRNA instead. For precursor orthologs, we required that the length of the aligned region was equal to or greater than 56% of the query length and that the gap-size was no greater than 10.

For orthologs identified via the mature miRNA, we required that the length of the aligned region was equal to or greater than 80% of the query length and that the number of non-matching bases is was no greater than three. Furthermore, the first two to eight bases had to be conserved.

Estimating the scale of pri-miRNAs

We used the annotation produced by Saini et al. [14] to examine the distance between putative miRNA transcription start sites/transcription ends and their associated pre-miRNAs. We found that 93% of the TSSs were within 50 kb upstream of the pre-miRNAs and 90% of the miRNA transcription ends lay within 20 kb downstream of the pre-miRNAs (Figure 1). These results were used to define the general span of miRNA transcripts and extend the analysis window accordingly.

For the analysis of miRNAs with primary transcript predictions, i.e. the study of CGIs and bivalent promoters, we applied a search window that was the same size as the predictions from the study mentioned above [14]. Otherwise, the search window was 50 kb upstream and 20 kb downstream of the pre-miRNAs in the CGI analysis and only the region 50 kb upstream of the pre-miRNA in the bivalent promoter analysis. If the search window overlapped known/annotated protein-coding gene regions (the gene itself plus 1 kb up- and downstream of it) it was truncated to exclude the overlapping gene(s).

Analysis of enhancer density

Enhancer (HCNE or p300 binding site) densities were scanned using a 300 kb sliding window as described in Engström et al.[29]. For the density of the region around a sequence feature of interest (TSS for protein-coding genes), we used the HCNE density value in a 300 kb window centered on the feature.

We constructed 10000 non-coding/200 protein-coding control sets by randomly selecting the same number of coding and non-coding regions in the same genome as the conserved ST miRNAs. They had distributions of length, intergenic/intragenic ratio and gene density equal to that of the conserved ST miRNAs. For the comparison with GRB target genes, we used the whole set of GRB target genes.

Separation of ST miRNAs based on HCNE or p300 density

In the analysis of p300 density, we first ordered the human:mouse conserved ST miRNAs by human:mouse HCNE densities and mouse:human conserved ST miRNAs by p300 binding densities. Then we evenly divided human and mouse conserved ST miRNAs into two groups. The top half is located in regions enriched in HCNE/p300 binding (HCNE/p300 enriched miRNAs), whereas the remaining ST miRNAs were classified as being located in regions that are depleted in HCNE/p300 binding (HCNE/p300 poor miRNAs).

For predicting ST miRNAs that are located in regions of high HCNE density on a genome-wide scale, we ordered the conserved ST miRNAs by their HCNE densities, using HCNEs conserved between human and one of the seven vertebrate species shown in Table S4 and dividing them into HCNE enriched and HCNE poor ST miRNAs in the same way as mentioned above.

Identification of bivalent domains and bivalent promoters

Bivalent domains are regions that are enriched in both H3K4me3 and H3K27me3 histone modifications at the same developmental stage and the same tissue. We only predicted bivalent promoters for intergenic miRNAs since it was impossible to distinguish whether the features were associated with the ST miRNA or the host gene. Mikkelsen et al. [28] made the mouse histone modification information available and the human histone modification information was originally from Pan et al. [37]. We predicted bivalent promoters of protein-coding genes and miRNAs using their respective definitions.

Comparison of CpG island features

We compared the CpG island features between HCNE enriched and HCNE poor miRNAs as well as protein-coding GRB target genes, bystander genes and two other control sets, that were comprised of known transcription factors and CGI-associated genes. These gene sets were taken from the previous study of GRBs by Akalin et al. [25]. Contrary to this prior work, we only used the CpG-to-gene ratio in our analyses and not CGI length and count, since the length of the putative primary ST miRNAs is much shorter than that of protein-coding genes (p-value = 1e-03; one-sided Wilcoxon test). Furthermore, we focused only on miRNAs that have primary transcript predictions in order to ensure a high level of accuracy.