Original research••

Identifying the potential transcriptional regulatory network in Hirschsprung disease by integrated analysis of microarray datasets

•,,,,,,,,,,,.

...

Abstract

Objective Hirschsprung disease (HSCR) is one of the common neurocristopathies in children, which is associated with at least 20 genes and involves a complex regulatory mechanism. Transcriptional regulatory network (TRN) has been commonly reported in regulating gene expression and enteric nervous system development but remains to be investigated in HSCR. This study aimed to identify the potential TRN implicated in the pathogenesis and diagnosis of HSCR.

Methods Based on three microarray datasets from the Gene Expression Omnibus database, the multiMiR package was used to investigate the microRNA (miRNA)–target interactions, followed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. Then, we collected transcription factors (TFs) from the TransmiR database to construct the TF–miRNA–mRNA regulatory network and used cytoHubba to identify the key modules. Finally, the receiver operating characteristic (ROC) curve was determined and the integrated diagnostic models were established based on machine learning by the support vector machine method.

Results We identified 58 hub differentially expressed microRNAs (DEMis) and 16 differentially expressed mRNAs (DEMs). The robust target genes of DEMis and DEMs mainly enriched in several GO/KEGG terms, including neurogenesis, cell–substrate adhesion, PI3K–Akt, Ras/mitogen-activated protein kinase and Rho/ROCK signaling. Moreover, 2 TFs (TP53 and TWIST1), 4 miRNAs (has-miR-107, has-miR-10b-5p, has-miR-659-3p, and has-miR-371a-5p), and 4 mRNAs (PIM3, CHUK, F2RL1, and CA1) were identified to construct the TF–miRNA–mRNA regulatory network. ROC analysis revealed a strong diagnostic value of the key TRN regulons (all area under the curve values were more than 0.8).

Conclusion This study suggests a potential role of the TF–miRNA–mRNA network that can help enrich the connotation of HSCR pathogenesis and diagnosis and provide new horizons for treatment.

What is already known on this topic

Hirschsprung disease (HSCR) is one of the common neurocristopathies in children that involves a complex pathogenesis. It is difficult to develop early diagnosis of HSCR, and surgery commonly gives rise to medical complications, especially fatal enterocolitis (about 35% after surgery).

What this study adds

A potential transcription factor–microRNA–mRNA regulatory network was identified as for the key regulons of which the receiver operating characteristic analysis revealed a strong diagnostic value in HSCR.

How this study might affect research, practice or policy

This study suggests a transcriptional regulatory network implicated in the pathogenesis and diagnosis of HSCR, which also provides new horizons and targets for treatment.

Introduction

Hirschsprung disease (HSCR) is one of the common neurocristopathies in children, which is characterized by aganglionosis.1 2 HSCR is primarily treated by surgery to eliminate the aganglionic bowel while commonly giving rise to medical complications, especially fatal enterocolitis (about 35% after surgery),2–4 stool leakage, anastomotic stricture, anastomotic leak with abscess, and chronic constipation. Therefore, detailed pathogenesis and effective alternatives should be developed.

At present, it is well known that the pathogenesis of HSCR is the dysfunction of enteric neural crest-derived precursors migrating through the bowel in a rostral-to-caudal direction from week 3 to week 8 of human gestation.2 Emerging studies have reported the effects of enteric neural crest-derived cell (ENCC) transplantation for treating the HSCR model.5–7 However, because of the limited proliferation, migration and large-scale apoptosis during transplantation, ENCC transplantation often tends to be an insufficient cure for HSCR.1 Although researchers have tried the ENCCs treated with cytokines, drugs, and signaling pathway regulators to optimize cell transplantation, it failed to completely repair the enteric nervous system (ENS).8 9 As supposed, HSCR is associated with at least 20 genes of more than seven chromosomal loci, involving a complex regulatory to ENCCs, but not single genetic factors.2 10 11 Therefore, it is necessary to explore more details of the gene expression regulatory in HSCR.

Previous studies have shown that microRNAs (miRNAs) bind on the 5′ untranslated regions of mRNAs through partial complementarity and reduce gene expression by restraining mRNA translation and/or facilitating mRNA degradation.12 Many miRNAs have been reported to be related to HSCR,13–15 such as miRNA-206,16miR-146b-5p,17 and miR-181a.18 Like the functional genes, miRNA expression is regulated by transcription factors (TFs). Transcriptional regulatory network (TRN), demonstrating the relationship of TF–miRNA–mRNA, commonly plays roles in the regulation of gene expression and cell biological function,19–21 and has been reported in ENS development,22 neural stem cell phenotype,20 and cancer pathogenesis.23 However, the role of TRN in HSCR remains to be investigated.

In this study, we performed integrated analysis of three microarray datasets from the Gene Expression Omnibus (GEO) database, based on which a potential TF–miRNA–mRNA network was constructed. Receiver operating characteristic (ROC) analysis based on the support vector machine (SVM) method revealed a strong diagnostic value of the key TRN regulons, which can help enrich the connotation of HSCR pathogenesis and diagnosis and provide new horizons for treatment.

Materials and methods

Microarray datasets and processing

The mRNA and miRNA expression profiles of patients with HSCR were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which was searched using the following terms: “Hirschsprung disease” AND “microarray” AND “Homo sapiens”. The following eligibility criteria were used to include or exclude datasets and samples: (1) the dataset contained at least three patients with HSCR and three controls; (2) the colons from HSCR and normal subjects were used for microarray analysis; and (3) raw data were available in the GEO database. Detailed information of the microarray datasets is listed in table 1.

Table 1

Characteristics of three microarray datasets included in the study

