Aberrant Expressions and Variant Screening of SEMA3D in Indonesian Hirschsprung Patients

Background: The semaphorin 3D (SEMA3D) gene has been implicated in the pathogenesis of Hirschsprung disease (HSCR), a complex genetic disorder characterized by the loss of ganglion cells in varying lengths of gastrointestinal tract. We wished to investigate the role of SEMA3D variants, both rare and common variants, as well as its mRNA expression in Indonesian HSCR patients. Methods: Sanger sequencing was performed in 54 HSCR patients to find a pathogenic variant in SEMA3D. Next, we determined SEMA3D expression in 18 HSCR patients and 13 anorectal malformation colons as controls by quantitative real-time polymerase chain reaction (qPCR). Results: No rare variant was found in the SEMA3D gene, except one common variant in exon 17, p.Lys701Gln (rs7800072). The risk allele (C) frequency at rs7800072 among HSCR patients (23%) was similar to those reported for the 1,000 Genomes (27%) and ExAC (28%) East Asian ancestry controls (p = 0.49 and 0.41, respectively). A significant difference in SEMA3D expression was observed between groups (p = 0.04). Furthermore, qPCR revealed that SEMA3D expression was strongly up-regulated (5.5-fold) in the ganglionic colon of HSCR patients compared to control colon (ΔCT 10.8 ± 2.1 vs. 13.3 ± 3.9; p = 0.025). Conclusions: We report the first study of aberrant SEMA3D expressions in HSCR patients and suggest further understanding into the contribution of aberrant SEMA3D expression in the development of HSCR. In addition, this study is the first comprehensive analysis of SEMA3D variants in the Asian ancestry.


