Regulatory Network Analysis of MicroRNAs and Genes in Neuroblastoma

Neuroblastoma (NB) is the fourth most common malignancy of childhood. And it is the most common intraabdominal malignancy of infancy and the most common extracranial solid tumor of childhood. Transcription factors (TFs) are one kind of proteins, which can bind to specific DNA sequences to control the transcription of genetic information from DNA to messenger RNA (Karin, 1990; Latchman, 1997). They have the ability to regulate gene expression at the transcription level alone or with other proteins. MicroRNAs (miRNAs) are small (about 22 nucleotides) non-coding RNAs, which are involved in various biological processes including cell proliferation, differentiation and apoptosis. Therefore, they can regulate the expressions of genes at the posttranscriptional level and become important cellular components (Chen et al., 2007). The genes, where miRNAs are located inside, are named as host genes. Rodriguez et al. (2004) used to


Regulatory Network Analysis of MicroRNAs and Genes in Neuroblastoma
Li Wang, Xiang-Jiu Che*, Ning Wang, Jie Li, Ming-Hui Zhu suggest that miRNAs can be transcribed parallelly with their host transcripts, and the two different transcription classes of miRNAs, exonic miRNAs and intronic miRNAs, identified in their study may require slightly different mechanisms of biogenesis. Intronic miRNAs are closely associated with their host genes. And Baskerville et al. (2005) identified that they were usually equally expressed in biological progression, acting as a potential partner in order to achieve biological function and affect the alteration of pathways (Cao et al., 2010). According to the above results, it suggested that miRNAs and their host genes together or separately were capable of contributing to cancer progression. In this study, we considered host genes were also differentially-expressed genes if their miRNAs involved in cancer progression.
The miRNAs, their host genes and TFs involved in the study are significant regulators for gene expression. And it's critical for the regulation of miRNAs by TFs. More and more evidences indicate that aberrant regulation of miRNAs by TFs can cause phenotypic variation and diseases. Our study's experimental data shows that differentially-expressed genes and miRNAs play key roles in NB's development, metastasis and therapy. MYCN, for example, which is a differentially-expressed gene targeting eight differentially-expressed miRNAs, may be helpful in predicting whether NB cells possibly proliferate and spread to nearby tissues.
Although numerous NB experiments have identified many differentially-expressed genes and miRNAs, none of them tried to find out the great picture about the regulatory network of them, which has an entire connection among them and perhaps help people further identify the key pathways controlling mechanisms of NB. Our study's main aim was to determine the associations of genes, miRNAs, host genes and TFs involved in NB. Therefore, these extracted transcriptional relations were regarded as a penetration point to build the regulatory network of the genes and miRNAs in NB. There were three networks used to be considered: differentially-expressed network, related network and global network. However, the global network was much more complex and difficult to find clear pathways associated with NB, we decided to use the other two networks for further study. NB formation was partially revealed by the differentially expressed network, and we identified the topology network in NB progression. The work of comparing and analyzing similarities and differences among the three networks was great useful to identify the key nodes and pathways. We will study further into the mechanism of NB with increasingly comprehensive data in the future, and hope this study will be helpful to pertinent research into the pathogenesis and treatment of NB.

Material collection and data processing
The two experimentally validated data of miRNA targets and miRNA-target interactions in this study were extracted from TarBase 5.0 (Sethupathy et al., 2006) and miRTarBase (Hsu et al., 2011) separately. As a comprehensive dataset of experimentally supported animal miRNA target, TarBase 5.0 can provide us large validated data of miRNA targets. And it is great help to build up the interactions network with miRNA-target interactions extracted from miRTarBase. As we know, there are various different symbols representing the same miRNAs and genes, so it's necessary to systemize the symbolic representation method. Official symbols from the National Center for Biotechnology Information (NCBI) database were used. The complete data were considered as set D1.
The interactions of human TFs and miRNAs were extracted from TransmiR (Wang et al., 2010). TransmiR, which is free for academic usage, is a database for transcription factor-microRNA regulation. The complete data were considered as set D2.
NCBI and miRBase (Griffiths-Jones et al., 2008) supported this study strongly with providing experimentally data of host genes of miRNAs and confirmed human miRNAs separately. The complete data were considered as set D3.
The whole differentially-expressed genes were partly extracted from Cancer Genetic Web and partly manually obtained from relevant literature on human neuroblastoma. Similarly, the related mutated genes of NB were manually obtained from SCI studies as well. In addition, we used the P-match method (Chekmenec et al., 2005) to extract some popular TFs. These TFs were also considered as related genes. As a new tool, P-match was designed for identifying transcription factor (TF) binding sites in DNA sequences by combining pattern matching and weight matrix approaches. It can provide higher accuracy of recognition than any other methods alone. In this study, we used the P-match method to identify the TF-binding sites (TFBSs) mapped onto the promoter region sequences, which were in the 1, 000-nt promoter region sequences targeted by differentially-expressed miRNAs and downloaded from UCSC database. Using the matrix library and sets of aligned known TFBSs collected in TRANSFAC, the P-match is able to provide the possibility to search for a large variety of different TFBSs. In order to indicate the matrix with a cut-off value allowing a false-negative rate of 50%, while requiring the false-positive rate to reduce below a certain threshold, we decided to use the vertebrate matrix and restricted high-quality criterion. Therefore, the matrices defined as low-quality matrices can produce the highest number of false-positive matches. The complete differentially-expressed genes and related genes were considered as set D4.
Differentially-expressed miRNAs of NB extracted from mir2Disease (Jiang et al., 2009) and extra differentiallyexpressed miRNAs manually obtained from recently related SCI studies make sure to provide a comprehensive resource for the study. The complete data were considered as set D5.

