ChIP-Seq-Based Approach in Mouse Enteric Precursor Cells Reveals New Potential Genes with a Role in Enteric Nervous System Development and Hirschsprung Disease

Hirschsprung disease (HSCR) is a neurocristopathy characterized by intestinal aganglionosis which is attributed to a failure in neural crest cell (NCC) development during the embryonic stage. The colonization of the intestine by NCCs is a process finely controlled by a wide and complex gene regulatory system. Several genes have been associated with HSCR, but many aspects still remain poorly understood. The present study is focused on deciphering the PAX6 interaction network during enteric nervous system (ENS) formation. A combined experimental and computational approach was performed to identify PAX6 direct targets, as well as gene networks shared among such targets as potential susceptibility factors for HSCR. As a result, genes related to PAX6 either directly (RABGGTB and BRD3) or indirectly (TGFB1, HRAS, and GRB2) were identified as putative genes associated with HSCR. Interestingly, GRB2 is involved in the RET/GDNF/GFRA1 signaling pathway, one of the main pathways implicated in the disease. Our findings represent a new contribution to advance in the knowledge of the genetic basis of HSCR. The investigation of the role of these genes could help to elucidate their implication in HSCR onset.


Introduction
Hirschsprung disease (HSCR, OMIM 142623) or aganglionic megacolon, is a disorder characterized by intestinal aganglionosis in a variable segment of the distal gut, whose main manifestation is constipation or intestinal obstruction [1]. The cause of HSCR is attributed to a defective proliferation, differentiation, survival, and/or migration of enteric precursor cells (EPCs) originating from the neural crest between the seventh and twelfth weeks of gestation. HSCR can be classified into three types based on the length of the aganglionic tract, including short-segment HSCR (S-HSCR), long-segment HSCR (L-HSCR), and total colonic/intestinal aganglionosis (TCA/TIA), which represents around 80%, 15%, and 5% of the cases, respectively. In addition, HSCR can appear either sporadically or within a familial aggregation and may be associated with other developmental defects [2]. Alterations in the specific gene expression patterns required at different stages of the enteric nervous system (ENS) formation can lead to HSCR [3]. To date, a large number of different genes have been related to HSCR. In most cases, these genes encode molecules that are widely involved in the major signaling pathways associated with HSCR, the most important being the RET/GDNF/GFRA1 and EDN3/EDNRB signaling pathways [4]. Nevertheless, the genetic cause of the disease in a large percentage of HSCR patients remains unknown; thus, the search for additional genes and mechanisms to better outline HSCR etiology is still required. HSCR, especially the short segment forms, is a disorder with a complex genetic basis that could be explained by the heterogeneity of ENS development, which is regulated by an ever-increasing range of molecules. In accordance with this, the HSCR phenotype is often the result of the presence of various genetic variants, which contribute in an additive manner [5]. Despite the fact that some families show an either dominant or recessive pattern, most of the cases show a complex and non-Mendelian mode of inheritance with variable penetrance [2,5].
One of the best-known mechanisms that controls gene expression during ENS development is the regulation carried out by transcription factors. In this sense, several studies have identified transcription factors expressed during the different stages of the ENS development [6,7]. A few transcription factors have been previously associated with HSCR onset, such as ZFXH1B, SOX10, PHOX2B, and PAX6 [8][9][10][11][12], which modulate expression of different downstream effectors at the developmental stage where EPCs migrate along the gut. Specifically, the expression of PAX6 both at the mRNA and the protein levels in EPCs was previously detected by our group [12,13]. Reduced PAX6 expression levels were detected in HSCR patients and different mechanisms that might lead to this downregulation have been proposed. This study suggested that PAX6 may have a role in the regulation of the transcriptional network involved in ENS development [12]. Therefore, the aims of the present study were to clarify the specific gene expression pattern established by PAX6 for correct ENS development, as well as to identify new susceptibility genes for HSCR. With this purpose, we performed a chromatin immunoprecipitation coupled with a massively parallel sequencing (ChIP-seq) assay to identify PAX6 target genes in mouse EPCs, as well as a differential expression assay in EPCs from HSCR patients versus controls to evaluate their implication in the onset of the disease. In addition, we have searched for rare sequences that may be associated with this pathology within the coding sequence of the resulting genes that we found to be expressed in human EPCs.

