Time Series miRNA-mRNA integrated analysis reveals critical miRNAs and targets in macrophage polarization

Polarization of macrophages is regulated through complex signaling networks. Correlating miRNA and mRNA expression over time after macrophage polarization has not yet been investigated. We used paired RNA-Seq and miRNA-Seq experiments to measure the mRNA and miRNA expression in bone marrow-derived macrophages over a time-series of 8 hours. Bioinformatics analysis identified 31 differentially expressed miRNAs between M1 and M2 polarized macrophages. The top 4 M1 miRNAs (miR-155-3p, miR-155-5p, miR-147-3p and miR-9-5p) and top 4 M2 miRNAs (miR-27a-5p, let-7c-1-3p, miR-23a-5p and miR-23b-5p) were validated by qPCR. Interestingly, M1 specific miRNAs could be categorized to early- and late-response groups, in which three new miRNAs miR-1931, miR-3473e and miR-5128 were validated as early-response miRNAs. M1 polarization led to the enrichment of genes involved in immune responses and signal transduction, whereas M2 polarization enriched genes involved in cell cycle and metabolic processes. C2H2 zinc-finger family members are key targets of DE miRNAs. The integrative analysis between miRNAs and mRNAs demonstrates the regulations of miRNAs on nearly four thousand differentially expressed genes and most of the biological pathways enriched in macrophage polarization. In summary, this study elucidates the expression profiles of miRNAs and their potential targetomes during macrophage polarization.

which generates nitric oxide from L-arginine, is up-regulated along with the cytokines TNF-α and IL-1β 6 . In M2a macrophages, CD206 (mannose receptor C type 1) along with three mouse-specific genes, Ym1, Fizz1 and Arg1 are up-regulated 7,8 . However, the regulatory mechanisms that control gene expression involved in phenotypic responses are not fully understood.
Polarized M1 or M2 macrophages behave very differently in various diseases. Tumor-associated macrophages display an M2 phenotype and favor tumor growth and angiogenesis in neoplastic tissues 9,10 . Macrophage foam cells induced by western diet (high cholesterol and high fat diet) also present an M2-like anti-inflammatory phenotype 11 , contrary to previous perception that cholesterol overloading would have led to a series of M1 type pro-inflammatory responses [12][13][14] . This M2-like phenotype is linked to the accumulation of desmosterol, a natural oxysterol and endogenous ligand of Liver X Receptor (LXR) in the LXR pathway, which is well-known for the anti-inflammatory effect and the maintenance of lipid homeostasis in atherosclerosis [15][16][17] . However, it is observed in vivo that a pro-inflammatory response of the typical M1 phenotype is activated, substantially driving plaque development and progression, when macrophages within atherosclerotic lesions encounter modified cholesterol. It remains unclear exactly how the anti-inflammatory response from cholesterol overloading in the plaque macrophages is overwhelmed by the production of complex pro-inflammatory molecules that are typical of M1 phenotype.
To investigate the molecular mechanisms of macrophage polarization, we performed an integrative study of miRNA-mRNA transcriptomic dynamic changes in bone marrow derived macrophages (BMDM) under M1 and M2a conditions. We hypothesize that microRNAs (miRNAs) mediate the early events of polarity switch between M1 and M2a phenotypes, through a complex and dynamic miRNA-targeted mRNA interactome network. MiRNAs are a class of 19-24 nt non-coding RNAs that can regulate genome-wide gene expression by either destabilizing targeted mRNAs and/or inhibiting their protein translation. With advances in identification and functional annotation of miRNAs, it has been established that miRNAs play an important role in regulating immune response, especially in monocytes and macrophages [18][19][20][21] . Although several miRNA studies have been previously published investigating M1 and/or M2 macrophages 5,22,23 , these studies either only provided a "snap-shot" picture of miRNA profile at a single time point, or lacked the systematic integration of the miRNA-mRNA interactome, leading to limited discoveries.
To overcome these limitations, we designed our study using a time-series (1, 2, 4 and 8 hrs) to examine the expression profiles of paired miRNAs and mRNAs in polarized primary mouse bone marrow derived macrophages (BMDM), using next-generation sequencing platforms. Our aim was to identify: (1) the miRNAs that are robustly and differentially expressed (DE) between M1 and M2 macrophages; (2) the genes and biological functions that are targeted by DE miRNAs; and (3) whether miRNA-mRNA targetome interactions mediate the polarization of macrophage phenotypes. Our results provide detailed identification of miRNA-target regulatory networks during the early stage of macrophage polarization, which may aid in the development of future therapeutic applications for diseases involving macrophage activation.

