MIREyA: a computational approach to detect miRNA-directed gene activation

Introduction

Conventionally, microRNAs (miRNAs) are considered to suppress gene expression through RNA interference (RNAi) by binding complementarily to mRNAs, forming a RISC complex, and causing RNA degradation (reviewed in detail in).¹ However, recent studies provide evidence that some miRNAs act in the opposite way – stimulating gene activation. Numerous studies have demonstrated the ability of miRNAs to up-regulate genes by targeting their promoters.^2,3 Ago1 from miRNA-Ago complex associates with the Ccnb1 promoter and miR-744 induces enrichment of RNA Pol II and H3K4me3 at the Ccnb1 transcription start site.³ MiRNA let-7i interacts with the TATA-box of IL-2 gene and stimulates transcription initiation by contributing to the assembly of the pre-initiation complex.⁴ Relatively fewer miRNAs demonstrated the ability to unconventionally target and activate enhancers, thus inducing genes regulated by these enhancers. MiR-24-1 acts as a modulator of the chromatin state of an enhancer. Furthermore, it increases p300 and RNA Pol II binding at the enhancer locus. The miR-24-1 actually originates from the enhancer locus. However, some genes regulated by other enhancers are also expressed at higher levels when miR-24-1 is transfected, and the enhancers of induced genes contain a sequence similar to the seed of the miRNA.⁵ These observations suggest that other miRNAs might trigger enhancers and activate gene expression.

When miRNAs function as activators in a nucleus, different targeting mechanisms are possible: miRNA:DNA Watson-Crick duplex formation as well as miRNA:DNA Hoogsteen triple helix formation. Nuclear miRNA target prediction tools utilize are based on the idea that miRNA:DNA interaction requires an intact seed region.⁶ MicroPIR2 predicts targets in mouse and human promoter regions.⁷ Trident predicts miRNA:DNA Hoogsteen-type base pairings.⁸ Some tools designed to predict conventional miRNA targets may also be utilized to find nuclear activational targets, e.g. miRanda.⁹

In this work, we report MIREyA (MIRnas functioning through Enhancer Activation), a pipeline for detection of miRNAs and their gene targets up-regulated through triggering their enhancer in the provided expression dataset. We applied MIREyA in order to identify and characterize activational miRNAs in Mycobacterium tuberculosis (Mtb) infected macrophage dataset. MiRNAs are important regulators in macrophage responses during Mtb infection that may act as host immunity agents as well as a tool exploited by pathogen agents to manipulate host cell pathways.¹⁰ Yet, only the classical role of miRNA has been investigated so far in the context of Mtb infection.^11–19 Although multiple studies have shown the possibility of the activational role of miRNAs, this potential remains neglected in the study of miRNAs in bacterial infections. MIREyA found several miRNAs, which have not been shown to be functional in TB yet, suggesting it could be useful to find candidate activational miRNAs for further experimental validation. Mir155 has previously been shown to act as a negative regulator of essential mRNAs during TB,^20,21 but not as an activator.

Methods

Main steps of the algorithm

MIREyA aims to detect miRNAs with the potential to upregulate a gene via activation of its enhancer. It consists of three major steps:

• 1) The algorithm detects miRNA bound to the enhancers. This step is implemented with three different approaches described below.

• 2) The algorithm selects genes regulated by the enhancers selected in step 1). For this step the output from the first step is required as well as a table with enhancer:gene pairs where an enhancer is assumed to regulate the corresponding gene.

• 3) The algorithm calculates the Spearman’s correlation coefficient (SCC) between the expression levels of miRNAs and genes regulated by corresponding enhancers selected in step 1), and estimates the p-value of the SCC with a Benjamini-Hochberg correction using the number of the miRNA:gene pairs for one miRNA (FDR < 0.05). The input data is gene expression data with sample size ≥ 8 and the output of the previous step.

Figure 1 illustrates the workflow of our algorithm implemented in Python and R.

Figure 1. The outline of the implemented algorithm.

Use case

