Common genetic variants in GAL, GAP43 and NRSN1 and interaction networks confer susceptibility to Hirschsprung disease

Abstract Hirschsprung disease (HSCR) is a severe multifactorial genetic disorder. Microarray studies indicated GAL,GAP43 and NRSN1 might contribute to the altered risk in HSCR. Thus, we focused on genetic variations in GAL,GAP43 and NRSN1, and the gene‐gene interactions involved in HSCR susceptibility. We recruited a strategy combining case‐control study and MassArray system with interaction network analysis. For GAL,GAP43 and NRSN1, a total of 18 polymorphisms were assessed in 104 subjects with sporadic HSCR and 151 controls of Han Chinese origin. We found statistically significant differences between HSCR and control groups at 5 genetic variants. For each gene, the haplotypes combining all polymorphisms were the most significant. Based on SNPsyn, MDR and GeneMANIA analyses, we observed significant gene‐gene interactions among GAL,GAP43,NRSN1 and our previous identified RELN,GABRG2 and PTCH1. Our study for the first time indicates that genetic variants within GAL,GAP43 and NRSN1 and related gene‐gene interaction networks might be involved in the altered susceptibility to HSCR in the Han Chinese population, which might shed more light on HSCR pathogenesis.


| INTRODUCTION
Hirschsprung disease (HSCR) is a complex genetic disorder caused by congenital defect of the enteric nervous system (ENS) which is derived from neural crest cells (NCCs). HSCR affects approximately 1/5000 live births worldwide, and the highest incidence was observed in Asian population (2.8/10 000 live births). 1 Based on the extent of aganglionosis, the HSCR cases can be anatomically categorized into three subtypes: short segment HSCR (S-HSCR, 80% of cases) in which the aganglionic segment does not extend beyond the upper sigmoid, long segment HSCR (L-HSCR, 15% of cases) and total colonic aganglionosis (TCA, 5% of cases). 2 Importantly, HSCR shows a dramatic sex bias with at least 4 times more males affected than females in S-HSCR (male:female % 1:1 in L-HSCR) for causes that remain unclear. 3 Genetic factors or multiple gene interactions are crucial to the development of Hirschsprung disease as HSCR is a non-Mendelian disorder in nature with low sex-dependent penetrance and interfamilial variation. 4 It has been suggested that there are variations in penetrance and severity of aganglionosis between family members Yang Wang and Weihui Yan are contributed equally to this work. bearing mutations in HSCR genes. 5 Obviously, only a single homozygous null mutation of HSCR-related genes is insufficient to cause serious aganglionosis phenotype in HSCR. 3 Genetic variants of at least 15 genes so far have been implicated in HSCR aetiology, including RET (receptor tyrosine kinase), one of the major HSCR susceptibility genes, 6,7 whereas only~0.1% of the heritability in HSCR can be attributed to the mutations in these genes that account for about 50% of familial and 7%-35% of sporadic HSCR cases, 8 indicating more genes that might be involved in HSCR development.
On the other hand, recent genomewide association studies have revealed dozens of novel HSCR genes, which may facilitate the description of a complete landscape of genetic networks in HSCR.
Taking advantage of whole exome sequencing, several genes, including DENND3, FAT3 and AGL, were linked to HSCR pathogenesis. [9][10][11] NRG3 has recently been proved to be a new HSCR risk gene based on exome sequencing and genomewide copy number analysis, 2,12 which was further confirmed by our previous work. 13 In addition, genomewide association studies on HSCR trios and sporadic cases have uncovered the class 3 semaphorin gene cluster and certain large-scale chromosomal aberrations regarding HSCR aetiology. 14,15 Recent genomewide microarray analysis has reported the levels of GAL (galanin), GAP43 (growth-associated protein 43) and NRSN1 (neurensin 1) were significantly down-regulated in HSCR cases when compared to controls, indicating the possibility that all 3 genes might be associated with HSCR risk. 16 More importantly, joint gene-gene effects, such as RET and PHOX2B genes, might have a crucial impact on the development of HSCR. 17 Our previous study has proved the interactions among GABRG2, RELN and PTCH1 may contribute to altered susceptibility to HSCR. 13 Additionally, galanin-expressing GABA neurons in the lateral hypothalamus may have important implications for treatment strategies of psychiatric disorders. 18 In Ptch1 (+/À) mice that causes aberrant hedgehog signalling, reduced Gap43 expression leads to the Nos2-mediated medulloblastoma development. 19 Recently, it has been demonstrated that reelin blockade results in decreased levels of phospho-GAP43 in the superior colliculus, suggesting the interaction of reelin signalling and phospho-GAP43 might be involved in the development of neural circuits. 20 With all these lines of evidence and results, we aimed to explore whether genetic variants within GAL, GAP43 and NRSN1 might contribute to the altered susceptibility to HSCR, and based on the 18 polymorphisms involved in this study ( Figure 1A), we further assessed the interaction relationship among GAL, GAP43, NRSN1 and our previous identified GABRG2, RELN and PTCH1 genes.