Results
Summary of experimental design and quality assessment. To reveal the transcriptome of polarized macrophages, we sequenced paired mRNAs and small RNAs from M1 macrophages (induced by IFNγ at 10 ng/mL and LPS at 50 ng/mL) and M2a macrophages (induced by IL-4 at 10 ng/mL) at four time points (1, 2, 4 and 8 hours), with duplicates at each time point. Each duplicate was pooled from three C57BL/6 wild type mice. To ensure that the proper phenotypes were induced, we conducted qRT-PCR experiments and confirmed the expression levels of two M1 marker genes TNF-α and IL-1β as well as two M2 marker genes Arg1 and CD206 (Fig. 1).
Expression profiling analysis of miRNA-seq data. To analyze miRNA profiles in M1 and M2 macrophages, we conducted standard adapter removal and aligned the reads to mouse miRNA annotation file downloaded from miRbase (version 21) 24 , using miRDeep2 25 which enabled the discovery of novel miRNAs. We obtained an average of 7.8 million filtered reads per sample, with the average mapping rate of 77.6% (Table S1). Next we performed differential expression (DE) analysis of miRNAs between M1 and M2 macrophages, using two criteria (1) an average fold change between M1 vs. M2 over all time points greater than 1.5; and (2) statistical significance with "two-factor" (M1/M2 condition and time) matrix design using R package limma 26 , where the pairing information is considered among samples extracted at the same time.
To identify relationships between these 31 signature miRNAs, we calculated the correlations amongst them over the time course (Fig. 2B). Besides anti-correlations of M1 and M2 signature miRNAs as expected, we also discovered that M1 signature miRNAs can be divided into two sub-clusters. Sub-cluster 1 includes miR-3473b, miR-222-5p and miR-29b-1-5p, as well as the newer miRNAs miR-1931, miR-3473e, and miR-5128, all of which are tightly correlated. On the other hand, sub-cluster 2 is composed of the previously reported M1-specific miR-NAs and/or their complements, such as miR-155-3p, miR-9-5p/3p, miR-147-5p/3p, miR-125a-3p and let-7e-3p. Such expression trends show that sub-clusters 1 and 2 can be classified as early-and late-response miRNAs, respectively (Fig. S1, Table S3). Early-vs. late-response transcripts in M1 macrophages have been identified for mRNAs but not miRNAs 28 . Here, using a time-series design we were able to tease out these two types of miRNAs as well.
Expression profiling analysis of mRNA-Seqdata. To analyze mRNA profiles in M1 and M2 macrophages, we obtained an average of 30.8 million reads per sample after quality filtering and mapping rate of 96.6% (Table S4). To assess the quality of mRNA-Seq data globally, we performed PCA analysis of mRNA transcriptome in all M0, M1 and M2 samples (Fig. 3A). Consistent with earlier observations 29 , M2 is located closer to M0 on the 2D PCA plot, whereas M1 cells deviate from M0 dramatically even after 1 hour of activation (Fig. 3A). Moreover, the fold changes of specific M1 or M2 marker genes between mRNA-Seq and qRT-PCR measurements are highly correlated (correlation coefficients of 0.984-0.999). The M1 marker genes TNF-a and IL-1β and the M2 marker genes Arg1 and CD206 were used (Fig. 1), to confirm the high quality of the mRNA-Seq data.
Next we performed DE analysis on mRNA transcriptomes between M1 and M2 macrophages, similar to that for miRNAs. We identified 3937 genes as DE genes (Tables S4 and S5) over the time-course, which represents a To investigate the functional dynamic changes of DE genes, we conducted clustering analysis according to similar expression pattern using R package Mfuzz 30 , and then performed KEGG pathway analysis on each cluster. Considering expression patterns and biological functions, we separated the DE genes into 4 optimal clusters (Fig. 3C, Fig. S2 and Table S7). The genes in Clusters 1 and Cluster 2 have higher expressions in M1, compared to M2 macrophages. Cluster 1 has a decreasing trend of gene expression, with enriched pathways in signal transduction and immune response. Cluster 2 has an increasing trend and is also enriched with pathways for immune response and antigen processing and presentation. Although Clusters 1 and 2 have similar functions, their expression trends show that they are composed of early-and late-response genes, respectively. For example, the Jak-STAT signaling pathway is enriched both in Cluster 1 and 2, but it participates in both early-and late-responses by utilizing different subsets of genes. Pik3cb, Pik3r5 and interferon family members Ifna2, Ifnb1 and Ifna4 are early-response genes found in Cluster 1, while the Stat family members Stat1-Stat5, Jak2, Jak3, Akt2, Akt3, Pik3cb and Pik3cg are late-response genes found in Cluster 2.
Conversely, genes in Cluster 3 and Cluster 4 have higher expression in M2 compared to M1 polarized macrophages. In particular, Cluster 3 has high and stable expression in M2 macrophages, and is enriched in pathways involved in cell growth, DNA replication and pyrimidine metabolism. As such, it is postulated to be involved in cell proliferation and immune regulation 31,32 . On the other hand, Cluster 4 has increasing expression levels in M2 macrophages over time, and is enriched in the genes Cry2, Per1, Per2, Per3 and Nr1d1, which are involved in circadian rhythm. In summary, we have found significant differences of gene expression profiles between M1 and M2 macrophages. M1 is characterized by immune responses and signal transductions, whereas M2 is involved in basic survival mechanisms, such as cell growth and proliferation, and metabolic processes.