We applied MIREyA to three time-series (0, 4, 12, 24, 48, 96 hours) expression datasets (CAGE) of mouse bone marrow-derived macrophages infected with hypervirulent Beijing/W lineage Mycobacterium tuberculosis (Mtb) HN878 strain, 2-3 replicates per time point.²⁴ Each dataset corresponds to the time series after infection for macrophages of different phenotypes: not pre-stimulated (M0), interferon-γ stimulated (M1-polarized) and interleukin-4/interleukin-13 stimulated (M2-polarized). Only differentially expressed (DE) miRNAs and genes were considered. We obtained enhancer-gene interactome from²⁵ where an enhancer is predicted to regulate a gene if their expression levels correlate significantly and they belong to the same topologically associated domain (TAD).

We searched for candidate enhancers targeted by miRNAs with all three methods described previously and merged the results for further steps. CAGE enables one to estimate expression at the promoter level, while enhancers are associated with whole genes. As a proxy of gene expression, we used either an expression value of a promoter with the highest SCC (a) or summed up the expression values of all promoters of the gene (b). To reduce the number of false positive predictions, we selected among miRNA:enhancer duplexes only such cases where (1) a duplex was predicted both by miRanda and Needle-based approach; (2) an enhancer was associated with several genes since in the original paper on miRNA-activated enhancers⁵ one nuclear miRNA affected expression of multiple genes regulated by the triggered enhancer; (3) the miRNA-gene pair was obtained in both ways (a and b) to estimate expression of the gene. We also added miRNA-gene pairs with highly correlated expression levels (SCC ≥ 0.8) in any of four combinations of two methods to detect miRNA:enhancer interactions (seed match + Needle and miRanda) and two approaches to treat expression of different promoters (a and b). Among predicted miRNA:enhancer triplexes we selected cases where (1) an enhancer was associated with multiple genes; (2) both approaches to estimate gene expression (a and b) yielded this triplex.

Results

We applied MIREyA to three time-series datasets of Mtb-infected macrophages with 3 different phenotypes prior to infection. In M0 macrophages 10 miRNAs were differentially expressed (DE) in at least one time point compared to the state before infection (0h), in M1 there was no DE miRNA, while in M2 only Mir1956 was DE. Figure 2 and Table S1 (Extended data⁴³) represent detected miRNA-enhancer-gene trios for M0 macrophages (Extended data: Table S2 for M2 macrophages⁴³). We investigated the functions of the obtained genes and miRNAs and confirmed that many of them might be involved in the response to Mtb infection. We detected Mir155, which is vastly studied, and known to subvert autophagy in human dendritic cells²⁰ and to be a potential diagnostic marker of active tuberculosis.²¹ Other miRNAs and their targets which we consider promising for further investigation are summed in Tables 1 and 2. Mir22, Mir221, Mir222 are annotated with high confidence.²⁶ Expression of only some promoters of Klf6 and BC016423 genes correlates with Mir22 expression, but we included them, because they have predicted triplexes with three different enhancers.

Figure 2. Predicted regulation graphs of miRNAs via RNA:DNA duplex (A) or RNA:DNA triplex formation (B) with an enhancer for M0 macrophages.

The length of an edge between a miRNA and an enhancer or between an enhancer and a gene is in inverse proportion to the average correlation between the expression of the miRNA and the target gene or the enhancer and the gene respectively. Full coordinates of enhancers are available in Table S1.

