Analyzing the miRNA-Gene Networks to Mine the Important miRNAs under Skin of Human and Mouse

Genetic networks provide new mechanistic insights into the diversity of species morphology. In this study, we have integrated the MGI, GEO, and miRNA database to analyze the genetic regulatory networks under morphology difference of integument of humans and mice. We found that the gene expression network in the skin is highly divergent between human and mouse. The GO term of secretion was highly enriched, and this category was specific in human compared to mouse. These secretion genes might be involved in eccrine system evolution in human. In addition, total 62,637 miRNA binding target sites were predicted in human integument genes (IGs), while 26,280 miRNA binding target sites were predicted in mouse IGs. The interactions between miRNAs and IGs in human are more complex than those in mouse. Furthermore, hsa-miR-548, mmu-miR-466, and mmu-miR-467 have an enormous number of targets on IGs, which both have the role of inhibition of host immunity response. The pattern of distribution on the chromosome of these three miRNAs families is very different. The interaction of miRNA/IGs has added the new dimension in traditional gene regulation networks of skin. Our results are generating new insights into the gene networks basis of skin difference between human and mouse.


Introduction
The integument includes the skin, coat/hair, and nails [1]. The mammalian coat hair is one of the defining characteristics of mammals [2]. Enormous morphological variations are found among mammalian integument, especially in coat hair. Hair is present in differing degrees in all mammals [3]. In most mammals, the hair is abundant enough to cover the body and form a thick coat, while dolphins, naked mole-rats, and humans are among the most hairless of all mammals [4]. However, humans with the hypertrichosis syndrome have hair covering their faces, their eyelids, and their bodies [5,6]. Therefore, we provide a hypothesis that key genes related to hair coat formation exist in every mammalian species. Mammals have diverse coat morphology, due to changes in the gene regulatory pathway/network. Common genetic network variants have been shown to affect many complex traits, including integument morphology. Therefore, to understand the differences of molecular mechanism for integument phenotype between mouse and human, it is necessary to employ the scope of system biology. The comparison of gene coexpression networks and regulation networks between two species is often a useful method to identify the differences of critical biologic process associated with morphology diversity.
Mammalian Phenotype (MP) ontology [7] (http://www .informatics.jax.org/) is made up of 17,330 terms (as of 02/08/2016), and most of these terms were characterized 2 BioMed Research International from abnormal mouse phenotypes. The ontology database deposited 1627 mouse/human orthologous genes with phenotype annotation related with integument (MP:0010771). High-throughput experiments have produced many microarrays and next generation sequencing data that are collected in Gene Expression Omnibus (GEO) database [8]. Genetic networks have been widely used in biological research, bridging the gap between single genes and biological systems by investigating the relationships between different genes. miRNAs are small noncoding RNA molecules which play an essential role during skin development [9]. Variations in the interaction between miRNA and target gene are likely to influence the phenotypic differences between species. Here, we use Weighted Gene Co-Expression Network Analysis (WGCNA) to reveal shared and unique characters of the gene networks in human and mouse skin. We identify networks of coexpressed genes, which might be associated with biological functionally relevant coat morphology, and explore differences in these biological processes between two species. We also found that candidate miRNA may play a role in the regulation of gene expression networks in skin. These results provide a system-level insight into evolutionary changes of the integument.

Materials and Methods
In the present study, we emphasize the comparison between the gene network for integument of human and that of mouse and identification of candidate miRNAs relevant to the integument genetic pathway and constructed an interaction network between miRNAs and targeted genes. Figure S1 in Supplementary Material available online at http://dx.doi.org/ 10.1155/2016/5469371 shows pipeline for the study.

Selection of Genes and miRNAs for Human and Mouse.
Mouse/Human Orthology with Phenotype Annotations were downloaded from the MGI database [10]. 1627 of these genes are related to the integument phenotype (MP:0010771). To draw the chromosome location for integument genes (IGs), the R package org.Mm.eg.db [11] and org.Hs.eg.db [12] were used. We downloaded the gene expression data of skin for mouse and human from the Gene Expression Omnibus (GEO) [13] and mouse and human miRNAs from miRBASE (http://www.mirbase.org/cgi-bin/browse.pl/) [14]. Mouse and human 3 UTR sequences of 1627 integument genes were obtained from Ensemble (http://www.ensembl.org/) by biomaRt in R language [15].

Construction of Gene Networks.
To compare the gene networks of skin for human and mouse, we selected many relevant data from the GEO database [8], and then we filtered out all but a core collection of datasets that were similar enough for useful bioinformatic comparison. First, we downloaded all datasets that were run on same Affymetrix platform, one platform in human (GPL570) and one in mouse (GPL1261) (Table S1). Second, we selected only relevant experiments about skin samples on the platform. Third, we extracted IGs expression data from the chip matrix. Finally, a total of 1487 genes were analyzed in this study. Then, Weighted Gene Co-Expression Network Analysis (WGCNA) [16,17] was used for these expression datasets to create consensus networks between human and mouse, and the networks were visualized in Cytoscape 3.3.0 [18].

Functional Enrichments for IGs Network.
To classify the terms and group for IGs in gene ontology (GO) and KEGG pathway, we used the Cytoscape plug-in ClueGO to implement enrichment analysis. We set kappa score threshold to 0.5 and the value to 0.05; we used two-sided test, Bonferroni step down and GO term fusion. The other parameters of the software were set to default values [19].

Prediction of miRNA Targets for Integument Genes.
For prediction of miRNA targets for hub genes, miRanda [20] software version aug2010 available at http://cbio.mskcc.org/microrna data/miRanda-aug2010.tar.gz was employed to predict miRNAs regulated IGs. miRanda was running as the following command options: sc ≥ 180, en = 1 (no energy filtering), go = −9.0, ge = −4.0, and scale = 4.0. We set the option -strict to prevent gaps and noncanonical base pairing in target sites. Human/mouse miRNA sequences were used as query sequences, and the IGs sequences were used as [20]. We also predicted the targets with two miRNA prediction databases, miRTarBase [21] and miRDB [22], and selected the intersections targets to construct the miRNA-mRNA regulatory network.

Chromosome Location for Integument Genes of Human and Mouse.
We observed a uniform distribution of the 1627 integument genes across all chromosomes in human and mouse (Supplementary Figure S2). As can be seen in Figure  S1, in general, 1627 integument genes are dispersed along the chromosomes of human and mouse. These results provide insights into the biological underpinnings of integument phenotype and the pleiotropic connections between traits. It is very difficult to detect which chromosome is essential to integument phenotype because these integument genes will be uniformly distributed on all chromosomes.

Gene Expression
Network. By using WGCNA, we generated two consensus networks of the human and mouse (Figure 1) which show a network heat map plot of a gene network together with the corresponding hierarchical clustering dendrograms and the resulting modules. Then, the gene expression networks were visualized in Cytoscape. There are two groups in human skin consensus networks constructed from 560 genes (nodes) and 18391 interactions (edges), while only one group in mouse skin gene networks was constructed from 368 genes and 1757 interactions ( Figure 2). To evaluate the complexity of gene networks, we calculated the interaction degrees for gene expression network. The average interaction degrees are 32.8 (18391/560) and 4.8 (1757/368) for human and mouse, respectively.
We use the ClueGO to compare three clusters of genes from Group 1 and Group 2 of human and the mouse   consensus networks (Figure 2). 44 significantly overrepresented categories were identified ( Figure 3). Most of the positive regulations of cellular process were clustered into two independent modules. More than half of the common genes among three clusters were involved in regulation of cellular and metabolic processes, such as negative regulation of nucleobase-containing compound metabolic process and regulation of macromolecule biosynthetic process. Some of the most significant categories were response to stimulus, regulation of localization, negative regulation of biological process, cell differentiation, single-multicellular organism process, and response to organic substance. Analysis of the modules showed that the majority of the response to stimulus genes and response to organic substance were common in three clusters.
Furthermore, we found that four GO terms were highly enriched, and these terms were specific in human compared to mouse. These four GO terms could be divided into three categories ( Figure 4). The regulation of secretion category contains two GO terms: regulation of secretion and regulation of secretion by cell.     main miRNA in the total miRNA target file, we selected the top 10 miRNA families based on the target number of IGs ( Table 1). The hsa-miR-548 family has predicted 4643 target sites distributed on 541 IGs of the human. The mmu-miR-466 family and mmu-miR-467 family have predicted 1,704 target sites and 956 target sites, respectively, distributed on 375 and 310 IGs of the mouse. The hsa-miR-548 family is widely distributed in the whole human genome. However, both of these two miRNA families are located in an intron of Sfmbt2 gene on chromosome 2 of the mouse ( Figure 5). These three miRNA families have common target genes ( Figure 6) and might present some similarity function to skin morphogenesis.

Construction of a miRNA-mRNA Regulatory Network.
Two miRNA prediction databases miRTarBase and miRDB were used to verify the results of miRanda. The intersections targets of three algorithms were selected to construct the miRNA-mRNA regulatory network (Table S2). The resulting miRNA/target mRNA pairs and gene networks for human and mouse were visualized in Cytoscape, where edges from miRNA to genes represent a potential regulatory relationship and edges from gene to gene represent an expression correlation. The network consisted of two almost separate groups, which could be connected when miRNAs is added to the expression gene network of human (Figure 7). These miRNAtarget interactions were added to the network of the mouse,

Discussion
Recently, a lot of biology databases have been published. How to integrate these resources in system biology is a big challenge for analysis of certain biological problems [23]. In this study, we have integrated the MGI, GEO, and miRNA databases to analyze the genetic regulatory networks under morphology difference of integument of human and mouse. 1627 mouse/human orthologous genes related with integument phenotype were deposited in MGI database. These genes were widely scattered across the mouse genome and the human genome, which suggest that the integument is a quantity phenotype with polygenic determinism. We also have constructed the expression correlation networks of IGs with the gene expression matrix from the GEO database.
With the process of evolution, the organisms have increased in complexity. However, the organism's complexity is not simply associated with the number of its genes [24]. Now, the interaction degrees as a straightforward detector were used to evaluate the complexity of gene networks [25]. Our results revealed that the networks of human skin are much more complex than those of mouse which might be explained by the fact that the integument structure at the anatomical level is much more complex in human, compared with mouse [26]. Many researches reported that one gene could evolve new function to adapt to the change of environment during species divergence, from lower to higher species [27,28]. Erwin and Davidson reported that the reorganizations of gene networks could change the animal morphology and the basic of the evolution is regulatory changes within a gene regulation network [29]. These results also confirm our hypothesis that those key genes related with the hair coat exist in every mammalian genome and the diverse morphology in mammal just because of the difference in gene regulation networks.
Analysis of the GO terms category showed that the majority of the positive regulations of cellular process and response to stimulus genes were common in three clusters, which may be due to the perception function of skin to in vivo and in vitro environments [30,31]. In other words, these common categories in human and mouse skin were involved in protection [32,33], sensation [34], temperature regulation [34], immunity [35], exocrine, and endocrine. By comparing the GO terms for those genes in expression network, the term of regulation of secretion category was enriched in human skin. Human has eccrine glands and thermal apocrine glands. Apocrine glands are always associated with hair follicles, and eccrine sweat glands cover almost the entire body surface of human [36], while eccrine sweat glands in the mouse are found only on the footpads [37]. These results might provide an important insight into the evolution of thermal eccrine system in human. The secretion category of the genes might be involved in the eccrine system in human.
miRNA plays an essential role in the regulation of skin development and morphogenesis [38]. miRNA family is a group of miRNAs that share common seed sequences, which has a similar regulation function [39]. We found that the hsa-miR-548 family has the highest amount of target sites among the identified miRNAs in human and the mmu-miR-466 and mmu-miR-467 families are top two in the miRNAs list predicted in the mouse. Generally, increasing the number of target sites within a gene improves the efficiency of miRNA regulation [40]. More target sites provide more opportunity for recognition by miRNA and increase the kinetics and overall level of transcript regulation [41]. miR-548 is a primatespecific miRNA family, which has 69 members distributed in almost all human chromosomes [42]. The hsa-miR-548 family is involved in multiple biological processes, such as signaling pathways, immunity, and osteogenic differentiation, and some cancers [43][44][45][46][47]. hsa-miR-548 takes part in IFN signaling which responds to the virus and bacterial infections on the cell [46]. hsa-miR-548 also can turn down the host antiviral response by degradation of IFN-1 [43]. UVB irradiations can downregulate hsa-miR-548 of human epidermal melanocytes cell [48]. These results suggested that hsa-miR-548 might contribute to dynamic regulatory network of skin transcriptome of human. Compared to human, mmu-miR-466 and 467 families only located in intron 10 of Sfmbt2 genes. This area of intron 10 has a larger cluster of miRNAs, which is specifically present in mouse and rat [49]. Based on miRBase database (http://www.mirbase.org/) definition, clustered miRNAs are a group of miRNAs located within 10 Kb of distance on the same chromosome. In this study, mmu-miR-466, mmu-miR-467, and mmu-miR-669 clusters have one core promoter region and transcriptional start site shared with the Sfmbt2 gene. The expressions of mmu-miR-466 and mmu-miR-467 markedly waved during the hair follicle cycling in mouse [50] and were downregulated in melanoma of mouse by curcumin diet [51]. And histone deacetylation and metabolic oxidative stress can induce the activity of mmu-miR-466 [52]. Skin serves as a barrier between body and pathogens (disease causing organisms), which is part of the innate immune system [53]. Both of hsa- miRNA-548 and mmu-miR-466 and mmu-miR- TNFSF11 TPH2  COL19A1 KIF1A  BTN1A1 KSR1 MC5R  HOXB1  SLC6A3 AGER  SSTR4  ATOH1 LALBA  CAMK2A  P2RX3 RAC3  HCN1  HOXB13 OXT  APOA1  GRM1 SOX2  NKX2-5 CSN3  RXFP2  TRPM3  POU3F4  PRKCA CAMK4 PITX3  ZAP70 ASCL1  GATA4 STOML3  CDK5R2 PLD4  WNT7A  HAP1  POMC NOS2 FSHR  LMX1A CDH8  CLEC1B  PROKR2 POU3F2 CASR HTR2C  CD40LG  BRSK1  USH1C GHRHR  PRDM8 DRD2  GFI1B  SYT7 PTGER1  FOXB1  SLC5A7 CACNA1S  NTSR2   CYP19A1   ADORA1   SLC38A10  ATE1   OXTR  SEMA6C CSF2  HOXD1  VGF  GLRA3  LTC4S  ACADS   RARG   PORCN  ANO1   NR1H2  POLD1  HOMER3  ZFPM1   TOM1L2 RXRB  FGF6 SLC39A3  IL17RA  VHL  TYRO3  GABRA2 SLC12A1   E2F5   INPP5D   LATS1  ARX  TCF15  GNAO1   SOX11 MYOG  HOXD13 OPRD1 ALX4  P2RX6 FGF20 NAGS  OPRK1  RLBP1 GATA1  MTA2 NMUR1  CLRN1 PTCH2  TREX2 CDH23 FGA  AQP2 SRC  CCR10 CTSE  HFE2 PDX1  KCNIP3  TLL1  PTGIR  OPRM1  CHRNA4 NPB  SNAP25 CHST8  HTR1F  CITED1  EMX1 CNR2  TCF21  KCNJ11 ANK1 CLCN1  CLPS  SLC17A6  ADAMTS20 RXFP1  TMPRSS6 GBX1 AICDA  FOXP3  AMHR2 NTSR1   SLC6A4  CLDN6 FGF5  HTR3A  GPRC6A  GRIK1 FGF23 KISS1R KRT76  MKL1  SCN11A  NTRK1  RSPO2 RLN1  GALR2 AVPR2  GRIA1 IL4  GSC LMX1B  GPR3 TMC1 hsa-miR-3927-3p hsa-miR-3919  inhibit the host immunity response [54]. It is necessary to carry out the research on how these miRNAs contribute to integument morphogenesis in the future. When we added the relationship of miRNAs and their target genes by three prediction algorithms to the expression correlation networks, those two clusters in gene regulatory networks of human have been integrated, which add a new dimension to genetic networks under human integument. The changes of the hub gene in gene regulatory networks result in the evolutionary alternation and the morphological difference [29].

Conclusions
Genetic networks variants have been shown to affect many complex traits, including integument morphology. In this study, we try to compare the regulatory networks of miRNAs and IGs in the skin of human and mouse. We downloaded mRNA expression data in human and mouse skin from the GEO database to create within-species consensus networks by WGCNA. There is a big difference in consensus networks between human and mouse; human consensus networks are more complex compared to mouse. Two principal regulatory networks were found in human: one module contains 286 IGs specifically involved in secretion, whereas the other module contains 250 IGs which are enriched for cellular response to stress and catabolic process. The secretion category, which is specific category for human compared to mouse, of the genes might be involved in eccrine system in human. Then, total 62,637 miRNA binding target sites were predicted in human IGs, while 26,280 miRNA binding target sites were predicted in mouse IGs. The interactions between miRNAs and IGs of human are also more complex compared to mouse. To further detect the role of these miRNAs, miRNA/IGs specific regulatory networks were added in IGs expression correlated networks, which will advance the dimension to skin regulation networks. hsa-miR-548, mmu-miR-466, and mmu-miR-467 have an enormous number of targets on IGs, which both have the role of inhibition of host immunity response. The regulations of transcription factors to downstream genes and miRNA to transcription factors affect the spatial and temporal expression of genes during skin morphogenesis. Our results provide a new avenue to understand the genetic networks basis of skin difference between human and mouse.