Identification of the PAX6 Target Genes in Mouse Neurosphere-Like Bodies (NLBs)
EPCs grow in culture as cell aggregates known as neurosphere-like bodies (NLBs), which are a population of neural crest-derived multipotent stem cells present in postnatal gut. The current study was focused on the identification of PAX6 gene regulatory networks during ENS development to detect additional genes and mechanisms to better outline HSCR etiology. With this aim, a ChIP-seq assay was performed using primary culture NLBs from mouse. The reason for this is that a considerable starting quantity of DNA is needed for this technique and that only one sample for each patient was available.
Upon completion of sequencing, the percentage of reads that passed the internal quality filtering, defined as "pass filter reads" (PF Reads) was >90% (27.2 × 10 6 reads). A total of 16.8 × 10 6 reads mapped with the reference genome. As result, a total of 2247 peaks were detected as well as a total of 4250 genes as the nearest genes to the peaks. Therefore, a methodological approach based on the one described in Villalba-Benito et al. [14] was used for the purpose of selecting a handful of putative targets as well as the most reliable genes suitable to be HSCR candidate genes (pipeline described in the Methods Section). Application of the data analysis pipeline for ChIP-Seq allowed for the identification of 19 PAX6 target genes that resulted from 17 different PAX6 binding regions (Table 1). To further characterize these 17 PAX6 binding regions, we investigated the presence of possible novel PAX6 DNA-binding motifs into these sequences using the MEME-ChIP tool (motif discovery). This analysis revealed two possible biologically relevant motifs that were over-represented (Table 2).

Selection of Genes to Evaluate Their Implication in ENS Development and HSCR
To determine the involvement of PAX6 target genes in human ENS development, we searched for the human orthologs genes of the PAX6 target genes identified in mice by the Mouse Genome Informatics database (MGI). It should be noted that the human orthologs gene for Speer5-ps1, Gm15997, and Wdr95 still remain unknown. Therefore, we eliminated these genes from the analysis. Afterwards, gene interactions of the PAX6 target genes by the Ingenuity Pathway Analysis (IPA) tool were analyzed. Our main aim was to identify new potential susceptibility genes for HSCR and thus, after analyzing them, we decided to select 10 additional genes (TGFB1, VHL, APP, MYC, NTRK1, GRB2, HRAS, HTT, TNF, and HIST1H3A) based on their interaction with several PAX6 target genes (Table S1). Therefore, a total of 26 genes were considered to further assess their potential implication in ENS development and HSCR (16 PAX6 target genes and 10 related genes) ( Table 3). GTATRTGTGTGTGT

Selection of Genes to Evaluate Their Implication in ENS Development and HSCR
To determine the involvement of PAX6 target genes in human ENS development, we searched for the human orthologs genes of the PAX6 target genes identified in mice by the Mouse Genome Informatics database (MGI). It should be noted that the human orthologs gene for Speer5-ps1, Gm15997, and Wdr95 still remain unknown. Therefore, we eliminated these genes from the analysis. Afterwards, gene interactions of the PAX6 target genes by the Ingenuity Pathway Analysis (IPA) tool were analyzed. Our main aim was to identify new potential susceptibility genes for HSCR and thus, after analyzing them, we decided to select 10 additional genes (TGFB1, VHL, APP, MYC, NTRK1, GRB2, HRAS, HTT, TNF, and HIST1H3A) based on their interaction with several PAX6 target genes (Table S1). Therefore, a total of 26 genes were considered to further assess their potential implication in ENS development and HSCR (16 PAX6 target genes and 10 related genes) ( Table 3).