The construction of three networks
Differentially expressed genes and miRNAs are primary elements for building up these very complex transcriptional networks of NB because their differential expressions can be found in every process of cancer formation. Thus, the construction of the core regulation network was utilized by the following method: differentially expressed genes and miRNAs from set D4 and D5 were mapped onto set D1, D2 and D3, then the regulatory network of miRNAtargets, TF-miRNA and hostgene-miRNA were used to determine the core network.
As well as differentially expressed genes and miRNAs, many other NB-related genes and miRNAs actually are also involved in each progress of cancer formation. Thus, in order to demonstrate the regulatory network further, we combined the core network and more complex regulatory network as the NB-related network using the same method mentioned above.
In the third network, the pathways of TFs and miRNAs from related network may not be as crucial as the former two networks, but some certain genes and miRNAs are just not proved by validated experiments to be involved in the progression of NB. We mapped them onto set D1, D2 and D3 to extract the regulatory relations of TF-miRNA, miRNA-targets and hostgene-miRNA. After combining all these relations, we got the expanded global network.

Core regulation network of NB
Using the methods and datasets mentioned above, we succeed to construct a core transcriptional network that could show clearly the regulation mechanism of NB. There are three TFs, i.e., MYCN, PDGFA, TP73, and another 11 differentially expressed genes such as BCL2 and VEGFA, and also 37 differentially expressed miRNAs and their host genes in this network, which are regarded as essential regulatory elements. Figure 1 shows a few significant regulation pathways, MYCN, for example, is targeted by hsa-miR-34a and hsa-miR-29a, which regulates hsa-miR-106a, hsa-miR-92a, hsa-miR-335, hsa-miR-19b, hsa-miR-221, hsa-miR-18a, hsa-miR-19a, and hsa-miR-17. And hsa-miR-106a and hsa-miR-17 target the genes VEGFA, BCL2 and RBL2. Thus, we can see that hsa-miR-29a, by targeting MYCN, which regulates hsa-miR-17, can regulate hsa-miR-17 and affect on the expression of gene RBL2 at the end of the pathway. Besides, there are two different pathways that hsa-miR-29a can impact on the expression of gene VEGFA. One is that hsa-miR-29a targets MYCN, MYCN regulates hsa-miR-106a, and hsa-miR-106a targets gene VEGFA; The other is that hsa-miR-29a targets MYCN, MYCN regulates hsa-miR-17, and hsa-miR-17 targets gene VEGFA. As above, hsa-miR-34a also can use the same pathway to affect on gene RBL2. As a member of the MYC proto-oncogene family, MYCN is a transcription factor that controls expression of many target genes, which regulate fundamental cellular processes including proliferation, cell growth, protein synthesis, metabolism, apoptosis and differentiation (Westermark et al., 2011). Gene VEGFA (VEGF) is a specific endothelial cell mitogen that stimulates angiogenesis and plays a crucial role in tumor growth, and Langer et al. (2000) revealed that VEGFA affected neuroblastoma cell growth directly and could be an autocrine growth factor. van Golen et al. (2000) identified that BCL2 overexpression in neuroblastoma cells promotes cell survival by preventing mitochondrial membrane depolarization and caspase-3 activation, ultimately leading to increased tumor growth. In this network, we see the connection between MYCN and VEGFA, BCL2, and understand that gene MYCN can affect on the expression of gene VEGFA and BCL2 by regulating hsa-miR-17 in order to regulate cellular processes, such as proliferation, cell growth and protein synthesis. We also noticed that two TFs (TP73 and PDGFA) regulate hsa-miR-145 and hsa-miR-221 separately without any miRNAs targeting them. As we know, genes ALK and PHOX2B are also considered as differentially expressed genes that are involved in the pathogenesis of several different human cancers including NB (Di et al., 2011). And ALK tumors were identified to represent a specific subset of NB and usually develop in infants (Ogura et al., 2012). However, we did not find them having any connection with differentially expressed miRNAs of NB in this core network. This suggested that genes like ALK and PHOX2B might have other specific pathways without differentially expressed miRNAs already identified, which we have not known yet. Several pathways found in the present study partially revealed the regulatory mechanism of NB.