Functional analysis of pathways and leading edge genes targeted by DE miRNAs.
To identify the targetome of DE miRNAs, we intersected the predicted miRNA target results of targetScan with the mRNAs which have negative correlations (Spearman's correlation coefficients < − 0.5) with the miRNAs. As a result, we obtained 4464 putative targets for 31 DE miRNAs. To search for evidence of experimental validations on the 31 DE miRNAs, we checked them in two relatively comprehensive databases miRTarbase 33 and miRWalk 34 to search for evidence of experimental validation. Six miRNAs had experimentally validated targets from miRTarbase 33 and ten others had recorded targets from miRWalk 34 (Table S9). Among them, miR-9-5p and miR-155-5p were the miRNAs with the highest number (over 20) of experimentally validated targets. Fisher's exact tests on the overlapped targets between experimentally based databases and our computational prediction yielded p-values of 1.12e-38 and 3.45e-05 for miR-9-5p and miR-155-5p, respectively, indicating that our integration method yielded biological relevance.
To systematically reveal the biological processes that each miRNA is related to, we conducted Gene Set Enrichment Analysis (GSEA) on predicted targets per miRNA. We list the top ranked pathways according to the P-values (Table S10), and demonstrated the network of miRNAs and their leading edge genes of significantly targeted pathways (P < 0.05) in Fig. 4. Many miRNAs act on the same biological pathways. For example, miR-199a-5p, let-7c-1-3p and miR-26a-2-3p all target Cytokine-cytokine receptor interaction, with P-value less than 0.03. The synergy at the pathway level is largely attributed from the molecular level, where many miRNAs have common leading edge gene targets. For example, miR-155-5p, miR-155-3p, miR-199a-3p, miR-199b-3p and miR-1931 all target SMAD9. Smad9 is one of the R-smads that transduce the signals initiated from bone morphogenetic proteins in the TGF-beta (Transforming growth factor beta) pathway. The fact that multiple M1 specific miRNAs target SMAD9 may provide a molecular mechanism as to how the TGF-beta pathway is repressed in M1-polarized macrophages.

Network analysis of DE miRNAs and their targeted DE mRNAs.
Here we focused on direct and integrative targetome analysis of both DE miRNAs and targeted DE mRNAs, rather than the list of top pathways and leading edge genes targeted by individual miRNAs, as done in the previous section. Among the 3937 DE genes, 2095 (53.2%) are predicted targets of the 31 DE miRNAs. To quantitatively assess the effects of miRNAs on gene DE over time, we conducted hypergeometric tests. The results show that miRNA target genes are highly enriched among DE genes, as 29 out of 31 DE miRNAs have significantly enriched target genes (adjusted P-value < 0.05) in at least three out of four time points (Fig. 5). Moreover, among DE genes of Clusters 1 to 4, the miRNA targets account for 40.4%, 40.4%, 62.7% and 64.9% respectively, suggesting that miRNAs may play even more important roles in M2 phenotype induction compared to that of M1. Almost all 38 significant pathways in each cluster of DE genes comprise genes of miRNAs targets, except antigen processing and presentation in Cluster 2 ( Table 1). The largest percentage of miRNA targets comes from the p53 signaling pathway, where 11 out of 13 DE genes are targets of 19 miRNAs. These results provide strong evidence that miRNAs broadly and dominantly affect many of the biological processes associated with M1 and M2 macrophages.
Since the relationship between miRs and targets are complex, with one miRNA potentially targeting multiple genes and multiple miRNAs targeting one common gene, we next conducted network analysis to illustrate their complex interactions. The network consists of 10264 pairs of DE miRNAs and target genes, among which 1095 pairs (including 29 miRNAs and 864 gene targets) have more than 0.8 negative correlation that we used for further analysis. The network of miRNAs enriched in M1-polarized macrophages is composed of 22 miR-NAs and 634 genes. Of these target genes, 419 are down-regulated significantly in M1-polarized macrophages, as seen in Fig. 6A. The two hub genes targeted by the highest numbers of miRNAs are Eif2c4 and Tmem106b. Eif2c4 is presumably a target of 5 miRNAs including miR-9-5p and miR-199a-5p/3p. It encodes protein AGO4, a component of RISC that binds to short RNAs such as microRNAs (miRNAs) and represses the translation of mRNAs 36,37 . Consistent with this result, RISC was reported repressed early on during inflammation as a response to the translation of the proinflammatory cytokines during macrophage activation 38 . Tmem106b encodes Transmembrane Protein 106B, regulating lysosomal morphology and function. Previous research has shown that Tmem106a, another member of the same family, promoted an M1 phenotype through the activation of MAPK and NF-κ B pathways 39 . Additionally, among the DE targets in this highly correlated miRNA-target network, 5 immune response genes, 7 cell growth genes and 19 metabolic genes were identified. Their average fold changes in M1-versus M2-polarized macrophages is shown in Figs S3-S5. The network of miRNAs enriched in M2-polarized macrophages is composed of 7 miRNAs and 230 target genes, of which 143 are down-regulated significantly in M2-polarized macrophages (Fig. 6B). Interestingly, Cflar, Ptpn2 and Tnfsf15 are common targets of miR-27a-5p and miR-26a-2-3p; Pvr, Arid5a, Pik3r5 and Hbegf are common targets of miR-27a-5p and miR-23a-5p (in addition to Col4a3 which was identified earlier in the GSEA analysis), and Gpd2 is a common target of let-7c-1-3p and miR-25-5p. Consistent with the repression patterns in M2 in our study, these genes were largely shown to be expressed in M1-polarized macrophages. Pik3r5 (Phosphoinositide-3-Kinase Regulatory Subunit 5) and Ptpn2 (Protein Tyrosine Phosphatase Non-Receptor Type 2) are phosphatases repressed in M2. Their expression may promote M1 macrophages, since macrophage activation and polarization rely on phosphorylating and activating signaling proteins during the signaling transduction processes 40 . Other genes, such as Tnfsf15 (Tumor Necrosis Factor member 15), Hbegf (Heparin-Binding EGF-Like Growth Factor), and Cflar (CASP8 And FADD-Like Apoptosis Regulator) mainly located in the plasma membrane are involved in NF-κ B and MAPK signaling pathways, known to be active in M1 macrophages. Additionally, Pvr (poliovirus receptor) and Arid5a (AT Rich Interactive Domain 5 A) are known to be induced by LPS 41,42 . Gpd2 (Glycerol-3-Phosphate Dehydrogenase 2) is prevalent in mitochondria and possibly activates NF-κ B and MAPK by producing ROS 43 . Besides the hub genes, 12 immune response genes are identified to be repressed in M2 (Fig. S3). Overall, our results reveal that miR-27a-5p, miR-26a-2-3p, let-7c-1-3p, miR-25-5p and miR-23a-5p boost M2 phenotype through repressing multiple genes related to NF-κ B and MAPK signaling pathways.

Discussion
The aim of this study was to fine map the expression of miRNAs in macrophage polarization, and provide the first interactome of miRNAs and their targetome. To achieve this, we sequenced paired mRNAs and miRNAs from primary mouse bone marrow derived M1-and M2-polarized macrophages with up to 8 hours of induced polarization. We identified 31 robust DE miRNAs that provide signatures of M1 and M2-polarized macrophages. These miRNAs exhibited substantial regulation during macrophage polarization, correlating with a majority of the nearly four thousand DE genes identified, as well as nearly all of the biological pathways analyzed. Our time series analysis has revealed that not only mRNAs, but also miRNAs, can be categorized into earlyand late-response elements. The early-and late-response miRNAs present in M1 macrophages exhibit strong positive correlations between members within the same group. On the other hand, early-and late-response mRNAs are present in both M1-and M2-polarized macrophages. During M1 polarization, the early-and late-response DE mRNAs are involved in immune response/signal transduction through different functional gene subsets. In M2 macrophages, early-response DE mRNAs are involved in cell growth/metabolism, while the late-response DE mRNAs have enriched function relating to circadian rhythm. M1 macrophages are involved in proinflammatory activation while M2 macrophages engage in cell proliferation and metabolism.
We also observed several discrepancies between our study and others. For example, miR-27a-5p was reported higher in M1-polarized macrophages by Graff JW et al. 22 , however our study shows that miR-27a-5p is significantly up-regulated in M2 macrophages. The differences could be due to several reasons, such as the time and dosage difference of the treatment. Although both studies used identical conditions to induce an M1 phenotype, we used a lower dose of IL-4 at 10 ng/mL, compared to their 20 ng/mL to induce the M2a phenotype. Also, we profiled the gene expression in earlier time points (up to 8 hrs), while Graff JW et al. investigated the prolonged treatment effect after 72 hrs 22 . Thus it is possible that there is an initial surge of miR-27a-5p M2 but then its expression is higher in M1 after prolonged treatment.
One of the most interesting findings of this study is that a relatively large number of C2H2 zinc-finger genes among the transcription factors targeted by miRNAs, largely due to sequence similarities in their 3′ UTR regions. The function of these zinc-finger proteins in DNA binding and immune response regulation 46,47 , were previously reported to be post-transcriptionally regulated by miR-23, miR-181, miR-188, and miR-199 through seed-match in conserved coding regions 48,49 . Here we have found that they are targeted by 29 out of 31 DE miRNAs. The functions of C2H2 family are diverse, including both transcriptional activators and suppressors, although suppressor functions are more often found 50 . The majority of the C2H2 genes were shown to be repressed following LPS stimulation 51 , consistent with our results. Future research may elucidate the exact mechanism of miRNAs regulation of the C2H2 zinc-finger family.
In order to validate the 3 newly identified M1-early response miRNAs (miR-1931, miR-3473e, and miR-5128) and test their functional relevance in macrophage polarization, we performed experiments using mimics and inhibitors of each miR in primary BMDM. Our discovery that Arg1 is a bona fide target of miR-5128 implicates that miR-5128 inhibition may be beneficial in treating pathologies involving macrophage inflammation. While Arg1 is a predicted direct target of miR-5128 according to Targetscan 52 , IL-10 does not show complementary binding sites for miR-5128. Thus, it is likely that IL-10 is an indirect target of miR-5128. Although further investigation is needed to more thoroughly identify the regulatory mechanism of the 3 new miRNAs on target genes, our findings set the groundwork for further exploration of the biological functions of novel miRNAs in macrophage polarization. Based on our findings, regulation of miR-5128 appears to be a promising measure for inducing an anti-inflammatory macrophage phenotype. Other important points of discussion concern the complex spectrums of M1 and M2 macrophage phenotypes. We have followed the experimental conditions widely accepted in the field of macrophage biology and well documented by others 53 , with IFNγ and LPS for M1 polarization and IL-4 for M2a phenotype induction. IFNγ and LPS may induce slightly different but largely overlapping changes in gene expression in M1 29 . One obvious caveat to our study is that the in vitro induction of macrophage polarization cannot be extrapolated to in vivo conditions. Despite the limitation of our work (similar to others published within the field), we believe that our data and bioinformatics analysis are the most complete thus far in terms of the involvement of miRNAs and their targets in macrophage polarization.
In summary, we have deciphered functional miRNAs in the polarization of macrophages through the repressed targets, including immune response genes, NF-κ B signaling genes and phosphatase genes as well as cell cycle and metabolic genes. Additionally, our research has shown that C2H2 zinc-finger family members are the majority of targeted transcriptional factors by DE miRNAs.

Materials and Methods
The authors confirm that all methods were carried out in accordance with relevant guidelines and regulations, and that all experimental protocols were approved by University of Hawaii's Institutional Biosafety Committee. M1 and M2 primary macrophage cell culture and RNA extraction. Bone marrow cells were extracted from both femurs and tibias of 6-8 week old male C57/BL6J mice. The isolated cells were strained through a 40 μ m cell strainer, then plated at 25 × 10 6 cells per 15-cm plate for macrophage differentiation in DMEM supplemented with 10% FBS, 20% L929-conditioned media and 1% penicillin/streptomycin. The use of L-929 conditioned medium as a source of M-CSF in primary murine BMDM culture has been standard protocol and has been described in detail previously 53 . Adherent macrophages were replated on day 7 into 6-well plates at 3 × 10 6 cells per well in medium without L929 conditioned media supplementation. After overnight incubation, the macrophages were treated with 1) IFNγ (10 ng/mL) + LPS (50 ng/mL) or 2) IL-4 (10 ng/mL) in duplicate for 1 hr, 2 hr, 4 hr or 8 hr, in order to induce M1 or M2 phenotype, similar to previous published studies 22,29 . Untreated macrophages in duplicate were used as controls and collected at the 1hr time point. Cells were collected and centrifuged at 400 × G for 5 minutes. Cell pellets were used for total RNA extraction with an miRNeasy mini kit (Qiagen), designed to purify genomic total RNAs including mature miRNAs. The quality and quantity of the nucleotide samples were assessed using the Agilent Bioanalyzer and Nanodrop technologies (University of Hawaii Genomics Shared Resources Facility).
Quantitative Real-Time PCR (qRT-PCR) validation of marker genes. In order to verify the expression data generated by mRNA-seq, we performed qRT-PCR experiments for polarization marker genes. Total RNA was isolated with TRIzol (Invitrogen). cDNA was synthesized from 1 μ g of total RNA in a 20-μ l reaction with qScript cDNA Synthesis kit (Quanta Biosciences). qRT-PCR for the following genes: TNFα , IL-1β , Arg1 and CD206 was performed for each sample in triplicate and Gapdh served as internal control. qRT-PCR was performed using the 2x FastStart Universal SYBR Green Master Mix (Roche) and primers were purchased from Integrated DNA Technologies (IDT). Samples were run on an Applied Biosystems 7900HT Fast Real-Time PCR System. Read mapping. Small RNA sequencing and mRNA sequencing were done in Yale Keck Genomic Cores, both at 75 bp read-length with single-end reads. For RNA-seq data, raw reads were removed adapters and aligned to UCSC mouse mm10 genome using Tophat (v2.0.8b) 54 . Default parameters were used (2 base mismatches allowed in the read, and 0 allowed in the splice regions). HTSeq 55 was used to calculate the raw count numbers based on the mm10 gff file. For miRNA-seq data, raw reads were calculated by miRDeep2 25 by mapping to mouse miRNA annotation file (with 1280 miRNAs) downloaded from miRbase (version 21) 24 .
Differential expression analysis. The raw counts of mRNA-and miRNA-seq were normalized between samples using DEseq. A total of 16531 genes were detected in M1 or M2 macrophage. In miRNA annotation file, A total of 790 miRNAs were detected in M1 or M2 macrophage, among which 507 miRNAs had maximum count numbers over 10 reads after normalization. The 16531 genes and 507 miRNAs were subject to differential expression analysis using the R package limma 26 , with two-factor design to consider both time and sample conditions as two factors. For RNA-seq data, genes were considered DE when the averaged absolute fold-changes at four time points were greater than 2 and the adjusted p-values were less than 0.05. For miRNA-seq data, miRNAs were considered DE when the averaged absolute fold-change at four time points were greater than 1.5 and the adjusted combined p-values less than 0.05.
Time series clustering of differentially expressed genes. The R package Mfuzz (v2.30) 30 was used to cluster DE genes along time series. Mfuzz is a soft-clustering method based on fuzzy c-means algorithm. Average expression values at each time point were log2 transformed as the input to generate clusters based on the expression trend. The genes were assigned to each cluster with membership values over 0.5 and 98.86% (3892/3937) belonging to a unique cluster.
MiRNA target prediction. MiRNA target prediction was based on Targetscan 52 . Mouse UTRs were downloaded from the targetscan website, and mature sequences of miRNAs were downloaded from miRBase 24 . Pre-computed miRNA targets were downloaded from the targetscan website. For the miRNAs that were not in the pre-computed summary file, the targets were predicted using a perl script on the Targetscan website. Additionally, we used the correlation criteria to refine the predicted miRNA targets. we calculated the spearman correlation coefficient (SCC) in each miRNA-target pair across all samples, and chose miRNA-target pairs predicted by tar-getScan with SCC < − 0.5 and adjusted P-values of SCC statistical test < 0.05.
Network and function enrichment analysis. DAVID (v6.7) 56 was used for KEGG analysis on DE genes, with Benjamini-Hochberg adjustment for p-values. Gene Set Enrichment Analysis (GSEA) 57 was used to obtain significantly enriched pathways and the leading genes targeted by each miRNA. Cytoscape 3.0 58 was used to display the miRNA-target network. The protein-protein interaction downloaded from NIA Mouse Protein-Protein Interaction Database 59 was used to infer the network among targets, with specification of mouse TF data downloaded from AnimalTFDB 60 . High-degree targets were obtained through network analysis in Cytoscape.
Transfection with microRNA inhibitors and mimics for functional validation. Primary mouse bone marrow derived macrophages were isolated and differentiated as described above and 3 × 106 cells per well of a 6-well plate were seeded for experiments overnight in DMEM supplemented with 10% FBS and 1% pen/strep. miRCURY LNA TM miRNA inhibitors or mimics for mmu-miR-1931, mmu-miR-3473e, or mmu-miR-5128 were purchased from Exiqon. X-tremeGENE siRNA transfection reagent (Roche) was used to transfect the cultured macrophages with the miRNA mimics or inhibitors at a final concentration of 100 nM in 2 mL of fresh media for 24 hours. Following transfection, the cells were washed once with PBS and miRNA inhibitor-transfected cells were treated with IFNγ (10 ng/mL) + LPS (50 ng/mL) for 2 hours, while miRNA mimic-transfected cells were treated with IL-4 (10 ng/mL) for 4 hours.