Selection of Genes to Evaluate Their Implication in ENS Development and HSCR
To determine the involvement of PAX6 target genes in human ENS development, we searched for the human orthologs genes of the PAX6 target genes identified in mice by the Mouse Genome Informatics database (MGI). It should be noted that the human orthologs gene for Speer5-ps1, Gm15997, and Wdr95 still remain unknown. Therefore, we eliminated these genes from the analysis. Afterwards, gene interactions of the PAX6 target genes by the Ingenuity Pathway Analysis (IPA) tool were analyzed. Our main aim was to identify new potential susceptibility genes for HSCR and thus, after analyzing them, we decided to select 10 additional genes (TGFB1, VHL, APP, MYC, NTRK1, GRB2, HRAS, HTT, TNF, and HIST1H3A) based on their interaction with several PAX6 target genes (Table S1). Therefore, a total of 26 genes were considered to further assess their potential implication in ENS development and HSCR (16 PAX6 target genes and 10 related genes) ( Table 3).

Gene Expression Patterns in NLBs From HSCR Patients Versus NLBs From Controls
To identify if the selected genes were involved in human ENS development, their gene expression profiles were analyzed using NLBs from the human gut. We observed expression of 11 genes, 5 of them being PAX6 target genes (ACADM, ATXN1, BRD3, COL4A2, RABGGTB) and the remaining genes being related to them (APP, GRB2, HRAS, HTT, MYC, TGFB1) ( Figure S1). It is worth mentioning the identification of Edn3 as a target gene of PAX6 since EDN3 is a member of one of the main pathways implicated in HSCR. We did not detect expression in human NLBs, which agrees with the expression pattern of EDN3 in the surrounding mesenchyme of ENS. However, our interest in this study was not the study of well-known genes related to the disease, rather than detecting new genes involved in such pathology. Thus, we focused on following up the study of new candidate genes.
To assess the possible role of these genes in the pathogenesis of HSCR, their expression levels in NLBs from HSCR patients versus controls were compared. We identified 5 genes with statistically significant different expression levels (GRB2, BRD3, HRAS, RABGGTB, and TGFB1; p-value ≤ 0.05). Specifically, HRAS, RABGGTB, and TGFB1 were upregulated, whereas GRB2 and BRD3 were downregulated in HSCR-NLBs ( Figure 1).  Genes were hierarchically clustered by Pearson correlation coefficient using average linkage. The color scale, representing ΔCt, is shown on the right side. Green indicates genes with relatively decreased expression in HSCR, whereas red indicates genes with relatively increased expression in HSCR compared with the controls. * p value < 0.05, ** p value < 0.01 and *** p value < 0.001.

Identification and Assessment of Susceptibility Rare Variants for HSCR in the Genes Expressed in Human ENS
The sequences of the genes that showed expression in human NLBs (ACADM, ATXN1, BRD3, COL4A2, RABGGTB, APP, GRB2, HRAS, HTT, MYC, TGFB1) were analyzed to identify potential pathogenic variants related to HSCR, through the analysis of the whole exome sequencing (WES) data from 56 HSCR patients [15]. As a result, a set of rare variants were identified in ATXN1, RABGGTB, COL4A2, HTT, APP, and GRB2. All these rare variants were inherited from one of their parents and occasionally they appeared in an unaffected sibling; therefore, it seems that these rare variants do not play an important role in HSCR susceptibility. Such variants are available under request.

Functional Networks Shared among the Possible Susceptibility Genes for HSCR
To determine the functional association networks among the genes identified as putative susceptibility genes in HSCR, we used GeneMANIA software ( Figure S2). As a result, all the genes interacted with each other through physical interactions, co-expression, and pathway. In addition, gene function prediction showed that TGFB1, HRAS, and GRB2 contribute to the ErbB or epidermal

Identification and Assessment of Susceptibility Rare Variants for HSCR in the Genes Expressed in Human ENS
The sequences of the genes that showed expression in human NLBs (ACADM, ATXN1, BRD3, COL4A2, RABGGTB, APP, GRB2, HRAS, HTT, MYC, TGFB1) were analyzed to identify potential pathogenic variants related to HSCR, through the analysis of the whole exome sequencing (WES) data from 56 HSCR patients [15]. As a result, a set of rare variants were identified in ATXN1, RABGGTB, COL4A2, HTT, APP, and GRB2. All these rare variants were inherited from one of their parents and occasionally they appeared in an unaffected sibling; therefore, it seems that these rare variants do not play an important role in HSCR susceptibility. Such variants are available under request.