We reconstructed regulatory networks for these miRNAs based on miRNA-enhancer-gene trios and investigated their potential role in the response to Mtb infection. Mir22-activated gene network is highly likely to be involved in Mtb response. Klf6, a potential target of Mir22, is a transcription factor essential for macrophage motility²⁷ and plays an important role in the regulation of macrophage polarization promoting M1 phenotype cooperatively with NF-κB.²⁸ Nfkb1 gene, a putative target of miRNAs Mir221, Mir222 and Mir155, encodes a subunit of NF-κB protein complex, a master transcription factor in macrophage immune responses. The human ortholog of Ube2d3, potentially regulated by the same miRNAs, facilitates polyubiquitination of NFKBIA (a member of the NF-kappa-B inhibitor family) stimulating its subsequent degradation.²⁹ A detected target of both Mir22 and Mir221, Peli1 regulates the NF-κB activity negatively and attenuates the induction of proinflammatory cytokines in T-cells.³⁰ Cxcl1 and Cxcl2, detected targets of both Mir22 and Mir221, are chemoattractants for neutrophils contributing to tissue inflammation.³¹ Malt1, potentially up-regulated by Mir155, is known to activate NF-κB in lymphocytes.³² Among targets of Mir1956 detected in M2 macrophages dataset we found the Ccrl2 gene encoding a chemokine receptor-like protein which is expressed at high levels in primary neutrophils and primary monocytes. Another Mir1956 target, Dotl1, is an H3K79 methyltransferase which facilitates the expression of IL-6 and IFN-β in macrophages.³³ Cd14 leads to NF-κB activation and inflammatory response,³⁴ Cd14 KO mice infected with Mtb are protected due to reduced inflammatory responses at the chronic stage.³⁵ Rab20 plays a role in the maturation and acidification of phagosomes and the fusion of phagosomes with lysosomes during mycobacterial infection.³⁶ Ticam1 is involved in native immunity against pathogens: it interacts specifically with toll-like receptor 3, activates NF-κB.³⁷ Tnfaip3 (A20) is an important regulatory protein that down-regulates NF-κB activity.³⁸

We further investigated protein networks associated with the selected miRNAs. Proteins regulated by miRNAs Mir22, Mir155, Mir221, Mir222 are involved in regulation of NF-κB – a vital orchestrator of the response of the innate immune cells to pathogens³⁹ (Figure 3A). Cxcl1 and Cxcl2 regulated by Mir22 and Mir155 (Figure 3B) are chemokines which signal through CXC receptor 2 to attract neutrophils to the place of inflammation, which is essential to control tissue infection.³¹

Figure 3. MiRNA-enhancer-gene trios detected for M0 macrophages include genes with related functions and form subnetworks A (related to NF-kB) and B (both genes encode chemokines which signal through CXCR2).

Relationships are depicted with arrows of a different colour: green – activation, red – repression. Dashed arrows indicate predicted relationships, solid – known from published studies and described previously. Full coordinates of enhancers are available in Table S1.

Discussion

MIREyA aims to easily find candidate activating miRNAs which trigger an enhancer and genes up-regulated by the enhancer. New emerging methods for studying RNA–DNA interactome detect previously unknown miRNAs bound to chromatin^40,41 which might be promising for further experimental investigation in terms of understanding complex gene regulation networks in detail. Although high-throughput data on RNA:DNA interactions in several cell types but not macrophages are available now, we could not confirm MiRNAs detected with MIREyA for Mtb infection using either RADICL-seq⁴¹ or GRID-seq,⁴² which is unsurprising since the ncRNA:DNA interactions are highly cell-type-specific.⁴¹

Although the implemented algorithm considers multiple aspects of the suggested mechanism of gene activation, several factors remain unaccounted for. To run MIREyA a user requires a priori knowledge of enhancer-gene interactions. A common fast computational approach to determine gene-enhancer pairs is to use genes and enhancers co-localized in the genome, ignoring the long-distance spatial interactions. One of solutions is the approach suggested in:²⁵ to calculate the correlation between expression of genes and enhancers which belong to the same TAD.

Although the exact mechanism of RNA activation of enhancers remains unclear, we do know that miRNA as a mediator facilitates epigenetic modifications in the enhancer region. At this stage, MIREyA does not consider chromatin availability or any other epigenetic information. In order to reduce the number of false positive predictions, nuclear localisation of specific predicted miRNAs should be validated experimentally.

Despite the discussed limitations, using MIREyA we detected several promising miRNA candidates. We suggest that MIREyA provides a promising approach to select miRNAs which up-regulate genes by triggering enhancers for further experimental validation.

Conclusion

Our method extends the study of activational miRNAs and provides a basis for further research. The use case on Mtb-infected macrophages demonstrates the possibility of existence of novel miRNAs up-regulating gene expression.

Data availability

Software availability

Pipeline available from: https://github.com/veania/MIREyA.

Archived pipeline as at time of publication: https://doi.org/10.5281/zenodo.4635578.⁴⁴

License: Open Software License 3.0