INTRODUCTION
Hirschsprung disease (HSCR: MIM# 142623) is a complex genetic disorder, characterized by the lack of ganglion cells in the bowel, resulting in a functional obstruction during infancy (1). HSCR is categorized into the following types: short-segment, long-segment, and total colonic aganglionosis (1,2). HSCR incidence differs among ethnic groups, with 1.5, 2.1, and 2.8 cases per 10,000 live births in European, African, and Asian ancestry cases, respectively (1,2). There are at least 17 genes responsible for the development of HSCR, with most of them being members of the RET and EDNRB signaling pathways (1,2). Two genetic risk factors are the RET rs2435357 and rs2506030 variants (3,4). Our recent studies showed that the RET rs2435357 and rs2506030 risk alleles have higher frequency in Indonesian ancestry cases as compared with European ancestry cases (5,6), which might relate to the higher incidence of HSCR in Indonesia (3.1 cases per 10,000 live births) than other populations (7).
The third signaling pathway of HSCR pathogenesis includes class 3 semaphorins (SEMA3s), involving SEMA3D (4,8,9) SEMA3D has been implicated in the development of HSCR and contributes to risk through both common and rare variants in European ancestries (4,8,9), as evidenced by (1) the detection of both common and rare SEMA3D variants in HSCR patients; (2) the expression of SEMA3D in the human, mouse, and zebrafish intestines and, particularly, the enteric nervous system (ENS); and (3) the joint effect of Ret and Sema3d loss of function in an aganglionosis animal model. However, our recent study showed that the effect of SEMA3 rs11766001 common variant on HSCR depends on the ethnic background (10). In addition, the allele frequencies of common variants might differ among Asians, since the North Asians, Han Chinese, Japanese, and Southeast Asians can be distinguished based on their Y chromosome variants (11). Moreover, alterations in the expression of specific genes have been implicated in the development of HSCR (12)(13)(14)(15). Therefore, we wished to investigate the role of SEMA3D variants, both rare and common variants, as well as its mRNA expression in Indonesian HSCR patients.
All parents signed a written informed consent form before participating in this study. The Institutional Review Board of the Faculty of Medicine, Public Health, and Nursing, Universitas Gadjah Mada/Dr. Sardjito Hospital gave approval for this study (KE/FK/1356/EC/2015). All experiments were performed in accordance with relevant guidelines and regulations.

Polymerase Chain Reaction (PCR) and DNA Sequencing
A QIAamp DNA Extraction Kit (QIAGEN, Hilden, Germany) was used to extract genomic DNA from whole blood from each individual, according to the manufacturer's instructions. We stored the extracted DNA samples at −20 • C until analysis. PCR was conducted using a Swift Maxi thermal cycler (Esco Micro Pte. Ltd., Singapore), followed by Sanger sequencing analysis to

DNA Genotyping
DNA genotyping was performed using Sanger sequencing analysis. The SEMA3D rs7800072:A>C (chr7: g. 84,628,989A>C) variant was identified during the Sanger sequencing analysis to find a rare variant in Indonesian HSCR patients. The risk allele (C) was determined according to the 1,000 Genomes Project and ExAC population databases (17,18).

RNA Extraction and Quantitative Real-Time PCR (qPCR)
The ganglionic and aganglionic intestinal specimens were collected at pull-through surgeries from 18 HSCR patients, while control intestinal specimens were obtained at colostomy closure from 13 anorectal malformation (ARM) patients ( Table 2).
Total RNA was isolated from colonic specimens using the total RNA Mini Kit (Tissue) (Geneaid Biotech Ltd., New Taipei City, Taiwan). RNA was quantified by a NanoDrop 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The OD260/280 ratios ranged from 1.8 to 2.0, indicating high RNA purity.
SEMA3D gene expression was quantified using the BioRad CFX Real-Time PCR System (California, USA) and the SensiFAST TM SYBR R No-ROX One-Step Kit (Bioline, Meridian Bioscience, Memphis, USA) using the 5 ′ -CAACGCAGCCTG ATAAACAA-3 ′ (forward) and 5 ′ -TCTTTCATCTCTTGTGGGG AGT-3 ′ (reverse) (19) The primers were designed to bridge SEMA3D exon 8 and 9 junctions (19). Glyceraldehyde-3phosphate dehydrogenase (GAPDH), a housekeeping gene, was used as an endogenous control. The GAPDH primers were 5 ′ -GCACCGTCAAGGCTGAGAAC-3 ′ (forward) and 5 ′ -TGGTGAAGACGCCAGTGGA-3 ′ (reverse). qPCR reactions contained SensiFAST TM SYBR R No-ROX One-Step mix (2×) 10 µL, RiboSafe RNase Inhibitor 0.4 µL, reverse transcriptase 0.2 µL, forward primer (10 µM) 0.8 µL, reverse primer (10 µM) 0.8 µL, and total RNA 50 ng, with final volume of 20 µL. qPCR was performed for 10 minutes (min) at 45 • C for reverse transcription process, followed by 2 min at 95 • C and 39 cycles for 5 s at 95 • C, 10 s at 58 • C and, 5 s at 72 • C, and 1 cycle for 5 s at 65 • C, according to the manufacturer's instructions. We performed the gel electrophoresis for the qPCR of SEMA3D and GAPDH (Supplemental Figure 1). The Livak method was utilized to determine the SEMA3D mRNA expression level (20). This method is designed to calculate a relative gene expression and referred to as the C T method. The (log) expression is proportional to the negative C T value (the lower the C T , the higher the expression) (20).

Statistical Analysis
SEMA3D expression was described as a mean value ± SD. The Kolmogorov-Smirnov test was used to determine the data distribution, and a one-way ANOVA was utilized to assess statistical differences between groups. A chi-square test was used to establish p-value for the case-control association analysis for SEMA3D rs7800072 variant. A p < 0.05 was considered significant.

RESULTS
Most of our HSCR patients were male (70%) and short-segment aganglionosis type (98%), with the mean age at diagnosis being 34.6 ± 44.5 months ( Table 1). Among the 54 HSCR patients, 49, 3, and 2 children underwent a definitive surgery, a colostomy, and a full-thickness rectal biopsy waiting a pull-through procedure, respectively, with the most common definitive surgery conducted being transanal endorectal pull-through (43%) ( Table 1); and our controls were six males and seven females, with a mean age during stoma closure of 47 ± 45.1 months.

DISCUSSION
In this study, we have performed an in-depth genetic and gene expression study of the SEMA3D in Indonesian HSCR patients. We did not detect any rare variants in this gene, although previous studies in European ancestry cases have identified rare coding variants in SEMA3D associated with HSCR. Our results indicate that the association of such variants to the disease may be restricted to specific ethnic groups. It may still be the case that other semaphorin 3 genes could play a role in HSCR pathogenesis (4,9). A previous study showed the existence of other genetic factors conferring risk to HSCR in specific ethnic populations (3). For example, there were two different RET haplotypes involving the enhancer mutation that were overtransmitted to the HSCR offspring in the European sample, while in the Chinese sample, only one of those haplotypes was present (3). It might be speculated that the enhancer mutation arose in one haplotype, which after the Asian-Caucasian split, rearranged to also give the other haplotype, but exclusively in  the Caucasian (22). Another example is the NRG1 rs7835688 genetic marker, which has been originally discovered in Chinese HSCR patients (23) and has been shown in other Asian ancestry cases (5,24), but is rare and shows no effect in the Caucasian population (8,22). Makhmudi et al. (25) also demonstrated that the MTHFR c.677C>T is a genetic risk factor for Indonesian gastroschisis, but not seen in Caucasians (26,27). A recent study showed that loss of Sema3d in null homozygotes mice had no impact on the intestinal transcriptome (28). The authors were unable to find evidence for Ret and Sema3d interaction affecting survival, presence of myenteric plexus, or intestine transcriptome (28). Furthermore, Tang et al. revealed that the effects of RET and NRG1 variants are universal across Caucasian and Asian ancestries, but the impact of SEMA3 variant was restricted to Caucasian ancestries (29). It should be noted that our screening method did not cover the promoter and enrich regions of CpG of SEMA3D.
In addition, the observed SEMA3D rs7800072 variant frequency in Indonesian HSCR patients is similar to those reported for the 1,000 Genomes Project and ExAC East Asian ancestry controls (17,18). Therefore, we might conclude that this common variant does not have a role on the development of HSCR in Indonesia. ENS development is a complex process, regulated by a large range of molecules and signaling pathways. This strictly controlled process needs the correct regulation of ENS-specific gene expression. The expressions of several genes implicated in HSCR development have been shown to be regulated by epigenetic mechanisms. RET expression was regulated by retinoic acid through DNA methylation on this 5 ′ -CG-3 ′rich enhancer region (30), while a significantly lower level of EDNRB methylation was also detected in HSCR patients (31). Interestingly, our study clearly demonstrates that SEMA3D gene expression was strongly up-regulated in the ganglionic intestines of HSCR patients as compared to controls. To the best of our knowledge, our study is the first report of the aberrant CI, confidence interval; C T , cycle threshold; SD, standard deviation; *p < 0.05 is considered statistically significant. expressions of SEMA3D in HSCR patients. Another novelty in our findings is we conducted the study on patients with Indonesian ancestry [vs. European (4) ancestry cases]. It has been shown that reduced sema3d expression by morpholino resulted in severe effects on the intestine and its innervation of zebrafish (4). In addition, sema3d bound to nrp1a to facilitate the axonal guidance and contributed to peripheral axon outgrowth interdependently with dpysl3 (32, 33) Binding of Sema3d to neuropilin and plexin receptors introduced biochemical responses in specific neurons and stimulated the neuron migration (34,35). Therefore, we might hypothesize that the aberrant expressions of SEMA3D will have an impact in our HSCR patients by affecting the neuronal guidance during ENS development.
It has been shown that some HSCR patients have persistent bowel symptoms, such as constipation, soiling, and enterocolitis, after an appropriate pull-through procedure. Most HSCR patients with persistent bowel symptoms do not have any identifiable etiology for their ongoing bowel dysmotility (36). The current hypothesis is that the aberrant expression of some genes in the ganglionic colon of HSCR patients includes SK3, Cx26, ChAT, and nNOS (37)(38)(39)(40). Our study presented the altered SEMA3D expressions in the ganglionic colon of HSCR patients. Therefore, we might hypothesize that the aberrant SEMA3D expressions in the ganglionic colon involve in the pathogenesis of persistent bowel symptoms in HSCR patients following a properly performed pull-through surgery.
Our results should be interpreted with some caution, however, because they are based on overall means without accounting for other factors generating variation in the data such as gender, age, and degree of aganglionosis. Further studies of the methylation pattern of the SEMA3D gene are required to investigate whether the aberrant expression of SEMA3D in HSCR patients is due to abnormal DNA methylation. It is also necessary to compare the SEMA3D protein expression level between HSCR patients and controls and to identify its location in the colon tissue to prove the increased SEMA3D expressions. Unfortunately, we do not have any data on the methylation pattern, protein expression, and immunohistochemistry of SEMA3D due to resource limitations in our institution. Furthermore, our study utilized ARM patients' colon as controls. It has been shown that most ARM patients have abnormal colonic motility (41) and low expression of interstitial cells of Cajal marker (42). Therefore, further research with more proper control colon (e.g., autopsy bowel material from healthy infants or trauma patients) is necessary to better confirm the role of SEMA3D expression in the pathogenesis of HSCR. Noteworthy, however, our small sample size is a limitation of our study and suggests that a multicenter larger sample population needs to be studied to clarify our findings.
In conclusion, we report the first study of aberrant SEMA3D expressions in HSCR patients and suggest further understanding into the contribution of aberrant SEMA3D expression in the development of HSCR. In addition, this study is the first comprehensive analysis of SEMA3D variants in the Asian ancestry.

DATA AVAILABILITY STATEMENT
All data generated or analyzed during this study are included in the submission. The raw data are available from the corresponding author on reasonable request.

ETHICS STATEMENT
The Institutional Review Board of the Faculty of Medicine, Public Health and Nursing, Universitas Gadjah Mada/Dr. Sardjito Hospital gave approval for this study (KE/FK/1356/EC/2015). All parents signed a written informed consent form before participating in this study.

AUTHOR CONTRIBUTIONS
G and KI conceived the study. G drafted the manuscript. KI and JV critically revised the manuscript for important intellectual content. AK, NB, and FR performed total RNA extraction and qPCR. HH, AM, MF, and DY conducted the experimental PCRbased work for Sanger sequencing. All authors have read and approved the manuscript and agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

ACKNOWLEDGMENTS
We are grateful to Program Bantuan Seminar Luar Negeri, Ditjen Penguatan Riset dan Pengembangan, Kemenristekdikti, Indonesia for providing the conference grant to present the abstract at the 52nd Annual Meeting of PAPS in Christchurch, New Zealand at March 10-14, 2019. We would like to thank the patients, their families, and all those who were involved and provided excellent technical support and assistance during the study. We also thank Professor Aravinda Chakravarti (New York University, School of Medicine) for critically reading the manuscript and his suggestions.