| Subjects
The subject group involved in this study consisted of 104 cases with HSCR (84 male and 20 female) and 151 normal controls (86 male and 65 female). The mean ages of HSCR group and control group were 1.14 AE 1.83 years and 1.66 AE 1.05 years. The characteristics of the study subjects can be found in our previous study. 21 All the participants in the study were of Han Chinese origin and were recruited from the residents who were biologically unrelated.
F I G U R E 1 Distribution and representative mass spectra of the genetic variants in the present study. A-C, The 18 genetic variants distributed in GAL, GAP43 and NRSN1. Red lines indicate the studied SNPs; blue lines and arrows represent the exons located in theGAL, GAP43 and NRSN1; D, Representative mass spectra of the 8 polymorphisms in GAP43. Blue dotted lines indicate the presence of the studied alleles; red dotted lines represent no allele detected; grey dotted lines denote the unrelated peaks Diagnosis of HSCR was confirmed by the histological examination of either surgical resection material or biopsy for the absence of ganglion cells. The HSCR group included 86 subjects of S-HSCR (short segment HSCR), 15 subjects of L-HSCR (long segment HSCR) and 3 subjects of TCA (total colonic aganglionosis). Controls were randomly enrolled from the subjects with no history of chronic constipation.
The protocol of our study was reviewed and approved by the ethics committee of Xin Hua Hospital. Informed consent was obtained from parents of all participants after the procedure had been fully explained. All experiments were conducted in accordance with the tenets of the Declaration of Helsinki. DNA extraction was performed according to standard procedures with QIAamp DNA blood midi kit (Qiagen, Valencia, CA).

| Genotyping and quality control
Genotyping was carried out using the MassARRAY iPLEX Gold technology (Sequenom, San Diego, CA). Briefly, PCR and iPLEX singlebase extension primers (SBE) were designed taking advantage of the Assay Design Suite of Sequenom. The whole process consisted of the PCR amplification, the shrimp alkaline phosphatase (SAP) and the primer extension reactions using iPLEX Gold assay (Sequenom) that discriminates sequence differences at the single nucleotide level.
Mass signals for the different alleles were captured by MALDI-TOFbased system with high accuracy. Raw data from the assays were processed with Typer Version 4.0 (Sequenom).
We recruited the following criteria as a measure of acceptable genotyping: (1) 30 sample duplicates and 4 blank wells were involved in each 384-well plate; (2) concordance rate for the duplicates ≥ 99.5%; (3) call rate for the blank wells <5% in each 384-well plate; (4) call rate > 95% for each 384-well plate; and (5) overall call rate by individual or by marker > 95%. The data for any marker or individual failing the criteria were excluded from further analyses.

