### Overview of finding co-regulating miRNAs

Steps for identifying co-regulating miRNAs are presented in Figure 1. Co-regulating miRNAs correspond to a group of miRNAs silencing together target genes. Sequence complementarity-based computational prediction can consider static target relation only but not cooperative binding of miRNAs that may occur in a condition-specific manner. One needs to systematically investigate expression patterns between miRNAs and mRNAs across different conditions. We used miRNA-mRNA co-expression profiles to select condition-specific modules of cooperative miRNAs at the level of functional expression. Using miRNAs which share common targets, we constructed MRAN and extracted clusters of tightly co-regulating miRNAs from MRAN based on the network analytic approach.

### Datasets

We first used miRNA-mRNA co-transcriptomic datasets obtained from the same samples by ourselves in our previous studies (Table 1). Circadian rhythm dataset [12] consists of triplicated mRNA (Illumina Mouse-6 Expression BeadChip) and duplicated miRNA (Ambion mirVana™ miRNA Bioarray V2) expression profiles from 12 time points with four-hour interval across two complete circadian cycles, *i.e.*, 48 hours, resulting 36 and 24 hybridizations for mRNA and miRNA, respectively. One can download the dataset at Gene Expression Omnibus (GSE11516). Radiation dataset [36] consists of triplicated mRNA (Affymetrix Human Gene 1.0 ST Array) and duplicated miRNA (Ambion mirVana™ miRNA Bioarray V2) expression profiles from 6 time points, resulting 18 and 12 hybridizations for mRNA and miRNA, respectively. Two lung cancer cell lines, H460 and H1299 were irradiated at 2 Gy, and harvested after 0, 2, 4, 8, 12 and 24 hours to examine the expressions. Actinomycin D dataset (unpublished data) is created to explore the bio-molecular mechanism of miRNA decay. It consists of triplicated mRNA (Affymetrix Mouse Gene 1.0 ST Array) and duplicated miRNA (Ambion mirVana™ miRNA Bioarray V2) expression profiles from 7 time points, resulting 21 and 14 hybridizations for mRNA and miRNA, respectively. N2a mouse neuroblastoma cells were treated with Actinomycin D to block transcription at 0, 1, 2, 4, 8, 12 and 24 h time points. We also included external dataset from 57 prostate adenocarcinomas (GSE7055) [37] (Table 1).

Normalizations of the miRNA and mRNA expression profiles were performed separately. For datasets using Affymetrix GeneChips, we used RMA (Robust Multi-Array) normalization method [38]. For datasets using Ambion miRNA chips and OSU-CCC chips, we used vsn transformation after background subtraction [39] and applied quantile normalization method [40]. For datasets using Illumina chips, we used chip-wise method using the rank invariant algorithm [41].

### Target overlap score and MRAN construction