Related network of NB
The NB related network consists of differentially expressed genes and miRNAs, related genes and miRNAs, targets of miRNAs and host genes of miRNAs. Thus, it included the core network and thus more complex to analyze. Figure 1 shows three TFs, i.e., MYCN, FDGFA and TP73 in the core network. There are 14 TFs in Figure  2, which include the three essential TFs (Figure 1) and regulate 53 mutated miRNAs. MYC, MYCN and TP53, among the TFs in this network, regulate more miRNA expression, therefore they will be more potentially influential than other TFs. Figure 2 shows many additional pathways that are possible to affect the progression of NB. Related gene TP53 regulates hsa-miR-107, hsa-miR-125b, hsa-miR-143, hsa-miR-145, hsa-miR-192, hsa-miR-200b, hsa-miR-29a, hsa-miR-29c, hsa-miR-34a, hsa-miR-34c. And these 11 miRNAs target 39 relevant genes including VEGFA, BCL2, MYCN and TP53. Naves et al. (2013) revealed that mutated gene TP53 in neuroblastoma cells could be linked with the switch between apoptotic response and cell death by autophagy in response to hypoxia mimetic stress, and it can be partially revealed from our network. MYC, another important TF, regulates 17 relevant miRNAs, such as hsa-let-7b, hsa-miR-106a, hsa-miR-34a, hsa-miR-92a and hsa-miR-17. These miRNAs are all involved in the core regulatory network and therefore related gene MYC should have a high potential of being more influential on the progression of NB. We can also see that there are several feedback loops, i.e. PTEN and hsa-miR-19a, PTEN and hsa-miR-21, MYC and hsa-miR-34a, MYC and hsa-miR-17, TP53 and hsa-miR-125b, AKT1 and hsa-miR-125b, and IFNB1 and hsa-miR-145.

Figure 1. Core Regulation Network of Differentially Expressed Genes and miRNAs in Neuroblastoma
The associated network not only expands additional topology associations of differentially expressed elements but also makes it more clear and helpful to understand the mechanisms of NB.

Global network of NB
The global network has additional TFs, targets, miRNAs and host genes of miRNAs compared to the related network. It shows more comprehensive regulatory relations including all relations in D1, D2 and D3. Therefore it includes both the core regulatory network and the related network.

Host genes and miRNAs in NB related network
Although differentially expressed genes don't include miRNAs' host genes, these host genes and their differentially expressed miRNAs can still affect on the progression of NB. Figure 3 presents several pathways of host genes and miRNAs. For example, MIR34A includes hsa-miR-34a and hsa-miR-34a targets BCL2, CD44, MAGEA3, MAP2K1, MYCN and VEGFA. Several miRNAs can be located in one host gene, for example, hsa-miR-17, hsa-miR-19a and hsa-miR-18a were all located in host gene MIR17HG. Meanwhile several host genes  This result suggested that investing the associations of host genes and their miRNAs, especially when their miRNAs were differentially expressed, was helpful to understand the pathogenesis of NB.

Regulatory pathway of differentially expressed genes
In this study, the upstream and downstream information of important elements: differentially-expressed genes, differentially-expressed miRNAs and TFs were extracted to describe the NB network more clearly.
In order to compare and analyze the regulatory pathways of differentially expressed genes, we extracted the successor and precursor nodes of the differentially expressed genes from the three-level networks. Differentially-expressed gene MYCN, for example, we built Table 1 to show its predecessors and successors in the three networks and help us understand its regulatory associations. As a notable tumor suppressor, MYCN has significant features in these networks. Table 1 shows that two miRNAs (hsa-miR-29a, hsa-miR-34a) target gene MYCN and MYCN regulates another eight miRNAs in the core network. Thus, we realized that these two miRNAs indirectly had an effect on the expression of the eight miRNAs via MYCN. In Figure 1, TFs like PDGFA and TP73 only have successors in three networks. In addition, other genes like BCL2, VEGFA and IGF1R, which were not TFs, were investigated. These genes only were targeted by specific miRNAs and did not regulate any miRNA.

