We used the human metabolic pathways from KEGG database (ftp://ftp.genome.jp/pub/kegg/xml/organisms/hsa) to construct the metabolic network. We identified all the reactions and their associated enzymes and genes. We used a reaction-centric model to represent the metabolic network. The network is presented as a directed graph in which nodes represent enzymes or reactions. For two nodes (reactions) A and B, we linked A → B if a product metabolite of A is a substrate metabolite of B.

The most commonly used molecules such as ATP, NAD, CO₂ in biochemical reactions were excluded in the network. For each node in the metabolic network, we identified all the genes associated with it. Human metabolic network contains 1,549 nodes and 4,740 directed links. In the network, over 70% of the nodes are linked to form the largest network component. The largest component of the human metabolic network contains 1,099 nodes and 4,196 links. In this study, we used only the largest contiguous component of the network to perform analyses.

We used genome-wide computationally predicted miRNA target genes in human genome from a recent study [9]. Updated miRNA targets from the prediction were downloaded from authors’ website.

In a reaction-centric network, nodes represent enzymes or reactions. We classified the nodes into five categories based on their structural features in the network. The nodes that have no incoming links are called upstream nodes (UPNs, Fig. 1, Supplementary Data 1), where the metabolites are often up-taken from the extracellular space, are used to produce downstream metabolites of the network. Therefore, the UPNs are the most upstream reactions that many other downstream reactions and metabolites rely on. On the other hand, the nodes that have no outgoing links are called output nodes or downstream nodes (DSNs, Fig. 1, Supplementary Data 1), where many pathway output metabolites such as amino acids of the network are produced. These metabolites are then used for various cellular activities such as acting as cell membrane components or regulating signaling pathways. A cut vertex or a cut point (CP, Fig. 1, and Supplementary Data 1) is such a bottleneck node that its deletion will disconnect at least one component from the network. CPs are in crucial network positions and become bottlenecks of the network, and therefore control metabolic flows from a part to another in the network. The removal of one node from a network results in more network components than in the original network, we then regarded that node as a CP. To identify CPs, we examined each node in the network to check whether its removal results in more network components. In the network, highly connected nodes (hubs, Fig. 1, Supplementary Data 1) are the reactions that share metabolites with many other reactions. The metabolites of these highly connected nodes are also major suppliers for the peripheral region of the network where many output metabolites are produced. In this study, we defined the top 5% of the highly connected network nodes, which have higher degrees (in- and out-degrees together) as hubs. Besides these four types of nodes described above, all other nodes in the network are called intermediate nodes (ITNs, Fig. 1, and Supplementary Data 1).

Fig. (1) Node types with distinct network structural features.

To test the statistical significance of observations, we performed randomization tests. A more detailed explanation of randomization tests was previously described by Wang and Purisima [10]. Briefly description of the randomization tests was the followings:

To test the statistical significance of miRNA regulation of a multiple-gene-node, we calculated the fraction of miRNA targets of the genes which are associated with the multiple-gene-node. We then randomly shuffled the miRNA targets (genes) 5,000 times among the metabolic genes and calculated the P values. We further corrected the P values by applying False Recovery Rate (FDR). FDR corrected P values (<0.25) were considered as significant.

To test the statistical significance of the node types which are either enriched or less-enriched with miRNA targets, we calculated the fraction of miRNA targets of each node type, and then randomly shuffled the miRNA targets (nodes) 5,000 times among the network nodes to calculate the P values.

To test the statistical significance of the linear network motifs which are enriched with miRNA targets, we calculated the fraction of the linear motifs in which all the nodes are miRNA targets, and then randomly shuffled the miRNA targets (nodes) 10,000 times among the network nodes to determine the P values.

Metabolic flows are arranged as linear or branching reactions. To test the statistical significance of the miRNA regulatory patterns on metabolic branches, we enumerated all the possible miRNA targeting patterns for each branch type (see Table 1) then calculated the fraction of each pattern. Finally, we randomly shuffled the miRNA targets (nodes) 10,000 times among the network nodes to determine the P values.

Table 1

Patterns of miRNA Targeting in Convergent and Divergent Units

We obtained human metabolic pathways from the KEGG database. We conducted the analysis by considering two situations: (a) miRNA-target enrichment in each individual pathway; (b) in a network context, a pathway that contains miRNA-targeted downstream nodes (output nodes, DSNs) or miRNA-targeted CPs who are the parents of the pathway’s DSNs. To determine the statistical significance of the enrichment of miRNA targets of a pathway, for each pathway we first calculated the number of the reactions, which are miRNA-targeted nodes that have been determined in this study, and then randomly swapped the miRNA targets among the metabolic reactions (nodes). In each swap, we recalculated the fraction of the miRNA targets in each pathway and compared to that of each pathway. The P values were corrected by FDR. FDR corrected P values (<0.25) were considered as significant. To examine the pathways that are regulated by miRNA in a network context, we considered: (a) if a DSN of a pathway is the target of miRNA, we classified that pathway as miRNA-regulated; (b) when a miRNA targets a CP which is a parent of a DSN of a pathway, we also classified that pathway as miRNA-regulated.

To further explore miRNA regulatory principles, we asked if miRNA targets are enriched in certain node types, i.e. hubs, UPNs, DSNs CPs and ITNs. To address this question, we calculated the fractions of miRNA targets for each node type and performed randomization tests. To perform the randomization tests, we randomly assigned 238 targets onto the network nodes and examined the statistical significance of the miRNA targets in each node type (see Methodology). As shown in Table 2, miRNA targets are less enriched in ITNs in the network (P=0.032, Table 2, Supplementary Data 4), suggesting that miRNAs avoid regulating ITNs. On the other hand, miRNA targets are significantly enriched in the hubs and CPs of the network (P=0.015 and P=0.006, respectively, Table 2, Supplementary Data 4). These results suggest that miRNAs could regulate metabolism globally by preferentially regulating hubs, and control the downstream metabolic flux of the CPs with a local regulation strategy by preferentially regulating CPs. Hubs are the reactions that are shared by many pathways. miRNA regulation of hubs would affect many pathways and therefore such a regulation has a global effect on metabolism. For example, citrate synthase gene, encoding a major enzyme in citrate cycle (TCA cycle), could be regulated by a set of miRNAs: miR-152, miR-148a, miR-148b, miR-299-5p, miR-19b, miR-122a, miR-421, miR-494 and miR-19a. When these miRNAs regulate the citrate synthase gene, 78 pathways such as purine metabolism, pentose phosphate pathway, fatty acid biosynthesis, which represent lipid, carbon, nucleotide and amino acid metabolism, could be affected (Supplementary Data 5). This example illustrates that miRNA could regulate metabolism globally. On the other hand, GANAB, encoding glucosidase, a cut point in the network, could be regulated by miR-133a, miR-133b, miR-199a and miR-199b. Once GANAB is regulated by these miRNAs, the ouputs of three pathways, biosynthesis of steroids, pyrimidine metabolism and glycine, serine and threonine metabolism could be affected at the same time. In this way, a miRNA could regulate a few pathways in a coordinated manner.

Table 2

Fractions of miRNA Targets of Different Node Types

Taken together, miRNAs preferentially regulate network hubs and cut point nodes but avoid regulating the intermediate nodes. These results suggest that miRNAs have strategies to regulate metabolism by preferentially regulating certain types of nodes. In cellular signaling networks, we found that miRNAs avoid regulating the upstream nodes, i.e., ligands and receptors [5], while in metabolic networks, miRNAs avoid regulating intermediate nodes. In this study, we found that miRNAs preferentially regulate metabolic network hubs, which is in agreement with the previous studies that miRNAs preferentially regulate hubs in protein interaction networks and gene regulatory networks [6, 7]. It becomes a general principle that except in signaling networks, miRNAs preferentially regulate hubs of the cellular networks. In addition, miRNA preferentially regulates CPs for targeted tuning metabolic flows, while a similar local regulatory strategy of miRNA has been found in signaling networks, i.e., miRNA preferentially regulating the downstream components of adaptors in signaling networks [5].

Metabolic flows are arranged as linear or branching reactions. We first examined whether 2 and more consecutive linear reactions are enriched with miRNA targets. We found that 307 and 218 cases of the 2 and 3 consecutive linear reactions are all miRNA-targeted, respectively. Randomization tests showed that these cases are statistically significant (P<0.0002 and P<0.01, respectively, see Methodology). When performing the same analysis using the targets of individual miRNAs, no such statistically significant results were obtained. These results suggest that certain reaction regions (i.e., two or more consecutive linear reactions) are preferentially regulated by miRNA. For example, in N-Glycan biosynthetic pathway, fucosyltransferase 8 (FUT8), mannosyl (alpha-1,6-)-glyco-protein beta-1,2-N-acetylglucosaminyltransferase (MGAT2) and mannosidase (MAN2A1) are three enzymes that catalyze 3 consecutive reactions for N-Glycan biosynthesis. Their genes could be regulated by (miR-34c, miR-342, miR-449, miR-449b, miR-122a, miR-34a, miR-494 and miR-377), (miR-200b, miR-200c, miR-429, miR-181b, miR-181a, miR-495, miR-181c and miR-181d) and (miR-128b, miR-92b, miR-128a, miR-27a, miR-490, miR-32, miR-363, miR-500, miR-26b, miR-367, miR-25, miR-218, miR-92, miR-26a, miR-135b, miR-27b, miR-135a), respectively.

Branching metabolic flows could be convergent or divergent. To understand how miRNA regulate the branching metabolic flows, we extracted all the convergent or divergent units which are associated with three distinct reactions from the network. Convergent units could integrate metabolic flows while divergent units could allow a metabolic flow to diverge in two directions. We enumerated and extracted all the possible miRNA targeting patterns for the convergent or divergent units from the network and examined the statistical significance for each pattern. We found that both the convergent and divergent units without any miRNA targets are significantly enriched (Patterns 6 and 12 in Table 1, P<1.0 x 10^-4, randomization tests, see Methodology). In agreement with these observations, both the convergent and divergent units containing at least 2 miRNA-targeted nodes are significantly less-enriched (Patterns 1, 2, 5, 7, 8 and 11 in Table 1, randomization tests, see Methodology). Collectively, these results suggest that miRNAs preferentially regulate two or more consecutive linear metabolic reactions but avoid regulating metabolic branches.

To discover which metabolic pathways are regulated by miRNAs, we queried the KEGG human metabolic pathways using the reactions (nodes in the network), which are the miRNA targets, and conducted two types of analyses for determining the pathways that are significantly regulated by miRNA. We first performed an enrichment analysis for miRNA-targeted pathways using randomization tests. Using this method, we found that six pathways such as N-Glycan biosynthesis are enriched with miRNA targets (Supplementary Data 6, see Methodology). In the second type of analysis, we considered a pathway as miRNA-regulated pathway, if the pathway contains a miRNA-targeted reaction, which is a DSN or a CP that is the parent of the DSN in the network. Using this approach, we found that thirty pathways such as amino acid synthesis or degradation and other central metabolisms are extensively regulated by miRNA through targeting DSNs or CPs. Theses results suggest that miRNAs prefer to regulate DSNs or their parent-CPs. We merged the results from the two approaches, and found that in total 32 pathways are regulated by miRNAs (Table 3).

Table 3

miRNA Regulated Metabolic Pathways

As shown in Table 3, many biosynthetic pathways such as amino acid synthesis are extensively regulated by miRNAs. Glycan biosynthesis, pantothenate and CoA biosynthesis are also regulated by miRNAs. miRNAs also regulate certain lipid metabolism such as sphingolipid metabolism and glycerophospholipid metabolism. Similarly, miRNA regulation of central metabolic enzymes has been reported in Drosophila, i.e., miR-277 seems to be a regulatory switch for amino acid metabolism [20]. Amino acids are the building elements of proteins. Lipids are components of the plasma membrane or involved in cell signaling, while glycan and CoA are the major metabolic substrates for cell energy. These examples showed that miRNAs are predominantly involved in central metabolic pathways and thus play an important role in cell metabolism.