MRAN is defined as a combined network of cooperative miRNAs sharing condition-specific target genes. A node in MRAN corresponds to a miRNA and an edge corresponds to a condition-specific cooperativity between a miRNA pair. Condition specific cooperativity is defined as condition specific target sharing and determined by (static) significant sequence complementarity by computational target prediction algorithms provided by TargetScan and (functional) significant inverse expression pattern at a specific experimental condition (by Pearson's correlation coefficient). We used the Pearson's correlation coefficient (PCC) for measuring similarity/dissimilarity between expression patterns of miRNAs and genes. The PCC of a pair of genes commonly returns a real value in [+1, -1]. PCC(x, y) > 0 represents that x and y are positively correlated with the degree of correlation. On the other hand, PCC(x, y) < 0 represents that x and y are negatively correlated with a value |PCC(x, y)|. A positive value of the PCC indicates that miRNAs and genes are co-expressed and a negative value of the PCC indicates that opposite expression pattern exists between them. Considering the conservation of miRNA-target relationship, we retrieved and applied 50,339 and 50,349 human and mouse miRNA-target pairs for the 162 and 7,927 conserved miRNAs and mRNAs by using TargetScan target prediction database (release 4.2) [42].

Target overlap score between a pair of miRNAs is defined as Jaccard similarity coefficient, representing the fraction of shared condition specific targets. We define the target overlap score

*t*
_{
ij
} between miRNAs

*i* and

*j*as follows:

where *Targets(i)* represents to the set of targets of miRNAs *i*. We measured the statistical significance of the Jaccard similarity coefficient by using the exact randomization tests [43]. The *n* × *n* condition specific target overlap matrix *T* = [*t*
_{
ij
}] is transformed into an *n* × *n* adjacency matrix *A* = [*a*
_{
ij
}], which encodes the connection strengths between pairs of nodes based on the probabilities related to the Jaccard similarity coefficient. Since the networks considered here are undirected, *A* is a symmetric matrix with non-negative entries.

Figure 1 outlines the steps of the present study. Figure 1(B) conceptualizes (1) the miRNA-mRNA target relationship as a directed bi-partite graph at the sequence level. (2) Many condition-specific directed bi-partite graphs are emerging by combining inversely co-expression patterns between miRNA and mRNA pairs. (3) By means of applying target overlap scoring, many condition-specific MRANs are obtained as undirected weighted graph of miRNAs. (4) Many condition-specific MRANs are merged into a combined multiple graph to form global MRAN (see Figure 2(A)). Now one can extract cooperative miRNA modules (CMMs) by applying subnetwork detection algorithms and then characterize and evaluate them in terms of biological relevance.

### Extracting cooperative miRNA modules

In the present study, a subnetwork detection algorithm, called MCODE [12], was applied to detect coherent groups in the global MRAN. The MCODE is a graph theoretic clustering algorithm specifically designed to find complexes by identifying densely connected subgraphs in networks. MCODE algorithm consists of three stages: vertex weighting, complex prediction and an optional post-processing step. The weighting of nodes is based on the core clustering coefficient. Once the weights are computed, the algorithm traverses the weighted graph in a greedy fashion to isolate densely connected regions. The post-processing step filters or adds nodes based on connectivity criteria. MCODE has a parameter that specifies the size of clusters returned. All MCODE parameters are applied with default values. Subnetworks are filtered if they do not contain at least a 2-core (graph of minimum degree 2). This approach allows us to assign one miRNA to multiple clusters, considering a biological principle that miRNAs can have multiple functions.

### Regulators of cooperative miRNAs

To investigate the modular nature of global MRAN, we used predicted binding sites for all Position Specific Scoring Matrices (PSSMs) from TRANSFAC version 8.3, as they are defined by the UCSC hg17 genome assembly, in the tfbsConsSites (

http://genome.ucsc.edu/) and tfbsConsFactors. All RefSeq genes genomic locations were taken from hg17. In mouse, regulators of miRNAs were identified in the human genome in the regions orthologous to the mouse. For example, hsa-miR-101 in human and mmu-miR-101 in mouse are orthologous. To assess the statistical significance of the rate of co-occurrence of motif pairs, we used cumulative hypergeometric distribution to calculate the probability of obtaining the rate of co-occurrence,

*C*, equal to or higher than the observed rate of co-occurrence,

*c'*, by chance:

where *m*
_{
1
} and *m*
_{
2
} denote the numbers of miRNAs having each of the two motifs, *N* denotes the total number of miRNAs in the genome, and *i* the summation index.

In the co-occurrence graph of regulators on MRAN, nodes represent transcription factors (TFs), and edges statistical significance of co-occurrences. Hypergeometric *p*-values were used as weights (see Figure 2(B)).

### Functional coherence analysis of target genes

Plausible characteristics of an extracted subnetwork as a CMM might be functional coherence, meaningfulness as well as distinctness from other modules. We performed functional enrichment analysis by using KEGG (Kyoto Encyclopedia of Genes and Genomes) [44] and BioCarta (http://www.biocarta.com) pathways to test their functional coherence and meaningfulness. KEGG pathways mainly include cellular processes related to metabolism and biosynthesis. Those on BioCarta cover a wider variety of cellular processes including a large number of signal transduction and immune signalling pathways as well as metabolic and biosynthetic pathways. Pathways with *p*-value < 0.05 as revealed by the hypergeometric test were considered statistically significant in the present study. To visualize the distinctness of the cooperative modules in terms of functional enrichment, we created heat maps of the significances (see Figure 4).

### Sequence/structure analysis of the co-regulating miRNAs

Structural similarity of the members of CMMs was evaluated. T-COFFEE (version 5.05) [45] was used to calculate pairwise sequence similarity between miRNAs. It uses information from a pre-compiled library of different pairwise alignments including both local and global alignments. RNAdistance program of the Vienna RNA package (Version 1.7.1) [46] was applied to calculate structural distances. Sequence-structure based clustering using the LocARNA-RNAclust pipeline [47] was applied to pre-miRNA sequences. LocARNA uses a complex RNA energy model for simultaneous folding and sequence/structure alignment of the RNAs [47]. The resulting alignment scores can be used to cluster RNAs according to their sequential and structural similarities.