Regulatory pathway of differentially expressed miRNAs
Using the same method as the differentially expressed genes, we compared and analyzed the regulatory pathways of differentially expressed miRNAs. Taking hsa-miR-17 for example, we listed all its predecessors and successors in the three networks. Table 2 presents that MYCN regulates hsa-miR-17 to affect on the expression of genes BCL2, RBL2 and VEGFA. Besides, MYC and hsa-miR-17 in related network, CCND1 and hsa-miR-17, E2F1 and hsa-miR-17 in global network separately form self-adaptation relation. Hsa-miR-17 also indirectly affects other miRNAs via certain TFs. For example, hsa-miR-17 targets PTEN, which regulates hsa-miR-19a and hsa-miR-21.

Regulatory pathway of popular TFs
We totally found 14 TFs and noticed that there were several classes of TFs involved in the related network. There are six types of adjacent nodes in the first class of TF. And this class of TF has three types of predecessors and three types of successors. We have identified that MYCN, MYC, AKT1, IFNB1, PTEN, FOXO3 and TP53 all belong to this class of TF. Taking MYC as an example, Table 3 shows the relationship of MYC and miRNAs in all the three networks. There are four miRNAs (hsa-miR-145, hsa-miR-17, hsa-miR-21, and hsa-miR-34a) and 14 miRNAs separately being as its predecessors and successors in the core network. There are five miRNAs (hsa-miR-145, hsa-miR-17, hsa-miR-21, hsa-miR-34a and hsa-miR-34c) targeting MYC and 16 miRNAs being regulated by MYC in the related network. In the global network of NB, 24 miRNAs are identified to target MYC, and MYC regulates 57 miRNAs. We also found that hsa-miR-17 and MYC formed self-adapting associations in all the three networks, as well as hsa-miR-34a and MYC. Meanwhile, hsa-miR-17 and hsa-miR-34a are all differentially-expressed miRNAs in NB. According to this result, we knew that these upstream miRNAs indirectly affected the downstream miRNAs through TFs, which means TFs can be regarded as the bridge of connecting predecessors and successors, and affect on the expressions of miRNAs in the key pathways if they have changed.
There are five types of adjacent nodes in the second class of TF. And this class of TF has three types of predecessors and two types of successors. For example, CDKN1A has six miRNAs (hsa-miR-93, hsa-miR-17, hsa-miR-145, hsa-miR-125, hsa-miR-10b and hsa-miR-106a) as predecessors and no successors in the core network.
There are three types of adjacent nodes in the fouth class of TF. And we identified two kinds of the fouth class in the related network. One kind has one type of predecessors and two types of successors. For example, IL6 is only targeted by eight miRNAs in the global network, and regulates hsa-miR-152 in the related network and six miRNAs in the global network. The other kind has no predecessors and three types of successors. SRC, TNF, and IFNG all belong to this kind of class, which means no miRNAs target them but they regulate miRNAs in all the three networks.

Discussion
We also investigated some genes that had only three types of adjacent nodes (three types of predecessors). For example, BCL2 was only targeted by some miRNAs and did not regulate any miRNAs in all the three networks, which meant that genes like BCL2 might be the last nodes of the pathways.
There is a phenomenon needed to be changed that researchers used to spend much time only on one single cancer instead of combining different cancers with same types of miRNAs and genes, which could not be helpful to see the clearly networks between cancers. Actually, before our investigation on the networks of differentexpressed miRNAs, different-expressed genes and TFs, we had an idea that whether there existed a related networks between various human cancers, such as Breast Cancer and Prostate Cancer, in the level of miRNAs and genes. Fortunately, we found evidences. Liu et al. (2012) illustrated that hsa-miR-17 and hsa-miR-106a, which are significant different-expressed miRNAs in our research, were also differently expressed in Prostate Cancer (PC) and had profound influence on PC. In our both researches, the regulation of hsa-miR-17 targeting VEGFA gene is an important part of the key pathway in both NB and PC. And Tong et al. (2012) revealed that mutated gene AKT1, which also has a high potential of being influence on the mechanisms of NB, with other mutated genes play a key role on the key EGFR pathway in human Breast Cancer (BC). Furthermore, we found that this EGFR pathway might have a related network with MYC, which is as a related gene in NB and different-expressed gene in BC (Zhang et al., 2013). These indicated that the same different-expressed miRNAs or genes do have possibilities to influence on different human cancers, which also means it is possible that we can extract several miRNAs and genes having this feature and take some certain measurements to stop them differently expressing in order to helpfully prevent cancers if we sufficiently figure out the inter-relationship of different cancers and combine valid results together. We hope this method might be useful to prevent human cancers in the future. To achieve this goal, we will focus closely on similar researches and make most use of comprehensive information to analyze mechanisms of cancers in the future.