key: cord-0005522-se7vbrdi authors: Pan, Wei-Kang; Zhang, Ya-Fei; Yu, Hui; Gao, Ya; Zheng, Bai-Jun; Li, Peng; Xie, Chong; Ge, Xin title: Identifying key genes associated with Hirschsprung’s disease based on bioinformatics analysis of RNA-sequencing data date: 2017-01-25 journal: World J Pediatr DOI: 10.1007/s12519-017-0002-0 sha: 364cead576f72be17ef8bf8da931dfdca2aed5c6 doc_id: 5522 cord_uid: se7vbrdi BACKGROUND: Hirschsprung’s disease (HSCR) is a type of megacolon induced by deficiency or dysfunction of ganglion cells in the distal intestine and is associated with developmental disorders of the enteric nervous system. To explore the mechanisms of HSCR, we analyzed the RNA-sequencing data of the expansion and the narrow segments of colon tissues separated from children with HSCR. METHODS: RNA-sequencing of the expansion segments and the narrow segments of colon tissues isolated from children with HSCR was performed. After differentially expressed genes (DEGs) were identified using the edgeR package in R, functional and pathway enrichment analyses of DEGs were carried out using DAVID software. To further screen the key genes, protein-protein interaction (PPI) network and module analyses were conducted separately using Cytoscape software. RESULTS: A total of 117 DEGs were identified in the expansion segment samples, including 47 up-regulated and 70 down-regulated genes. Functional enrichment analysis suggested that FOS and DUSP1 were implicated in response to endogenous stimulus. In the PPI network analysis, FOS (degree=20), EGR1 (degree=16), ATF3 (degree=9), NOS1 (degree=8), CCL5 (degree=8), DUSP1 (degree=7), CXCL3 (degree=6), VIP (degree=6), FOSB (degree=5), and NOS2 (degree=4) had higher degrees, which could interact with other genes. In addition, two significant modules (module 1 and module 2) were identified from the PPI network. CONCLUSIONS: Several genes (including FOS, EGR1, ATF3, NOS1, CCL5, DUSP1, CXCL3, VIP, FOSB, and NOS2) might be involved in the development of HSCR through their effect on the nervous system. H irschsprung's disease (HSCR), also called Hirschsprung disease, congenital aganglionic megacolon or congenital megacolon, is a type of megacolon induced by deficiency or dysfunction of ganglion cells in the distal intestine. [1] In HSCR, the affected colon segment cannot relax the bowels; therefore, HSCR is characterized by fecal retention, intestinal obstruction, constipation and abdominal distention. [2] HSCR usually occurs in children, especially boys, and seriously impacts quality of life and, possibly, mortality. [3] In Japan, the incidence rate of HSCR is approximately one in 5000 births. [4] Treatment strategies of HSCR include resection of the abnormal colon segment and reanastomosis. However, some children still have bowel dysfunction after treatment and need a lifetime of intestinal fistulas. [5] Thus, it is necessary to investigate the mechanisms of HSCR and to develop novel targeted treatment strategies for this disease. To date, several studies have explored the key genes involved in HSCR. For instance, the RET gene, which encodes a tyrosine-kinase receptor, plays an important role in the development of HSCR by affecting the enteric ganglia. [6] Previous studies have demonstrated that neuregulin 1 (NRG1) can serve as an HSCR susceptibility locus via mediating the development of the enteric ganglia precursors. [7, 8] The ablation of glial cell line-derived neurotrophic factor (GDNF) family receptor α 1 (GFRα 1 ), which is a high affi nity receptor for GDNF, can cause rapid and extensive death of enteric neurons and result in colon aganglionosis reminiscent of HSCR in mice during late gestation. [9] Neurotrophin 3 (NTF3) and neurotrophic tyrosine kinase receptor 3 (NTRK3) may function in the pathogenesis of HSCR by affecting the development of the enteric nervous system (ENS) through the NTF3/NTRK3 signaling pathway. [10, 11] Duan et al [12] found mutations of the endothelin 3 (EDN-3) and endothelin receptor type B (EDNRB) genes in short-segment HSCR in the Chinese population and determined that these genes may be correlated with the pathogenesis of HSCR. Cantrell et al [13] reported that interaction between sex determining region Y (SRY)-box 10 (Sox10) and endothelin receptor type B (EdnrB) affects the development of the ENS in mouse models and can affect the epistatic network, producing phenotypic variation among patients with aganglionosis. RNA-sequencing (or massively parallel cDNA sequencing), which can detect gene changes, gene fusion, and somatic mutation, is an accurate method for quantifying and characterizing transcriptomes. [14] Using RNA-sequencing data, Saeed et al [15] analyzed differentially expressed genes (DEGs) between abnormal and normal bowel segments in patients with HSCR relative to controls and found that reelin (RELN), galanin/GMAP prepropeptide (GAL), growth associated protein 43 (GAP43), Neurensin-1 (NRSN1), and gamma 2 subunit of the gamma-aminobutyric acid receptor type A (GABRG2) correlate with the pathogenesis of HSCR. It is proposed that different analysis processes and thresholds can generate different results. [16] To gain a more comprehensive understanding of the pathogenesis of HSCR, RNA-sequencing data were used to analyze DEGs between expansion and narrow segment samples. In addition, enrichment analysis, protein-protein interaction (PPI) network construction and module analysis were carried out separately to screen critical genes in HSCR. This study was approved by the Institutional Review Board of the Second Affiliated Hospital of Xi'an Jiaotong University. In addition, all participants provided their informed consent. Using electronic colonoscopy, expansion segments and narrow segments of colon tissues were excised from children with HSCR who were treated in the Second Affiliated Hospital of Xi'an Jiaotong University from September 2014 to July 2015. The expansion segments and the narrow segments were randomly mixed in their respective groups to form 2 sequencing samples. All the patients were diagnosed with HSCR by barium enema and anorectal manometry evaluation before surgical treatment. Before surgery, a suction rectal biopsy was obtained to verify aganglionosis. After surgery, pathological staining was performed to establish a definitive diagnosis. The isolated tissues were cut into small pieces (<0.5 cm) and soaked in 5 times the volume of RNAlater (Ambion, Foster City, CA, USA) at 4°C overnight then stored at -80°C. Total RNA was extracted from the expansion and the narrow segments of colon tissues using TRIzol reagent (Invitrogen, Shanghai, China) according to the manufacturer's instructions. The integrity and purity of total RNA were detected using 2% agarose gel electrophoresis and spectrophotometry (Merinton, Beijing, China), respectively. Finally, RNA-seq library construction and high-throughput sequencing were separately carried out using a NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (New England Biolabs, Ipswich, MA, USA) and an Illumina Hiseq 2500 v4 100PE System (Illumina, San Diego, CA, USA). The NGSQC Toolkit [17] was used to preprocess the pairedend sequencing data and filter out low quality sequences. Using TopHat2 software, [18] sequence alignment of the preprocessed data was performed with the hg19 human genome sequences from the UCSC website (http:// genome.ucsc.edu) as references. [19] Subsequently, the expression levels of genes were calculated using Cuffdiffs software [20] to obtain a gene expression matrix of samples. The "edgeR" package [21] in R was used to screen for DEGs between the expansion and the narrow segment samples. The thresholds were P value <0.05 and |log 2 fold change (FC)| >1. Gene Ontology (GO, http://www.geneontology.org/) seeks to develop and use a series of controlled, structured vocabularies for annotating sequences, genes and their products. [22] The Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www. genome.jp/kegg/) aims to integrate biological systems information of genomic, chemical and systemic functional aspects. [23] Using DAVID software, GO functional and KEGG pathway enrichment analyses of DEGs were performed. Terms involving more than 2 genes and with a P value <0.05 were selected. The Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org/) database gives a critical evaluation and integration of PPIs consisting of functional and physical associations. [24] Using the STRING database, [24] the PPI relationships among the DEGs were analyzed, with the required confidence (combined score) >0.4 as the threshold. After the PPIs were searched, the PPI network World J Pediatr, Vol 13 No 3 . June 15, 2017 . www.wjpch.com was constructed using Cytoscape software (http://www. cytoscape.org/). [25] In the PPI network, nodes represented proteins, and their degrees were equal to the number of edges linked to them. Moreover, module analysis was performed for the PPI network using the MCODE plugin [26] in Cytoscape with the default parameters. DEGs with a P value <0.05 and |log 2 FC| >1 between the expansion and the narrow segment samples were analyzed. Compared with narrow segment samples, 117 DEGs were identified in the expansion segment samples, including 47 up-regulated and 70 downregulated genes. Using the DAVID software, enrichment analyses were carried out on the up-regulated genes and the downregulated genes separately. The up-regulated genes were significantly enriched in 43 GO terms, and the top 10 GO terms are listed in Table 1 Table 2) . The PPI network of the DEGs had 59 nodes and 113 interactions (Fig. 1 ). In particular, FOS (degree=20), early growth response 1 (EGR1, degree=16), activating transcription factor 3 (ATF3, degree=9), nitric oxide synthase 1 (NOS1, degree=8), CCL5 (degree=8), DUSP1 (degree=7), CXC chemokine ligand 3 (CXCL3, degree=6), vasoactive intestinal peptide (VIP, degree=6), FBJ murine osteosarcoma viral oncogene homolog B (FOSB, degree=5), and nitric oxide synthase 2 (NOS2, degree=4) had higher degrees and could interact with others genes in the PPI network. Afterwards, module analysis was performed for the PPI network, and two significant modules (module 1 and module 2) were identified (Fig. 2) . The top 10 enriched functions among the genes in module 1 are shown in Table 3 and include G-protein coupled receptor protein signaling pathway (P<0.001), cell surface receptor linked signal transduction (P<0.001) and taxis (P<0.005), which involved CXCL3 and CCL5. The genes in module 2 were enriched in functions such as regulation of transcription from RNA polymerase II promoter (P<0.001), DNA-dependent regulation of transcription (P<0.001) and regulation of RNA metabolic process (P<0.001) ( Table 3 ). As the main component of the peripheral nervous system, the ENS is a complex network of glia and neurons within the intestinal wall, [27, 28] and a developmental disorder of the ENS has been correlated with HSCR. [29] In this study, 117 DEGs were identified in the expansion segment samples, including 47 up-regulated and 70 downregulated genes. In the PPI network, FOS (degree=20), EGR1 (degree=16), ATF3 (degree=9), NOS1 (degree=8), CCL5 (degree=8), DUSP1 (degree=7), CXCL3 (degree=6), VIP (degree=6), FOSB (degree=5), and NOS2 (degree=4) had higher degrees, which had interactions with other genes. Moreover, two significant modules (module 1 and module 2) were identified in the PPI network. Functional enrichment analysis showed that FOS and DUSP1 were implicated in response to endogenous stimulus. The c-fos gene is related to apoptosis and regeneration of damaged retinal ganglion cells (RGCs), indicating that the c-fos-mediated cascade precisely correlates with apoptosis and that regeneration of RGC plays a critical role in neuronal survival and regeneration. [30] As a member of the Fos family, FosB is up-regulated after an acute increase in central venous pressure in unanesthetized rats. [31] DUSP1 can specially regulate sensory input neurons of the telencephalon and thalamus, which may reduce stimulus-induced expression of immediate early genes such as EGR1. [32] These fi ndings suggest that FOS, FOSB, DUSP1 and EGR1 might be correlated with the development of HSCR. Enrichment analyses indicated that CCL5 was involved in the pathway of prion diseases, while CXCL3 and CCL5 were functionally enriched in taxis. Inhibition of CCL5 by neutralizing mAb decreased the severity of central nervous system (CNS) disease in a viral model of demyelination. [33] CC receptor 5 (CCR5) and its primary ligand CCL5 play critical roles in promoting the transport of leukocytes into the CNS of mice infected by mouse hepatitis virus (MHV). [34] [35] [36] CXCL3 has a low expression level in 12-O-tetradecanoyl phorbol-13-acetate inducible sequence 21 (Tis21)-null cerebellar granule neuron precursor cells (GCPs) of lesions and the external granular layer, and CXCL3 activation induced by Tis21 contributes to migration of GCPs. [37] Therefore, CCL5 and CXCL3 might be implicated in the pathogenesis of HSCR. In addition, as a member of the ATF/cyclic AMP responsive element binding family of transcription factors, ATF3 is activated in motoneurons and sensory neurons in the spinal cord after nerve injury and may be a promising marker of nerve injury in the nervous system. [38] ATF3 can promote neurite elongation of dorsal root ganglion, indicating that ATF3 helps nerve regeneration by improving the intrinsic growth state of injured neurons. [39] [40] [41] [42] VIP is a neurotransmitter in the mucosa, lamina propria, and myenteric plexus of the small and large intestine, which can respond to secretomotor refl exes and activate adenylate cyclase. [43] A previous study has shown that VIP functions in reducing the net intestinal fluid absorption in patients with chemical colitis, which can be completely inhibited by its antagonist [4Cl-d-Phe6, Leu17]VIP. [44] Nitric oxide synthases (NOS) are the predominant enzymes producing nitric oxide (NO) from conversion of L-arginine and generally need NADPH and oxygen co-substrates. [45] Additionally, NO produced by the enteric nervous system is conducive for maintaining gut homeostasis through regulation of vascular tone, gut motility, mucosal barrier function, blood supply, immunity and secretions. [46] These findings suggest that ATF3, VIP, NOS1 and NOS2 might also function in HSCR by affecting the nervous system. In conclusion, based on the bioinformatics analyses, 117 DEGs were identified in the expansion segment samples. Moreover, several genes (including FOS, EGR1, ATF3, NOS1, CCL5, DUSP1, CXCL3, VIP, FOSB, and NOS2) might be involved in the development of HSCR through their impact on the nervous system. However, these predicted results must be further confi rmed through experimental studies, such as quantitative real-time polymerase chain reaction (qRT-PCR) in the future. Pathogenesis of Hirschsprung's disease and its variants: recent progress Symptomatology, pathophysiology, diagnostic work-up, and treatment of Hirschsprung disease in infancy and childhood The developmental etiology and pathogenesis of Hirschsprung disease Hirschsprung's disease in Japan: analysis of 3852 patients based on a nationwide survey in 30 years Notch signaling is required for the maintenance of enteric neural crest progenitors Mutations of the RET gene in isolated and syndromic Hirschsprung's disease in human disclose major and modifi er alleles at a single locus Genome-wide association study identifi es NRG1 as a susceptibility locus for Hirschsprung's disease Association of genetic polymorphisms in the RET-protooncogene and NRG1 with Hirschsprung disease in Thai patients Conditional ablation of GFRα1 in postmigratory enteric neurons triggers unconventional neuronal death in the colon and causes a Hirschsprung's disease phenotype A novel point variant in NTRK3, R645C, suggests a role of this gene in the pathogenesis of Hirschsprung disease NTF-3, a gene involved in the enteric nervous system development, as a candidate gene for Hirschsprung disease Clinical relationship between EDN-3 gene, EDNRB gene and Hirschsprung's disease Interactions between Sox10 and EdnrB modulate penetrance and severity of aganglionosis in the Sox10Dom mouse model of Hirschsprung disease RNA sequencing: advances, challenges and opportunities Identification of novel genes in Hirschsprung disease pathway using whole genome expression study Learning dysregulated pathways in cancers from differential variability analysis NGS QC Toolkit: a toolkit for quality control of next generation sequencing data TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions The UCSC genome browser database Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks Bioinformatics and computational biology solutions using R and Bioconductor The Gene Ontology in 2010: extensions and refi nements KEGG for linking genomes to life and the environment STRING v10: protein-protein interaction networks, integrated over the tree of life A travel guide to Cytoscape plugins An automated method for finding molecular complexes in large protein interaction networks The enteric nervous system: normal functions and enteric neuropathies Developmental determinants of the independence and complexity of the enteric nervous system Hirschsprung disease: a developmental disorder of the enteric nervous system The role of c-fos in cell death and regeneration of retinal ganglion cells FosB expression in the central nervous system following isotonic volume expansion in unanesthetized rats The dusp1 immediate early gene is regulated by natural stimuli predominantly in sensory input neurons Antibody targeting of the CC chemokine ligand 5 results in diminished leukocyte infi ltration into the central nervous system and reduced neurologic disease in a viral model of multiple sclerosis Reduced macrophage infiltration and demyelination in mice lacking the chemokine receptor CCR5 following infection with a neurotropic coronavirus Functional expression of chemokine receptor CCR5 on CD4+ T cells during virus-induced central nervous system disease Functional analysis of the CC chemokine receptor 5 (CCR5) on virus-specific CD8+ T cells following coronavirus infection of the central nervous system Tis 21 knock-out enhances the frequency of medulloblastoma in Patched1 heterozygous mice by inhibiting the Cxcl3-dependent migration of cerebellar neurons Activating transcription factor 3 (ATF3) induction by axotomy in sensory and motoneurons: A novel neuronal marker of nerve injury ATF3 increases the intrinsic growth state of DRG neurons to enhance peripheral nerve regeneration Activating transcription factor 3 and the nervous system The transcription factor ATF-3 promotes neurite outgrowth ATF3 enhances c-Jun-mediated neurite sprouting Effects of vasoactive intestinal peptide, secretin, and related peptides on rat colonic transport and adenylate cyclase activity Inhibitory effect of experimental colitis on fluid absorption in rat jejunum: role of the enteric nervous system, VIP, and nitric oxide Enzymatic function of nitric oxide synthases Importance of NO and its related compounds in enteric nervous system regulation of gut homeostasis and disease susceptibility