Functional Networks Shared among the Possible Susceptibility Genes for HSCR
To determine the functional association networks among the genes identified as putative susceptibility genes in HSCR, we used GeneMANIA software ( Figure S2). As a result, all the genes interacted with each other through physical interactions, co-expression, and pathway. In addition, gene function prediction showed that TGFB1, HRAS, and GRB2 contribute to the ErbB or epidermal growth factor receptor (EGFR) signaling pathway. HRAS and TGFB1 also showed functions related to cell cycle and cell migration regulation.

GRB2 Implication in the Onset of HSCR
Based on our results and the fact that GRB2 encodes for the intracellular adaptor protein growth factor receptor-bound 2, which is involved in RET/GDNF/GFRA1 signaling, we decided to further investigate the association of GRB2 with HSCR onset. We performed a mutational screening by direct sequencing of the GRB2 gene, including its coding sequence, intron/exon boundaries and untranslated regions (UTR) in 267 isolated HSCR patients. As a result, we identified several GRB2 rare variants that alone do not seem to significantly contribute to the manifestation of the phenotype, ruling out, therefore, the presence of deleterious mutations in this gene as a mechanism leading to HSCR. The variants identified in GRB2 are available on request.
To investigate the function of GRB2 in the onset of HSCR, Grb2 expression was downregulated in NLBs from mouse. Grb2 expression was reduced by 30% compared to the control. It is important to mention that this reduction was very similar to GRB2 downregulation found in NLBs from HSCR patients. First, we analyzed the effect of Grb2-knockdown (LV-Grb2) on NLBs' growth, but no changes in size or number of NLBs were identified. Then, we performed an expression analysis of markers associated with the enteric precursor (p75), neural precursor (Nestin) and neuronal precursor cell (βlll-Tubulin). As a result, LV-Grb2 cultures showed a soft decline of Nestin + cells compared to both the uninfected group (control) and non-target control group (LV-off-target) conditions (12% and 14.9% respectively). Data are shown in Figure 2. growth factor receptor (EGFR) signaling pathway. HRAS and TGFB1 also showed functions related to cell cycle and cell migration regulation.

GRB2 Implication in the Onset of HSCR
Based on our results and the fact that GRB2 encodes for the intracellular adaptor protein growth factor receptor-bound 2, which is involved in RET/GDNF/GFRA1 signaling, we decided to further investigate the association of GRB2 with HSCR onset. We performed a mutational screening by direct sequencing of the GRB2 gene, including its coding sequence, intron/exon boundaries and untranslated regions (UTR) in 267 isolated HSCR patients. As a result, we identified several GRB2 rare variants that alone do not seem to significantly contribute to the manifestation of the phenotype, ruling out, therefore, the presence of deleterious mutations in this gene as a mechanism leading to HSCR. The variants identified in GRB2 are available on request.
To investigate the function of GRB2 in the onset of HSCR, Grb2 expression was downregulated in NLBs from mouse. Grb2 expression was reduced by 30% compared to the control. It is important to mention that this reduction was very similar to GRB2 downregulation found in NLBs from HSCR patients. First, we analyzed the effect of Grb2-knockdown (LV-Grb2) on NLBs' growth, but no changes in size or number of NLBs were identified. Then, we performed an expression analysis of markers associated with the enteric precursor (p75), neural precursor (Nestin) and neuronal precursor cell (βlll-Tubulin). As a result, LV-Grb2 cultures showed a soft decline of Nestin + cells compared to both the uninfected group (control) and non-target control group (LV-off-target) conditions (12% and 14.9% respectively). Data are shown in Figure 2.

Discussion
Neural crest-derived cells give rise to ENS through a complex process carefully regulated by interacting signals and a specific gene expression pattern. It is accepted that failures in these processes are responsible for the incomplete gut colonization by these cells during embryonic development, resulting in HSCR. A deeper knowledge about transcriptional programs required for EPCs development, might help us to identify new genes, pathways, and mechanisms involved in HSCR pathogenesis. Specifically, PAX6 was described by our group as a potential susceptibility gene for HSCR, suggesting that it could have a regulatory role on the transcriptional network during the ENS development [12]. Thus, we have focused on the identification of the gene expression pattern established by PAX6 and their potential implication in the ENS development and in HSCR.
In the current study, we identified 11 genes with a potential role during human ENS development based on their expression in human NLBs. ACADM, COL4A2, RABGGTB, ATXN1, and BRD3 were identified as PAX6 targets and HRAS, HTT, MYC, APP, TGFB1, and GRB2 as genes related to such targets. In addition, the approach carried out revealed a group of genes that might be involved in the onset of HSCR, being good candidates to further investigate their role during ENS development (GRB2, BRD3, HRAS, RABGGTB, and TGFB1). The workflow with the selection, prioritization, and further validation of the candidate genes is shown in Figure 3.