The probe sets were also downloaded from the GEO database, and probes matching with multiple gene symbols were eliminated, while the mean values were calculated for gene symbols corresponding to multiple probes. The differentially expressed microRNAs (DEMis) and the differentially expressed mRNAs (DEMs) between HSCR and control samples in each dataset were identified by the Linear Models for Microarray Data (limma) package V.3.46.024 with the cut-off criteria of |log2 fold change|>0.5 and p value of <0.05. The Venn diagram was used to obtain the common DEMs between the two mRNA microarray datasets.

Hub DEMi identification

The miRNA similarity database (MISIM V.2.0, http://www.lirmed.com/misim/)25 was searched to recognize hub DEMis according to the MISIM V.2.0 Tutorial (http://www.lirmed.com/misim/Help).

miRNA–target interaction investigation

The multiMiR package V.1.20.026 was used to investigate the miRNA–target interactions. This package is a collection of miRNAs/targets from 14 external resources, including three validated miRNA–target databases (miRecords, miRTarBase, and TarBase) and eight predicted micRNA–target databases (DIANA-microT, ElMMo, MicroCosm, miRanda, miRDB, PicTar, PITA, and TargetScan), and so on, which can be used to retrieve all the validated and predicted target genes of a given miRNA, and all the validated and predicted miRNA–target interactions between a set of given miRNAs and mRNAs. Meanwhile, the top ten ranked miRNA–target couples were identified by Maximal Clique Centrality (MCC) algorithm via Cytoscape software V.3.8.2.

Protein–protein interaction (PPI) network analysis

All the target genes of hub DEMis identified previously were uploaded to the STRING database V.11.5 (https://www.string-db.org/)27 to construct the PPI network. Confidence of >0.4 was set as the screening criteria. The PPI network was subsequently reconstructed and visualized by Cytoscape software V.3.8.2. The robust target genes were subsequently screened out using the cytoHubba plugin,28 which investigates the most important nodes in the PPI network with several topological analysis algorithms.

Robust Rank Aggregation (RRA) analysis

To minimize the bias and inconsistencies, we integrated the top 20 ranked genes in the PPI network calculated by eight different topological analysis algorithms (MCC, MNC, EPC, EcCentricity, DMNC, Degree, Closeness, and BottleNeck method), and the RRA package V.1.129 was adopted to identify the robust target genes. The score in the RRA analysis result indicated the ranking degree of each gene in the gene list, and the genes with a score of <0.05 were considered as the robust target genes.

Functional and pathway enrichment analyses

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to investigate the biological process, cellular component, molecular function, and involved pathways of selected genes, which were performed with the clusterProfiler R package V.4.6.0.30 The GO/KEGG terms with an adjusted p value of <0.05 were considered statistically significant and were visualized via the ggplot2 R package V.3.3.3.

TF–miRNA–mRNA regulatory network analysis

The TF–miRNA regulations database (TransmiR V.2.0, http://www.targetscan.org/vert_72/)31 was searched to collect TFs of given miRNAs. Only the validated TF–miRNA interactions were included to construct the TF–miRNA–mRNA regulatory network, in which the key TRN regulon module was identified by the CytoHubba plugin.28 Moreover, the potential TF–miRNA interactions were further analyzed in the University of California Santa Cruz (UCSC) genome browser (https://genome.ucsc.edu/), and the TF–mRNA correlation in the colon was further analyzed in the Gene Regulatory Network Database (GRNdb, http://www.grndb.com/).32

Diagnostic analysis of the key TRN regulons in HSCR

The ROC curve was obtained by GraphPad Prism software V.8.0.1 to assess the accuracy of each key TRN regulon as biomarkers in predicting HSCR. The machine learning based on the SVM method was used to establish an integrated diagnostic model followed by the ROC curve.

Statistical analysis

Statistical analysis was performed by GraphPad Prism software V.8.0.1. Normally distributed data were presented as means±standard deviation (SD), and two-tailed Student’s t-test was applied to compare differences between groups. Statistical significance was set at a p value of <0.05.

Results

Microarray datasets and the workflow of this study

The microarray datasets derived from patients with HSCR were obtained from the GEO database. Only the databases with the normal subjects for control were included for further analysis, including two mRNA microarray datasets (GSE96854 and GSE98502) and one miRNA microarray dataset (GSE77296). The workflow of the study is shown in figure 1. Detailed information of the three datasets is shown in table 1.

Request permission

Figure 1

The whole study workflow. DEM, differentially expressed mRNA; DEMi, differentially expressed microRNA; GEO, Gene Expression Omnibus; GRNdb, Gene Regulatory Network Database; miRNA, microRNA; MISIM V.2.0, miRNA Similarity Database V.2.0; PPI, protein–protein interaction; ROC, receiver operating characteristic; RRA, Robust Rank Aggregation; SVM, support vector machine; TransmiR V.2.0, Transcription Factor Micro-RNA Regulations Database V.2.0; UCSC, University of California Santa Cruz; TRN, transcriptional regulatory network.

Identification of hub DEMis in HSCR

The miRNA microarray dataset (GSE77296) was analyzed by the limma package to identify DEMis of the colon between patients with HSCR and healthy controls. When setting the cut-off criteria as follows: p value of <0.05 and |log2 fold change|>0.5, we obtained 104 DEMis (including 21 upregulated and 83 downregulated DEMis) (figure 2A). Then, we searched the miRNA similarity database (MISIM V.2.0, http://www.lirmed.com/misim/) to recognize hub DEMis, generating 12 upregulated (figure 2B) and 46 downregulated (figure 2C) miRNAs, all of which were illustrated as heatmap (figure 2D) and detailed in online supplemental table 1.

Request permission

Figure 2

Identification of hub DEMis in HSCR. (A) Volcano plot of miRNA microarray dataset GSE77296. The 21 upregulated miRNAs are marked in red; the 83 downregulated miRNAs are marked in blue; and the gray dots represent miRNAs with no significant difference. Network of miRNAs interaction were searched in the MISIM V.2.0 to recognize upregulated (B) and downregulated (C) hub DEMis. (D) Heatmap diagram of the hub DEMis. DEMi, differentially expressed microRNA; HSCR, Hirschsprung disease; miRNA, microRNA; MISIM V.2.0, miRNA Similarity Database V.2.0.

Investigation and functional annotation of the genes targeted by hub DEMis

The multiMiR package was used to investigate the genes targeted by hub DEMis. The target genes shared in three validated databases or at least six predicted databases were chosen in subsequent analysis, including 31 validated and 75 predicted target genes of upregulated miRNAs, while 25 validated and 102 predicted target genes of downregulated miRNAs (marked with a red box in figure 3).

Request permission

Figure 3

Investigation of the genes targeted by hub DEMis. The multiMiR package was used to investigate the miRNA–target interactions. Upset diagram of target genes of upregulated hub DEMis in validated (A) and predicted (B) miRNA–target databases. Upset diagram of target genes of downregulated hub DEMis in validated (C) and predicted (D) miRNA–target databases. The genes marked with red boxes were shared targets in three validated databases or at least six predicted databases. DEMi, differentially expressed microRNA; miRNA, microRNA.

After removing duplicates, 197 target genes were uploaded to the STRING database (http://string.embl.de/) to perform PPI analysis. Then, to hide the disconnected nodes, the Cytoscape software was adopted to visualize the network (figure 4A). Robust target genes were subsequently screened out using the cytoHubba plugin, which investigates the most important nodes in the PPI network with several topological analysis algorithms. To improve the positive rate, the RRA method was used to integrate the top 20 ranked genes calculated by eight different topological analysis algorithms (MCC, MNC, EPC, EcCentricity, DMNC, Degree, Closeness, and BottleNeck), and a total of 14 genes were obtained accordingly (figure 4B). The upset diagram of the top 20 ranked genes from the eight algorithms is shown in online supplemental figure 1. Finally, GO/KEGG functional analysis was performed to explore the biological classifications of robust target genes in HSCR by the clusterProfiler package (figure 4C). GO enrichment analyses showed that the significantly enriched terms were related to the following: neurogenesis; cell cycle, apoptosis, differentiation, aging, and cell–substrate adhesion; protein phosphorylation; protein kinase activity; cellular response to transforming growth factor beta stimulus and vascular endothelial growth factor stimulus; DNA-binding TF activity, etc. In the KEGG pathway analysis, the significantly enriched terms were PI3K–Akt, mitogen-activated protein kinase (MAPK) (ERK1/2), notch, relaxin, and HIF-1 signaling pathway. RAS/MAPK and PI3K–Akt had been reported as the key signaling pathways in neurogenesis and neuroprotection20 33–36 and were related to RET and RET-regulating pathways in HSCR.14

Request permission

Figure 4

Functional annotation of the genes targeted by hub DEMis. (A) Whole PPI network with all target genes of hub DEMis; the bigger dots and deeper represent the higher degree. (B) The lollipop chart shows all robust target genes identified by the RRA method; the bigger dots represent the higher rank. (C) GO/KEGG functional enrichment analysis. BP, biological process; CC, cellular component; DEMi, differentially expressed microRNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; PPI, protein–protein interaction; RRA, Robust Rank Aggregation.

Taken together, these data indicated that the hub DEMis and their target genes identified previously contributed to the pathogenesis of HSCR.

Investigation and functional annotation of DEMs in HSCR

We further analyzed the other two mRNA datasets (GSE96854 and GSE98502) to identify the DEMs in the colon between patients with HSCR and healthy controls. When setting the cut-off criteria as follows: p value of <0.05 and |log2 fold change|>0.5, we obtained 3998 DEMs (including 2253 upregulated and 1745 downregulated DEMs) in GSE96854 and 219 DEMs (including 147 upregulated and 72 downregulated DEMs) in GSE98502 (figure 5A). The common DEMs in the two datasets (including 11 upregulated and 5 downregulated genes) (figure 5B) are detailed in table 2, which were significantly enriched in the GO/KEGG terms of Rho protein signal transduction, Ras protein signal transduction, IKappaB kinase (IKK)/nuclear factor kappa B (NF-κB), and cytokine-mediated signaling pathway (interferon-gamma, interleukin-5, interleukin-10, etc) (figure 5C). Various studies have shown that Rho/ROCK,37–39 RAS/MAPK,20 33 40 and IKK/NF-κB20 41 signaling played crucial roles in neurogenesis, which suggests the significant roles of the common DEMs in HSCR pathogenesis.

Request permission

Figure 5

Investigation and functional annotation of the DEMs in HSCR. (A) Volcano plot of mRNA microarray datasets GSE96854 and GSE98502; the upregulated mRNAs are marked in red; the downregulated mRNAs are marked in blue; and the gray dots represent mRNAs with no significant difference. (B) Venn diagram demonstrates the common 11 upregulated and 5 downregulated DEMs. (C) GO/KEGG functional enrichment analysis. DEM, differentially expressed mRNA; GO, Gene Ontology; HSCR, Hirschsprung disease; KEGG, the Kyoto Encyclopedia of Genes and Genomes.

Table 2

Characteristics of the 16 common differentially expressed mRNAs

Analysis of TF–miRNA–mRNA regulatory network

For a robust miRNA–target interaction, we investigated the hub DEMis and the DEMs shared in two databases by the multiMiR package. A total of 34 miRNA–target couples were identified, including 7 upregulated miRNA-downregulated mRNAs (2 validated and 5 predicted miRNA–target couples) and 27 downregulated miRNA-upregulated mRNA interactions (8 validated and 19 predicted miRNA–target couples) (figure 6A), all of which are detailed in online supplemental table 2. The top 10 ranked miRNA–target couples were identified by the MCC algorithm (figure 6B).

Request permission

Figure 6

Investigation of miRNA–target interactions. (A) miRNA–target interactions. The miRNAs are marked as diamonds, and mRNAs are marked as ellipses; upregulated genes are marked in red, while the downregulated genes are marked in green; bigger nodes indicate the higher degree; red and green lines represent the validated and predicted miRNA–target couples, respectively. (B) Top 10 ranked miRNA–target couples identified by MCC algorithm. The miRNAs and mRNAs are marked as diamonds and ellipses, respectively; sizes and colors of nodes represent the degree in the network. miRNA, microRNA; MCC, Maximal Clique Centrality.

Then, we searched the TF–miRNA regulations database (TransmiR V.2.0, http://www.cuilab.cn/transmir) for the TFs that target the miRNAs in figure 6B. Only the validated TF–miRNA interactions were included to construct the TF–miRNA–mRNA regulatory network (figure 7A). The cytoHubba was used to identify the key modules (figure 7B), which included 2 TFs (TP53 and TWIST1), 4 miRNAs (has-miR-107, has-miR-10b-5p, has-miR-659-3p, and has-miR-371a-5p), and 4 mRNAs (PIM3, conserved helix–loop–helix ubiquitous kinase (CHUK), F2RL1, and CA1). Finally, the potential TF–miRNA interactions were further analyzed in the UCSC genome browser (https://genome.ucsc.edu/) (figure 7C), that is, the promoter region analysis of miRNA genes, showing that a higher level of H3K4me3 methylated modification represents the more reliable TF–target relationship. The correlation analysis of TFs and miRNA-targeted mRNAs in the colon was further analyzed in the GRNdb (http://www.grndb.com/) (figure 7D).

Request permission

Figure 7

Analysis of TF–miRNA–mRNA network. (A) The TF–miRNA–mRNA network. The miRNAs, mRNAs, and TFs are marked as diamonds, ellipses, and octagons, respectively; upregulated genes are marked in red, while the downregulated genes are marked in green; bigger nodes indicate the higher degree; red, gray, and green lines represent the validated, reported, and predicted connections, respectively. (B) Key modules of TF–miRNA–mRNA network identified by cytoHubba. The miRNAs, mRNAs, and TFs are marked as diamonds, ellipses, and octagons, respectively; sizes and colors of nodes represent the degree in the network. (C) Promoter region analysis of miRNA genes in UCSC genome browser (https://genome.ucsc.edu/). Higher level of H3K4me3 methylated modification represents the more reliable TF–target relationship. (D) Correlation analysis of TFs and miRNA-targeted mRNAs in the GRNdb (http://www.grndb.com/). GRNdb, Gene Regulatory Network Database; miRNA, microRNA; TF, transcription factor; UCSC, University of California Santa Cruz.

Diagnostic value of the key TRN regulons as biomarkers in HSCR

The gold standard for the diagnosis of HSCR is rectal mucosal aspiration biopsy and pathological diagnosis, which are commonly invasive and difficult to perform.1 2 The relative expression of the key TRN regulons (has-miR-107, has-miR-10b-5p, has-miR-659-3p, has-miR-371a-5p, PIM3, CHUK, F2RL1, and CA1) were visualized as boxplots (figure 8A,B). To investigate the diagnostic value of these regulons in HSCR, the ROC curve was used, which showed that all eight regulons had area under the curve (AUC) values more than 0.8, indicating a strong diagnostic value (figure 8C). For better diagnosis prediction, these eight regulons were integrated to establish a multimarker diagnosis model based on machine learning by the SVM method. The ROC curve showed that the multimarker models could effectively predict HSCR (AUC=1.00) (figure 8C).

Request permission

Figure 8

Diagnostic value of the key TRN regulons as biomarkers in HSCR (A,B) Relative expression of the eight key TRN regulons (four miRNAs and four mRNAs) in HSCR. (C) ROC curve of the eight key TRN regulons and the integrated diagnosis model based on machine learning by SVM method. AUC, area under the curve; HSCR, Hirschsprung disease; miRNA, microRNA; ROC, receiver operating characteristic; SVM, support vector machine; TRN, transcriptional regulatory network; 95% CI, 95% confidence interval.

Discussion

The transplantation of ENCCs to induce enteric neurogenesis is a potential radical strategy for HSCR while generating insufficient efficacy. It may due to the complex genes regulatory to ENCCs in children with HSCR.1 5 9 Although many genes have been identified to be associated with HSCR,2 11 such as RET, EDNRB, RARB, GATA2, and SOX10, which commonly regulate ENCCs during the development of ENS, how the TRN contributes to HSCR pathogenesis remains to be investigated. This study identified a potential TF–miRNA–mRNA network, including the key regulons of two TFs (TP53 and TWIST1), four miRNAs (has-miR-107, has-miR-10b-5p, has-miR-659-3p, and has-miR-371a-5p), and four mRNAs (PIM3, CHUK, F2RL1, and CA1), that can help enrich the connotation of HSCR pathogenesis and diagnosis and provide new horizons for treatment.

Many miRNAs have been reported to be related to HSCR,13–15 including miRNA-206/SDPR,16 42miR-146b-5p/RET,17 and miR-181a/RAP1B.18 In this study, we found that has-miR-107, has-miR-10b-5p, has-miR-659-3-p, and has-miR-371a-5p were related to HSCR and exerted good diagnostic value. As reported, has-miR-107 regulated Wnt/β-catenin signaling43 and attenuated neurotoxicity induced by 6-hydroxydopamine.44MiR-10b-5p contributed to neurodegenerative disease, diabetes with dysfunction of interstitial Cajal cells, and neuroprotection for hippocampal neuronal cells.45–48 In cancer diseases, miR-659–3p and miR-371a-5p could regulate tumor progression and were associated with chemotherapy resistance.49–53 Novel research has shown that specific miRNAs in serum or plasma exosomal were identified to have good diagnostic value in HSCR.54 55 As mentioned previously, the miRNAs identified in this study had AUC values of more than 0.8 and remained unclear so far in HSCR, which provided new cues for future biomarker study of HSCR treatment and diagnosis.

As reported, approximately 50% of familial and 20% of patients with sporadic HSCR had RET expression abnormalities; 5% of patients had EDNRB variations, while 4% of patients had SOX10 variations. It seems to be difficult to diagnose HSCR by any one of the known pathogenic genes due to the complex non-Mendelian inheritance. In this study, we constructed a potential TF–miRNA–mRNA network, of which a key module with four functional genes (PIM3, CHUK, F2RL1, and CA1) was identified. Based on the key regulons, we constructed a multimarker model by the SVM method, which had an AUC equal to 1 to effectively predict HSCR. It has been reported that PIM3, a proto-oncogene with serine/threonine kinase activity, could regulate cell migration and apoptosis via PI3K–AKT, p38, or Rho GTPase signaling,56–58 and was related to demyelinating disease.59 Inhibitor-κB kinase α, which is encoded by the CHUK gene, was recognized to regulate NF-κB activity60 61 and involved the differentiation of mouse embryonic neuroectoderm. F2RL1 was reported as the key protease-activated receptor to stimulate neuronal repair after ischemic injury.62 63 The GO/KEGG annotations of carbonic anhydrase 1 (CA1) were carbonate dehydratase activity, hydrolyase activity and interleukin-12 family signaling. At present, all the aforementioned genes were still unclear but relevant to neuropathies, especially HSCR.

As reported, the development and functional maturity of ENS is regulated by complex mechanisms, which largely depend on the potential of ‘seed’ ENCCs and their compatibility with the intestinal microenvironment ‘niche’.64 65 The genetic factors, such as gene mutations (including RET, EDNRB, RARB, GATA2, SOX10, PHOX2B, etc)2 65 and signaling pathway disorders (including PI3K–Akt, MAPK, IKK/NF-κB, Rho/ROCK, etc), determine the inborn developmental potential of ENCCs. Meanwhile, the critical role of intestinal microenvironment, such as glial cell line-derived neurotrophic factor, 5-hydroxytryptamine, semaphorins, neuregulin 1, the extracellular matrix molecules (collagen, laminin, proteoglycans, etc),65 postnatal intestinal flora colonization, and their metabolites,66 has been gradually recognized. Although the functional annotation of the TRN regulons mentioned previously appeared to be associated with the signaling pathways in neurogenesis and neuroprotection, which suggests the significant roles in HSCR pathogenesis, how the TRN regulons regulate the ENCCs and interact with these intestinal microenvironment niche remain to be further investigated.

In conclusion, this study provided a potential TF–miRNA–mRNA network based on integrated analysis of three microarray datasets. ROC analysis based on the SVM method revealed a strong diagnostic value of the key TRN regulons, which can help enrich the connotation of HSCR pathogenesis and diagnosis and provide new horizons for further study. However, due to the limited datasets of HSCR, an integrated model containing miRNAs and mRNA to predict HSCR was unavailable. Moreover, further validated experiments with cells and animals were extensible.

WX and HY contributed equally.
Contributors: WX conceptualized the study design; acquired, analyzed, and visualized the data; interpreted the results; and wrote the manuscript with all authors providing feedback for revision. HY acquired, analyzed, and visualized the data and interpreted the results. DC acquired, analyzed, and visualized the data. WP, WY, JM, and BZ contributed to the data analysis, visualization, and curation. WJ, YL and XC performed the supervision. YG wrote the manuscript, with all authors providing feedback for revision. DT, as the guarantor of the study, conceptualized the study design, interpreted the results, and wrote the manuscript with all authors providing feedback for revision. All authors read and approved the final report.
Funding: The study was supported by grants from the National Natural Science Foundation of China (numbers 82071692, 81770513, and 82170531); Xi’an Jiaotong University (number YXJLRH2022053); and the General Project of Shaanxi Science and Technology Department (number 2022SF-133/033).
Competing interests: None declared.
Provenance and peer review: Not commissioned; externally peer reviewed.
Supplemental material: This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.

Data availability statement

Data are available in a public, open access repository. Publicly available datasets (GSE96854, GSE98502, and GSE77296) were analyzed in this study. All the datasets can be found in the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/).

Ethics statements

Patient consent for publication:

Ethics approval:

Not applicable.

Acknowledgements

We sincerely thank the scientists who shared their data on the public database. We thank the creators of the packages multiMiR, limma, RRA, clusterProfiler, and ggplot2, and the databases including GEO, MISIM V.2.0, STRING V.11.5, TransmiR V.2.0, UCSC, and Gene Regulatory Network Database. We thank Dr Weifeng Hong, Yuzhou Xue, and Zhenlu Cai for reviewing, supporting, and providing valuable comments.

Pan W, Goldstein AM, Hotta R, et al. Opportunities for novel diagnostic and cell-based therapies for Hirschsprung disease. J Pediatr Surg 2022; 57:61–8.
doi:10.1016/j.jpedsurg.2021.10.049•Google Scholar
Heuckeroth RO. Hirschsprung disease-integrating basic science and clinical medicine to improve outcomes. Nat Rev Gastroenterol Hepatol 2018; 15:152–67.
doi:10.1038/nrgastro.2017.149•Google Scholar
Haricharan RN, Georgeson KE. Hirschsprung disease. Semin Pediatr Surg 2008; 17:266–75.
doi:10.1053/j.sempedsurg.2008.07.005•Google Scholar
Kim AC, Langer JC, Pastor AC, et al. Endorectal pull-through for Hirschsprung’s disease-a multicenter, long-term comparison of results: transanal vs transabdominal approach. J Pediatr Surg 2010; 45:1213–20.
doi:10.1016/j.jpedsurg.2010.02.087•Google Scholar
Burns AJ, Goldstein AM, Newgreen DF, et al. White paper on guidelines concerning enteric nervous system stem cell therapy for enteric neuropathies. Dev Biol 2016; 417:229–51.
doi:10.1016/j.ydbio.2016.04.001•Google Scholar
Workman MJ, Mahe MM, Trisno S, et al. Engineered human pluripotent-stem-cell-derived intestinal tissues with a functional enteric nervous system. Nat Med 2017; 23:49–59.
doi:10.1038/nm.4233•Google Scholar
Yu H, Zheng B-J, Pan W-K, et al. Combination of exogenous cell transplantation and 5-HT4 receptor agonism induce endogenous enteric neural crest-derived cells in a rat hypoganglionosis model. Exp Cell Res 2017; 351:36–42.
doi:10.1016/j.yexcr.2016.12.022•Google Scholar
Zhang Y, Seid K, Obermayr F, et al. Activation of wnt signaling increases numbers of enteric neurons derived from neonatal mouse and human progenitor cells. Gastroenterology 2017; 153:154–65.
doi:10.1053/j.gastro.2017.03.019•Google Scholar
Fattahi F, Steinbeck JA, Kriks S, et al. Deriving human ENS lineages for cell therapy and drug discovery in Hirschsprung disease. Nature 2016; 531:105–9.
doi:10.1038/nature16951•Google Scholar
Rogers JM. Search for the missing lncs: gene regulatory networks in neural crest development and long non-coding RNA biomarkers of Hirschsprung’s disease. Neurogastroenterol Motil 2016; 28:161–6.
doi:10.1111/nmo.12776•Google Scholar
Tilghman JM, Ling AY, Turner TN, et al. Molecular genetic anatomy and risk profile of Hirschsprung’s disease. N Engl J Med 2019; 380:1421–32.
doi:10.1056/NEJMoa1706594•Google Scholar
Gebert LFR, MacRae IJ. Regulation of microRNA function in animals. Nat Rev Mol Cell Biol 2019; 20:21–37.
doi:10.1038/s41580-018-0045-7•Google Scholar
Sergi CM, Caluseriu O, McColl H, et al. Hirschsprung’s disease: clinical dysmorphology, genes, micro-RNAs, and future perspectives. Pediatr Res 2017; 81:177–91.
doi:10.1038/pr.2016.202•Google Scholar
Li S, Wang S, Guo Z, et al. Mirna profiling reveals dysregulation of RET and RET-regulating pathways in Hirschsprung’s disease. PLoS One 2016; 11.
doi:10.1371/journal.pone.0150222•Google Scholar
Hong M, Li X, Li Y, et al. Hirschsprung’s disease: key microRNAs and target genes. Pediatr Res 2022; 92:737–47.
doi:10.1038/s41390-021-01872-1•Google Scholar
Budi NYP, Gunadi. Aberrant expressions of mirna-206 target, FN1, in multifactorial Hirschsprung disease. Orphanet J Rare Dis 2019; 14.
doi:10.1186/s13023-018-0973-5•Google Scholar
Xia R-P, Zhao F, Ma T-D, et al. Circ-ITCH overexpression promoted cell proliferation and migration in Hirschsprung disease through miR-146b-5p/RET axis. Pediatr Res 2022; 92:1008–16.
doi:10.1038/s41390-021-01860-5•Google Scholar
Chen G, Peng L, Zhu Z, et al. Lncrna AFAP1-AS functions as a competing endogenous RNA to regulate Rap1b expression by sponging MiR-181a in the HSCR. Int J Med Sci 2017; 14:1022–30.
doi:10.7150/ijms.18392•Google Scholar
Chatterjee S, Kapoor A, Akiyama JA, et al. Enhancer variants synergistically drive dysfunction of a gene regulatory network in Hirschsprung disease. Cell 2016; 167:355–68.
doi:10.1016/j.cell.2016.09.005•Google Scholar
Yun SJ, Byun K, Bhin J, et al. Transcriptional regulatory networks associated with self-renewal and differentiation of neural stem cells. J Cell Physiol 2010; 225:337–47.
doi:10.1002/jcp.22294•Google Scholar
Blais A, Dynlacht BD. Constructing transcriptional regulatory networks. Genes Dev 2005; 19:1499–511.
doi:10.1101/gad.1325605•Google Scholar
Memic F, Knoflach V, Morarach K, et al. Transcription and signaling regulators in developing neuronal subtypes of mouse and human enteric nervous system. Gastroenterology 2018; 154:624–36.
doi:10.1053/j.gastro.2017.10.005•Google Scholar
Anastasiadou E, Jacob LS, Slack FJ, et al. Non-coding RNA networks in cancer. Nat Rev Cancer 2018; 18:5–18.
doi:10.1038/nrc.2017.99•Google Scholar
Ritchie ME, Phipson B, Wu D, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015; 43.
doi:10.1093/nar/gkv007•Google Scholar
Li J, Zhang S, Wan Y, et al. MISIM v2.0: a web server for inferring microRNA functional similarity based on microrna-disease associations. Nucleic Acids Res 2019; 47:W536–41.
doi:10.1093/nar/gkz328•Google Scholar
Ru Y, Kechris KJ, Tabakoff B, et al. The multimir R package and database: integration of microRNA-target interactions along with their disease and drug associations. Nucleic Acids Res 2014; 42.
doi:10.1093/nar/gku631•Google Scholar
Szklarczyk D, Gable AL, Nastou KC, et al. The string database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res 2021; 49:D605–12.
doi:10.1093/nar/gkaa1074•Google Scholar
Chin C-H, Chen S-H, Wu H-H, et al. CytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol 2014; 8 Suppl 4.
doi:10.1186/1752-0509-8-S4-S11•Google Scholar
Kolde R, Laur S, Adler P, et al. Robust RANK aggregation for gene list integration and meta-analysis. Bioinformatics 2012; 28:573–80.
doi:10.1093/bioinformatics/btr709•Google Scholar
Wu T, Hu E, Xu S, et al. ClusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb) 2021; 2.
doi:10.1016/j.xinn.2021.100141•Google Scholar
Tong Z, Cui Q, Wang J, et al. TransmiR v2.0: an updated transcription factor-microrna regulation database. Nucleic Acids Res 2019; 47:D253–8.
doi:10.1093/nar/gky1023•Google Scholar
Fang L, Li Y, Ma L, et al. GRNdb: decoding the gene regulatory networks in diverse human and mouse conditions. Nucleic Acids Res 2021; 49:D97–103.
doi:10.1093/nar/gkaa995•Google Scholar
Wang J, Meng X, Feng C, et al. Benzophenone-3 induced abnormal development of enteric nervous system in zebrafish through MAPK/ERK signaling pathway. Chemosphere 2021; 280.
doi:10.1016/j.chemosphere.2021.130670•Google Scholar
Wang R, Sun Y, Huang H, et al. Rutin, a natural flavonoid protects PC12 cells against sodium nitroprusside-induced neurotoxicity through activating PI3K/Akt/mTOR and ERK1/2 pathway. Neurochem Res 2015; 40:1945–53.
doi:10.1007/s11064-015-1690-2•Google Scholar
Nakano N, Matsuda S, Ichimura M, et al. Pi3K/Akt signaling mediated by G protein‑coupled receptors is involved in neurodegenerative Parkinson’s disease (review). Int J Mol Med 2017; 39:253–60.
doi:10.3892/ijmm.2016.2833•Google Scholar
Furusho M, Ishii A, Bansal R, et al. Signaling by FGF receptor 2, not FGF receptor 1, regulates myelin thickness through activation of ERK1/2-MAPK, which promotes mTORC1 activity in an Akt-independent manner. J Neurosci 2017; 37:2931–46.
doi:10.1523/JNEUROSCI.3316-16.2017•Google Scholar
Koch JC, Tönges L, Barski E, et al. Rock2 is a major regulator of axonal degeneration, neuronal death and axonal regeneration in the CNS. Cell Death Dis 2014; 5.
doi:10.1038/cddis.2014.191•Google Scholar
Zhao Y, Ge X, Yu H, et al. Inhibition of rock signaling pathway accelerates enteric neural crest cell-based therapy after transplantation in a rat hypoganglionic model. Neurogastroenterol Motil 2020; 32.
doi:10.1111/nmo.13895•Google Scholar
Pernet V, Joly S, Jordi N, et al. Misguidance and modulation of axonal regeneration by STAT3 and Rho/ROCK signaling in the transparent optic nerve. Cell Death Dis 2013; 4.
doi:10.1038/cddis.2013.266•Google Scholar
Zhang Y, Smolen P, Baxter DA, et al. Biphasic regulation of p38 MAPK by serotonin contributes to the efficacy of stimulus protocols that induce long-term synaptic facilitation. ENeuro 2017; 4:2017.
doi:10.1523/ENEURO.0373-16.2017•Google Scholar
Long D, Liu M, Li H, et al. Dysbacteriosis induces abnormal neurogenesis via LPS in a pathway requiring NF-κb/IL-6. Pharmacol Res 2021; 167:105543.
doi:10.1016/j.phrs.2021.105543•Google Scholar
Sharan A, Zhu H, Xie H, et al. Down-regulation of mir-206 is associated with hirschsprung disease and suppresses cell migration and proliferation in cell models. Sci Rep 2015; 5.
doi:10.1038/srep09302•Google Scholar
Zhang Z, Wu S, Muhammad S, et al. MiR-103/107 promote ER stress-mediated apoptosis via targeting the wnt3a/β-catenin/ATF6 pathway in preadipocytes. J Lipid Res 2018; 59:843–53.
doi:10.1194/jlr.M082602•Google Scholar
Sun L, Zhang T, Xiu W, et al. MiR-107 overexpression attenuates neurotoxicity induced by 6-hydroxydopamine both in vitro and in vivo. Chem Biol Interact 2020; 315:108908.
doi:10.1016/j.cbi.2019.108908•Google Scholar
Wang L, Liu W, Zhang Y, et al. Dexmedetomidine had neuroprotective effects on hippocampal neuronal cells via targeting lncRNA SHNG16 mediated microrna-10b-5p/BDNF axis. Mol Cell Biochem 2020; 469:41–51.
doi:10.1007/s11010-020-03726-6•Google Scholar
Singh R, Ha SE, Wei L, et al. MiR-10b-5p rescues diabetes and gastrointestinal dysmotility. Gastroenterology 2021; 160:1662–78.
doi:10.1053/j.gastro.2020.12.062•Google Scholar
Ruan Z, Li Y, He R, et al. Inhibition of microrna-10b-5p up-regulates HOXD10 to attenuate Alzheimer’s disease in rats via the Rho/ROCK signalling pathway. J Drug Target 2021; 29:531–40.
doi:10.1080/1061186X.2020.1864739•Google Scholar
Müller S. In silico analysis of regulatory networks underlines the role of mir-10b-5p and its target BDNF in Huntington’s disease. Transl Neurodegener 2014; 3.
doi:10.1186/2047-9158-3-17•Google Scholar
Wang J, Zhang X, Zhang J, et al. Long noncoding RNA CRART16 confers 5-FU resistance in colorectal cancer cells by sponging mir-193b-5p. Cancer Cell Int 2021; 21.
doi:10.1186/s12935-021-02353-5•Google Scholar
Zhang X, Wen L, Chen S, et al. The novel long noncoding RNA CRART16 confers cetuximab resistance in colorectal cancer cells by enhancing ERBB3 expression via miR-371a-5p. Cancer Cell Int 2020; 20.
doi:10.1186/s12935-020-1155-9•Google Scholar
Gong Y, Wei ZR. MiR-659-3p inhibits osteosarcoma progression and metastasis by inhibiting cell proliferation and invasion via targeting SRPK1. BMC Cancer 2022; 22.
doi:10.1186/s12885-022-10029-0•Google Scholar
Li S, Fang Y, Qin H, et al. MiR-659-3p is involved in the regulation of the chemotherapy response of colorectal cancer via modulating the expression of SPHK1. Am J Cancer Res 2016; 6:1976–85.
Google Scholar
Villaruz LC, Huang G, Romkes M, et al. Microrna expression profiling predicts clinical outcome of carboplatin/paclitaxel-based therapy in metastatic melanoma treated on the ECOG-ACRIN trial E2603. Clin Epigenetics 2015; 7.
doi:10.1186/s13148-015-0092-2•Google Scholar
Lv X, Li Y, Li H, et al. Molecular function predictions and diagnostic value analysis of plasma exosomal miRNAs in Hirschsprung’s disease. Epigenomics 2020; 12:409–22.
doi:10.2217/epi-2019-0190•Google Scholar
Tang W, Li H, Tang J, et al. Specific serum microrna profile in the molecular diagnosis of hirschsprung’s disease. J Cell Mol Med 2014; 18:1580–7.
doi:10.1111/jcmm.12348•Google Scholar
Dang Y, Jiang N, Wang H, et al. Proto-oncogene serine/threonine kinase PIM3 promotes cell migration via modulating rho gtpase signaling. J Proteome Res 2020; 19:1298–309.
doi:10.1021/acs.jproteome.9b00821•Google Scholar
Liu D, He M, Yi B, et al. Pim-3 protects against cardiomyocyte apoptosis in anoxia/reoxygenation injury via p38-mediated signal pathway. Int J Biochem Cell Biol 2009; 41:2315–22.
doi:10.1016/j.biocel.2009.05.021•Google Scholar
Cheung CHY, Hsu C-L, Lin T-Y, et al. ZNF322A-mediated protein phosphorylation induces autophagosome formation through modulation of IRS1-AKT glucose uptake and HSP-elicited UPR in lung cancer. J Biomed Sci 2020; 27.
doi:10.1186/s12929-020-00668-5•Google Scholar
McCombe PA, Clark P, Frith JA, et al. Alpha-1 antitrypsin phenotypes in demyelinating disease: an association between demyelinating disease and the allele pim3. Ann Neurol 1985; 18:514–6.
doi:10.1002/ana.410180417•Google Scholar
Li X, Hu Y. Attribution of NF-κB activity to CHUK/ikkα-involved carcinogenesis. Cancers (Basel) 2021; 13.
doi:10.3390/cancers13061411•Google Scholar
Lüningschrör P, Kaltschmidt B, Kaltschmidt C, et al. Knockdown of IKK1/2 promotes differentiation of mouse embryonic stem cells into neuroectoderm at the expense of mesoderm. Stem Cell Rev Rep 2012; 8:1098–108.
doi:10.1007/s12015-012-9402-7•Google Scholar
Wang Y, Zhao Z, Rege SV, et al. 3K3A-activated protein C stimulates postischemic neuronal repair by human neural stem cells in mice. Nat Med 2016; 22:1050–5.
doi:10.1038/nm.4154•Google Scholar
Liu Y, Li H, Hu J, et al. Differential expression and distinct roles of proteinase-activated receptor 2 in microglia and neurons in neonatal mouse brain after hypoxia-ischemic injury. Mol Neurobiol 2022; 59:717–30.
doi:10.1007/s12035-021-02594-5•Google Scholar
Ji Y, Tam P-H, Tang C-M, et al. Roles of enteric neural stem cell niche and enteric nervous system development in Hirschsprung disease. Int J Mol Sci 2021; 22.
doi:10.3390/ijms22189659•Google Scholar
Nagy N, Goldstein AM. Enteric nervous system development: a crest cell’s journey from neural tube to colon. Semin Cell Dev Biol 2017; 66:94–106.
doi:10.1016/j.semcdb.2017.01.006•Google Scholar
Vicentini FA, Keenan CM, Wallace LE, et al. Intestinal microbiota shapes gut physiology and regulates enteric neurons and glia. Microbiome 2021; 9.
doi:10.1186/s40168-021-01165-z•Google Scholar

Overview

Abstract
Introduction
Materials and methods
Results
Discussion
References
Supplementary files
Footnotes
Publication history
Rapid Responses

Article metrics

Altmetric

Dimensions