| SNP-SNP interaction analysis
In this study, SNPsyn (http://snpsyn.biolab.si) 22 was employed to interrogate the SNP-SNP interaction networks regarding HSCR. The genotyping data of the studied genetic variants were processed with the SNPsyn software tool, using which we uncovered the SNP-SNP interaction networks and carried out the SNP pair selection that was mainly based on information gain (I), synergy (Syn) and false discovery rate (FDR). 22,23 Moreover, the multifactor dimensionality reduction (MDR) analysis was included to explore the gene-gene interactions. We used the MDR software version 3.0.2 to perform the MDR analysis and identified all risk factors in the best model maximizing testing accuracy and cross-validation consistency (CVC). 24 We further recruited the GeneMANIA database, including co-expression, co-localization and genetic interaction datasets, to assess the gene-gene interaction networks and to conduct a function prediction. 25

| Statistical analysis
We used SHEsis (http://analysis.bio-x.cn/myAnalysis.php) to calculate Hardy-Weinberg equilibrium, allelic and genotypic association, odds ratio (OR) and 95% confidence interval (CI) and to estimate allelic distribution and linkage disequilibrium (LD). 26 "D" was included as the standardized measure for all possible pairs of SNP loci. All the P values in this study were two-tailed, and the significance level was set at P = .05. Bonferroni correction was performed to correct the P values of genetic analysis, and Plink was enrolled to conduct the association analyses with dominant model and recessive model, and perform the adjustment for gender factor in the association analysis. 27 Additionally, haplotype distribution was estimated using the program UNPHASED, 28 and power calculations were conducted using the G*Power 3 program. 29

| RESULTS
In regard to the studied genetic variants, Hardy-Weinberg equilibrium tests were conducted in HSCR group and control group, respectively. Allele and genotype frequencies of the 18 markers are listed in Tables 1-3 and (5) NRSN1_rs3829810, Dom P = .013, Rec P = .049. PLINK was recruited in the adjustment for gender factor, and the findings in the 5 positive SNPs remained significant after correction. Figure 1B presents representative mass spectra of the original MassARRAY reactions in GAP43. Additionally, the frequencies of certain alleles and genotypes regarding the 5 positive markers were significantly higher in HSCR group compared to normal control group, such as the A allele and AA genotype of GAL rs1042577, the T allele and TT genotype of GAP43 rs283367, the G allele and GG genotype of GAP43 rs14360, the G allele and GG genotype of NRSN1 rs10946675, and the C allele and CC genotype of NRSN1 rs3829810.
We then performed LD and haplotype analyses of genetic variants in the 3 genes as haplotypes constructed from polymorphisms with strong LD will increase the statistical power for association with the disease. Figure S1 shows LD for each pair of SNPs in HSCR group and control group. Strong LD was observed in the following marker groups: (1) GAL, rs1546309-s1042577; (2) GAP43, rs2118604-rs12632276; and (3) NRSN1, rs4449613-rs6935378, rs4449613-rs3178 and rs6935378-rs3178. We thus interrogated the haplotype distributions for these markers in the later analysis.
We selected haplotypes with strong LD for presentation (Table S1). As there were significant frequency discrepancies between HSCR and control groups, several haplotypes were observed to be strongly associated with HSCR. Additionally, haplotype analysis of these 18 polymorphisms revealed some significant global P values (Table S2). For each gene, the haplotypes that combined all markers were the most significant (GAL, P = 6.78 9 10 À8 ; GAP43, P = 4.16 9 10 À12 ; NRSN1, P = .0095). We further included G*Power 3 program in the power calculations and found our sample size had >80% power to detect a significant association (P < .05) for alleles, genotypes and haplotypes when an effect size index of 0.24 (corresponding to a "weak" gene effect) was adopted. We further compared the SNP frequency of normal controls in our present study with the SNP frequency of CHB (Han Chinese in Beijing, China) in 1000 Genomes Project Phase3 database (http://asia.ense mbl.org), and no significant difference was observed between these 2 datasets (Table S3).
Moreover, SNPsyn software tool was enrolled to interrogate the SNP-SNP interactions among GAL, GAP43, NRSN1 and our previous studied GABRG2, RELN and PTCH1 genes. 13 We investigated the SNP-SNP interaction networks based on both information gain (I) and synergy (Syn), and recruited only the SNP pairs with significant scores (I, Syn) in the network analysis ( Figure 2). 22 In our study, significant scores were found at several SNP pairs, corresponding to GAL-GAP43, GAL-NRSN1, PTCH1-GABRG2-GAP43 group, etc.
( Figure 2). The positive SNPs associated with HSCR were also involved in the SNP-SNP interactions, such as GAL_rs1042577, We further employed the MDR strategy to explore the potential gene-gene interactions among GAL, GAP43, NRSN1, GABRG2, RELN and PTCH1 corresponding to the best interaction model ( Figure 3A-C, Table 4). As for HSCR risk prediction, the best single factor model was GAP43 (rs14360) (testing accuracy = 0.5813; CVC = 10/ 10), which was significantly associated with HSCR. GAL To interrogate the functional association networks among these 6 HSCR-related genes, we included the GeneMANIA online software in the present study using the parameters limited to co-expression, co-localization and genetic interactions ( Figure 3D). The 6 genes interacted with each other mainly through co-expression and genetic interactions, and only the interaction between GAP43 and GABRG2 was partly due to co-localization. Moreover, gene function prediction showed NRSN1, RELN and GABRG2 might contribute to neuron projection, neuron-neuron synaptic transmission and neuron part.  The best model was referred to as the one with the maximum testing accuracy and maximum cross-validation consistency (CVC).

| DISCUSSION
HSCR is a congenital intestinal obstruction characterized by a deficit in migration of enteric neural crest cells (ENCCs), or by a defect in proliferation, differentiation or survival of ENCCs once they reach the intestinal tract. 5 As a model of non-Mendelian genetic disorder, HSCR can be attributed to multiple gene-gene interactions that modulate the ability of ENCCs to populate the developing gut, and therefore, the synergistic effects of multiple hypomorphic mutations in HSCR-related genes could affect disease penetrance and severity. 3 However, a complete landscape of genetic networks in HSCR remains obscure. Our present study provided first evidence that genetic variants within GAL, GAP43 and NRSN1 might contribute to the altered susceptibility to HSCR, and the interaction networks among GAL, GAP43, NRSN1 and our previous identified GABRG2, RELN and PTCH1 genes might confer an increased risk in HSCR.
Our findings suggested a significant association of GAL (rs1042577) with the altered susceptibility to HSCR. As rs1042577 is located in the untranslated region, this genetic variant may exert an impact on the regulatory mechanisms of gene expression. 30 Moreover, we observed that the G allele and GG genotype of rs1042577 were less frequent in HSCR group compared to normal control group, which indicated that the G allele and GG genotype might be involved in a protective effect against HSCR, and yet the A allele and AA genotype of rs1042577 were more common in HSCR cases than in controls, implying that all might be the risk factors for HSCR. Galanin encoded by GAL is a neuroendocrine peptide, which is widely expressed in the central and peripheral nervous systems and also the gastrointestinal tract. 31 Additionally, galanin modulates transmitter release from myenteric neurons via inhibition of voltagedependent calcium channels mediated by G-protein-coupled receptors. 32 A recent genomewide study 16 suggested GAL as a candidate for HSCR due to the reduced level of GAL expression in HSCR group compared with control group, and our results further supported this opinion.
We further interrogated the association between GAP43 and HSCR, and found that 2 genetic markers (rs283367 and rs14360) within GAP43 gene presented a strong association with the HSCR risk. Saeed et al 16  On the other hand, we tried to assess the relationship between NRSN1 gene and HSCR risk. As rs10946675 and rs3829810, the 2 positive SNPs found within NRSN1, were located in the intronic and untranslated regions, respectively, these polymorphisms might be involved in the regulatory mechanisms of gene expression. Specifically, we noticed the G allele and GG genotype of rs10946675 and the C allele and CC genotype of rs3829810 were more frequent in HSCR group than in control group, indicating that all may contribute to the altered risk of HSCR. Neurensin 1 (NRSN1) is a neuron-specific protein comprising one microtubule-binding domain and several membrane domains, and it is particularly abundant in neuronal processes, such as the process of neurite extension. 35 NRSN1, as a key regulator, may function in neuronal organelle transport and in the conduction of nerve signals, therefore contributing to axonal regeneration and development. 36 In addition, the expression of NRSN1 has been proved to be down-regulated in HSCR cases compared with normal controls, 16 further supporting NRSN1 gene as a potential susceptibility gene to HSCR.
In the present study, the significance regarding haplotypes might contribute to the altered risk of HSCR as well (Table S1), as under certain conditions haplotype analysis may increase the power to detect disease loci compared with the single SNP analysis. 37 As for each of the 3 genes, the most significant haplotype involved all genetic variants in the corresponding gene (Table S1). Interestingly, certain significant haplotypes might be the protective factors in HSCR, such as GAL_T-C-T-G (rs1546309-rs3136540-rs3136541-rs1042577, P = 1.62 9 10 À5 , OR = 0.43, 95% CI 0.29-0.64).
As the interactions among GAL, GABA signalling, GAP43, PTCH1 and reelin might play a crucial role in the neural functions and related disease processes, [18][19][20] we utilized the SNPsyn platform to further interrogate the SNP-SNP interaction networks among GAL, GAP43, NRSN1 and our previous identified GABRG2, RELN and PTCH1 genes. On the other hand, synergistic combinations carry more information compared to the sum of information contained in individual SNPs and specifically may carry information in regard to the phenotypes. 38 In the present study, we found significant interaction networks in the GAL-GAP43-NRSN1 and GAL-GAP43-NRSN1-GABRG2-RELN-PTCH1 group, respectively ( Figure 2). All the 5 positive genetic variants observed within GAL, GAP43 and NRSN1 were involved in the interaction networks ( Figure 2C,D). Moreover, we noticed that certain SNP pairs involved in the networks were located within the same gene, such as rs283367-rs14360 (GAP43), indicating that the cis-regulation effect might facilitate this kind of SNP-SNP interaction. 39 We further recruited multifactor dimensionality reduction (MDR) method to evaluate the gene-gene interactions on the risk of HSCR using the data in regard to the 6 HSCR-associated genes. MDR was a nonparametric approach that does not require specification of a genetic model to detect gene-gene interactions without main gene effects. 40 Taking advantage of the MDR analysis, we assessed the best interaction model with the maximum testing accuracy and maximum CVC between all the genes involved in our study. Of note, all the positive SNPs within GAL, GAP43 and NRSN1 were included in the best models obtained from the MDR analysis (Table 4). Specifically, the best two-factor model, GAL (rs1042577)-PTCH1 (rs28485160), was also identified in the SNPsyn analysis ( Figure 2).
Compared with other models, the best four-factor model, GAP43 (rs14360)-GAP43 (rs283367)-NRSN1 (rs3829810)-PTCH1 (rs28701981), presented the most significant OR, raising the possibility that a multifactor model was more likely to facilitate the increased risk to HSCR.
By utilizing the GeneMANIA approach, we further explore the functional networks between these 6 genes involved in HSCR risk ( Figure 3D). These 6 genes functionally connected to each other via co-expression, co-localization and genetic interactions, and in partic-