Discussion
Neural crest-derived cells give rise to ENS through a complex process carefully regulated by interacting signals and a specific gene expression pattern. It is accepted that failures in these processes are responsible for the incomplete gut colonization by these cells during embryonic development, resulting in HSCR. A deeper knowledge about transcriptional programs required for EPCs development, might help us to identify new genes, pathways, and mechanisms involved in HSCR pathogenesis. Specifically, PAX6 was described by our group as a potential susceptibility gene for HSCR, suggesting that it could have a regulatory role on the transcriptional network during the ENS development [12]. Thus, we have focused on the identification of the gene expression pattern established by PAX6 and their potential implication in the ENS development and in HSCR.
In the current study, we identified 11 genes with a potential role during human ENS development based on their expression in human NLBs. ACADM, COL4A2, RABGGTB, ATXN1, and BRD3 were identified as PAX6 targets and HRAS, HTT, MYC, APP, TGFB1, and GRB2 as genes related to such targets. In addition, the approach carried out revealed a group of genes that might be involved in the onset of HSCR, being good candidates to further investigate their role during ENS development (GRB2, BRD3, HRAS, RABGGTB, and TGFB1). The workflow with the selection, prioritization, and further validation of the candidate genes is shown in Figure 3. Transcription factors bind to short 5-20 bp segments of DNA to regulate gene expression. Such regulatory sequences can be discovered computationally by a computational algorithm that looks for a short, conserved, and often repeated pattern called "the motif". Such motifs are likely to be the transcription factor binding sites [16]. We detected two potential biologically relevant sequences over-represented in the PAX6 binding sequences.
Our results about Edn3 as a PAX6 target allowed us to validate the experimental approach performed in the present study to identify new susceptibility genes for HSCR, as well as to reinforce the role of PAX6 during ENS development. This study has prompted us to examine in detail in future studies the regulatory mechanism of PAX6 over EDN3 since, given that transcription factors can inhibit gene expression [17][18][19], PAX6 might be playing an inhibitory role on EDN3.
The detection of GRB2 as a new candidate gene for HSCR should be mentioned. GRB2 encodes for the intracellular adaptor protein growth factor receptor-bound 2, which is involved in RET/GDNF/GFRA1 signaling. RET activation results in phosphorylation of key docking tyrosines that bind to several intracellular adaptor proteins such as GRB2, and therefore, this protein acts as a Transcription factors bind to short 5-20 bp segments of DNA to regulate gene expression. Such regulatory sequences can be discovered computationally by a computational algorithm that looks for a short, conserved, and often repeated pattern called "the motif". Such motifs are likely to be the transcription factor binding sites [16]. We detected two potential biologically relevant sequences over-represented in the PAX6 binding sequences.
Our results about Edn3 as a PAX6 target allowed us to validate the experimental approach performed in the present study to identify new susceptibility genes for HSCR, as well as to reinforce the role of PAX6 during ENS development. This study has prompted us to examine in detail in future studies the regulatory mechanism of PAX6 over EDN3 since, given that transcription factors can inhibit gene expression [17][18][19], PAX6 might be playing an inhibitory role on EDN3.
The detection of GRB2as a new candidate gene for HSCR should be mentioned. GRB2 encodes for the intracellular adaptor protein growth factor receptor-bound 2, which is involved in RET/GDNF/GFRA1 signaling. RET activation results in phosphorylation of key docking tyrosines that bind to several intracellular adaptor proteins such as GRB2, and therefore, this protein acts as a regulator of this signaling pathway during ENS development [20]. GRB2 can directly bind to RET Y1096 and indirectly to RET Y1062, leading to PI3K/AKT or MAPK signaling activation [21]. A variant at Y1062 caused distal colon aganglionosis reminiscent of HSCR in mice and it has also been described in HSCR patients [22,23]. Based on our data, Grb2 knockdown in ENS precursors resulted in a slight reduction in Nestin-expressing cells (neural precursor cell), suggesting that these enteric precursors could be in a more differentiated state as compared with the control. However, we did not observe any difference in the number of beta-tubulin III-expressing cells, suggesting that Grb2 has no effect on neuronal differentiation. Some studies have described a role of Grb2 that affects the cellular processes related to glial cells [24,25]. Therefore, although further studies are needed, it seems that Grb2 might have a role in glial differentiation.
It is also interesting to notice the identification of APP as a gene with a potential role during human ENS development. According to our study, this gene was previously identified to be highly expressed in enteric NCCs together with Ret and Sox10, as well as a key regulator of ENS development [29], and a potential gene associated with HSCR [30].
In our study, we have observed that these potential susceptibility genes for HSCR were functionally connected to each other. In particular, TGFB1, HRAS, and GRB2 were predicted to be involved in the ERBB signaling pathway, which plays an important role in the development of ENS [31]. Moreover, some genes implicated in this pathway (NRG1 and NRG3) have already been identified as susceptible genes for HSCR [32][33][34]. These functional connections further support the role of these candidate genes in ENS formation.
In conclusion, we deepened the understanding of the expression pattern regulated by PAX6 during the ENS development and determined a set of genes directly regulated by PAX6 through specific regulatory binding sequences. Moreover, we showed a group of genes with a potential implication in ENS development and in HSCR onset. Among them, we identified GRB2 and its direct association with HSCR for the first time. The application of the approach carried out in this study can be an effective method for the knowledge of the genetic background of other rare disorders with a complex genetic basis.

Generation of Enteric NLBs
The isolation of EPCs from postnatal gut represents a powerful tool for the study of ENS development. In culture, these cells form aggregates in suspension known as NLBs, which contain neural crest-derived stem cells and their progeny. Enteric NLBs from human and CD-1 mice (P7) postnatal ganglionic gut were obtained and cultured as previously described by Torroglosa et al. [13]. Primary NLBs without subcultures were used, since secondary and tertiary NLBs acquire different properties compared with primary NLBs [35].

ChIP-seq Assay
ChIP-seq assay was performed in EPCs derived from mouse-NLBs (10 6 ). ChIP-seq and sequencing data analysis were carried out as previously described in detail in Villalba-Benito et al. [14]. A 10% of samples volume/reaction was kept as controls (input) before immunoprecipitation. Mouse IgG (Diagenode) was used as negative control (IgG) for ChIP assay in parallel with Anti-PAX6 antibodies (Millipore, Burlington, MA, USA) at 3 µg. A total of 3 pools of libraries, containing 16 ChIP, 16 inputs, and 4 IgG samples were sequenced in the MiSeq System (Illumina, San Diego, CA, USA). Bowtie was used to align ChIP-seq reads to mouse genome mm9. Model-based Analysis of ChIP-Seq (MACS) was used for calling the genomic regions significantly enriched (peaks). Overlapping peaks from IgG libraries with respect to input were eliminated from the analysis in order to ensure the quality of the results, and to obtain the specific target regions of PAX6. Regarding sequencing data analysis, all the FASTQ files from the ChIPs within the same pool were merged before using the MACS software. Then, BETA minus (http://cistrome.org/ap/root), PAVIS (http://manticore. niehs.nih.gov/pavis2/), and PAPST (https://github.com/paulbible/papst) tools were used to identify the genes located nearest to the regions represented by those peaks. Finally, we selected the genes obtained after the analyses of at least 2 of the 3 software in at least 2 of the 3 pools ( Figure 4).

qRT-PCR
The expression pattern was analyzed in NLB cultures through the TaqMan® Array Plate (Life Technologies, Carlsbad, CA, USA) in the Applied Biosystems 7900HT system (Life Technologies, Carlsbad, CA, USA). Total RNA was isolated using RNeasy Micro kit (Qiagen, Hilden, Germany). Synthesis of cDNA was performed using PrimeScript RT Master Mix (Takara, Kusatsu, Japan). The results were analyzed using the RQ Manager Software (Life Technologies, Carlsbad, CA, USA) based on the comparative cycle threshold (Ct) (ΔΔCt) method. GAPDH was used as an endogenous control gene. Following the software recommendations, we considered expression when Ct values were less

qRT-PCR
The expression pattern was analyzed in NLB cultures through the TaqMan®Array Plate (Life Technologies, Carlsbad, CA, USA) in the Applied Biosystems 7900HT system (Life Technologies, Carlsbad, CA, USA). Total RNA was isolated using RNeasy Micro kit (Qiagen, Hilden, Germany). Synthesis of cDNA was performed using PrimeScript RT Master Mix (Takara, Kusatsu, Japan). The results were analyzed using the RQ Manager Software (Life Technologies, Carlsbad, CA, USA) based on the comparative cycle threshold (Ct) (∆∆Ct) method. GAPDH was used as an endogenous control gene. Following the software recommendations, we considered expression when Ct values were less than 35.

Grb2 Small Hairpin RNA-Expressing Lentiviral Vectors and Infection of NLB Cultures
Lentiviral vectors that constitutively expressed a small hairpin RNA (shRNA) directed against Grb2 mRNA were used to knock-down their expression in mouse NLBs. The ShRNA non-target control was used as the control (Sigma Aldrich, St. Louis, MO, USA). Cultures were divided into 3 groups: Uninfected group (control), shRNA non-target control group (LV-off-target), and infected group (LV-Grb2). Both LV and LV-off-target conditions were infected with a combination of two different shRNAs at a multiplicity of infection of 2. Grb2 expression was evaluated through SYBR Green method (Bio-Rad, Hercules, CA, USA).

Grb2-Knockdown Study in NLB Cultures
In order to evaluate the Grb2-knockdown effect on the neurogenesis and cell growth, the expression of Nestin, βIII-Tubulin, and p75 markers, as well as both the size and number of NLBs were analyzed as previously described in Torroglosa et al. [36].

Sequence Variants Analysis
Potential susceptibility variants for HSCR were searched in the gene sequences with expression in human NLBs using WES data from 56 HSCR patients [15]. Mutational screening of GRB2 was performed in our series of 267 isolated HSCR patients by Sanger sequencing (187 S-HSCR, 27 L-HSCR, 19 TCA, and 55 with phenotype not available) using an ABI PrismH3730 Genetic Analyzer and the SeqScape v2.5 software (Life Technologies, Carlsbad, CA, USA).

Databases and in Silico Tools
Several tools and databases were used in order to perform a comprehensive analysis of the results. The human orthologs of the PAX6 target genes selected in mouse were searched using MGI (http://www.informatics.jax.org). Gene interactions of our target genes as well as their functions were obtained through the IPA tool and GeneMANIA database (http://genemania.org/) [37]. For de novo discovery of motifs in a set of sequences, we used the web application MEME-ChIP (http://meme-suite. org/tools/meme-chip) [38].

Statistical Analyses for Gene Expression Studies
Data were presented as the mean ± SEM (standard error mean) of values obtained from at least three experiments. Comparisons between values obtained in control NLBs and HSCR-NLBs were analyzed using the Student's t-test. Differences were considered significant when the p-value ≤ 0.05.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/21/23/ 9061/s1. Table S1: Additional genes selected because of their interaction with more than one of the PAX6 target genes showed by IPA tool; Figure S1: Average Ct values of candidate genes expressed in NLBs cultures from human. The upper limit of the Ct was set to be 35; Figure S2: Functional networks shared among the potential susceptibility genes for HSCR provided by GeneMANIA database. Funding: This work was supported by Instituto de Salud Carlos III through the projects "PI16/0142" and "PI19/01550" (co-funded by the European Regional Development Fund/European Social Fund, "A way to make Europe"/"Investing in your future") and by the Regional Ministry of Health and Family of the Regional Government of Andalusia through the project PEER-0470-2019. L.V.-B. was supported by a fellowship associated with the CTS-7447 project, which has been funded by the Regional Ministry of Innovation, Science and Enterprise of the Regional Government of Andalusia.