https://orcid.org/0000-0002-5076-1698
Abstract
Disruptions in gene regulatory networks (GRNs), driven by multiple deleterious variants, potentially underlie complex traits and diseases. Hirschsprung disease (HSCR), a multifactorial disorder of enteric nervous system (ENS) development, is associated with at least 24 genes and seven chromosomal loci, with RET and EDNRB as its major genes. We previously demonstrated that RET transcription in the ENS is controlled by an extensive GRN involving the transcription factors (TFs) RARB, GATA2 and SOX10 and other HSCR genes. We now demonstrate, using human and mouse cellular and animal models, that EDNRB is transcriptionally regulated in the ENS by GATA2, SOX10 and NKX2.5 TFs. Significantly, RET and EDNRB expression is regulated by their shared use of GATA2 and SOX10, and in turn, these TFs are controlled by EDNRB and RET in a dose-dependent manner. This study expands the ENS development GRN to include both RET and EDNRB, uncovers the mechanistic basis for RET–EDNRB epistasis and emphasizes how functionally different genes associated with a complex disorder can be united through a common GRN.
Introduction
Although complex traits and disorders have now been shown to comprise small genetic effects at many loci, by genome-wide association (GWAS) and sequencing studies, there is no unifying hypothesis explaining why a particular constellation of genes affects a specific trait (1,2). Some studies attempt to search for enrichment of different functions or pathways among the candidate genes at GWAS loci, but these analyses explain only a minority of the genes or show enrichment for very general protein functions (3,4). Boyle and colleagues have suggested instead that functional networks are shallow, each gene being transcriptionally connected to many genes and being broadly expressed across many cell types (4). Thus, GWAS can fail to distinguish between “core” genes affecting the biology of a trait and “peripheral” genes whose activity are altered as a collateral consequence. A significant task is, therefore, to uncover the core gene regulatory network (GRN) underlying a trait or disease (2). GRNs are modular, comprised of a small number of subcircuit classes and are conserved (5,6), and as we have previously demonstrated, the majority of the genes in core GRNs can be enriched in human disease (7).
In a recent study (7), we deciphered the core enteric GRN underlying Hirschsprung disease (HSCR; congenital intestinal aganglionosis), which is a severe, common (~15/100 000 live births) developmental disorder of the enteric nervous system (ENS) in which the gut fails to be innervated. Genetic studies have identified rare coding variants in 24 genes along with common enhancer variation in two genes and large copy number changes at seven loci to explain the genetic basis of 72% of HSCR cases (8). The features common to these HSCR genes are that they are expressed in the ENS and appear to regulate early differentiation of enteric neural crest cells (ENCCs) into enteric neurons and their subsequent proliferation and migration along the gut. Among these, the two most HSCR genes contributing to the greatest susceptibility encode RET, a receptor tyrosine kinase, and EDNRB, a G protein coupled receptor (9). These two genes are epistatic, in human patients as well as in mouse models of aganglionosis (10–12). The molecular mechanism of this interaction is unknown but important to decipher for understanding the role of early ENS development in HSCR and how EDNRB is related to the RET GRN.
Genetic control of gene expression is mediated through multiple cis-regulatory elements (CREs) and enhancers, which bind specific transcription factors (TFs) organized within a topologically associating domain (TAD) (13,14). Starting with a core TAD and gut-specific epigenomic marks, we have previously deciphered the RET GRN and demonstrated that disruption of RET enhancer function leads to variable loss of RET function and disruption of its enteric GRN (7). Given the strong epistasis between RET and EDNRB, we, therefore, asked whether genetic control of EDNRB shared components with the RET GRN. Here, we first identify five enhancers for enteric EDNRB gene expression and show that two of these bind GATA2 and NKX2.5; further, the EDNRB promoter binds SOX10. Second, by using cellular and animal models, we demonstrate that GATA2 and SOX10 are subject to transcriptional feedback by EDNRB in a dosage-sensitive manner. Thus, common transcriptional control of RET and EDNRB by GATA2 and SOX10 and common feedback on these two TFs is the most parsimonious explanation of RET–EDNRB epistasis. Shared genetic control is, however, insufficient for RET and EDNRB gene expression because these genes are uniquely controlled by RARB and NKX2.5, respectively (7). This study demonstrates how two major genes in ENS development are functionally regulated via a single GRN and how dose-dependent interactions within this GRN are fundamental to HSCR susceptibility.
Results
EDNRB gut enhancers
To identify cis-regulatory elements, we searched for DNase I hypersensitive (DHS), H3K4me1 and H3K27ac sites within a 341 kb (Chr13:78359528-78 701 413) TAD containing EDNRB and common to seven human cell lines (GM12878, HMEC, HUVEC, IMR90, K562, KBM7 and NHEK) (14). Candidate EDNRB human fetal gut enhancers were identified at 11 regions using NIH Epigenomics Roadmap Consortium data from a 108-day-old human fetal gut (15) (Fig. 1A). Among these, eight regions were marked by only one feature, one region by two and two regions by all three features (Table 1). For in vitro functional tests of these elements, we cloned them individually into a pGL4.23 luciferase vector with a minimal TATA-box of the β-globin gene and transfected them into the human neural crest-derived neuroblastoma cell line SK-N-SH. We have previously demonstrated that SK-N-SH expresses EDNRB and all other components of the RET GRN and is a model system for studies of transcriptional regulation in the ENS (7). For regions with multiple features, we used the overlap between their peaks to define the region of interest (Table 1). As shown in Figure 1B, the genomic elements E1 (3.5-fold; P = 2 × 10−3), E4 (1.8-fold; P = 1 × 10−3), E7 (8.2-fold; P = 1.9 × 10−11), E9 (4.8-fold; P = 3.1 × 10−8) and E10 (3.5-fold; P = 7.1 × 10−5) had significantly higher relative luciferase activity when compared to the basal promoter-only vector and represent five candidate fetal gut EDNRB enhancers. Element E10 overlaps the EDNRB transcriptional start site and appears to contain the potential promoter. Of these, the strongest enhancer activity was driven by E7, which overlaps a DHS peak and H3K27ac and H3K4me1 epigenetic marks.
Identifying putative enhancers at the EDNRB locus. (A) A 350 kb genomic segment comprising the core TAD in eight human cell lines containing 11 putative CREs defined using ENCODE DHS and H3K4me1 or H3K27ac marks from a 108-day-old human fetal intestine. The locations of all ENCODE TF ChIP-seq sites and the EDNRB gene are indicated. (B) In vitro luciferase assays in SK-N-SH cells show significant enhancer activity in 5 of the 11 putative CREs (marked in green) compared to a promoter-only control. Error bars represent SEs of the mean (*P < 0.01, **P < 0.001) for three independent biological replicates for each luciferase assay.
Eleven potential regulatory elements (hg19 coordinates) within a 341 kb EDNRB TAD marked by H3Kme1, H3K27ac and DHS marks in human fetal large intestine
| Element (size) | hg19 coordinates | Epigenomic marks |
|---|---|---|
| E1 (1587 bp) | chr13: 78353421-78 355 008 | DHS, H3K4me1, H3K27ac |
| E2 (654 bp) | chr13: 78 394 428-78 395 082 | H3K4me1 |
| E3 (1853 bp) | chr13: 78 406 337-78 408 190 | DHS |
| E4 (777 bp) | chr13: 78 410 031-78 410 808 | H3K27ac |
| E5 (1217 bp) | chr13: 78 427 391-78 428 608 | DHS |
| E6 (1253 bp) | chr13: 78 434 192-78 435 445 | DHS |
| E7 (1054 bp) | chr13: 78439541-78 440 595 | DHS, H3K4me1, H3K27ac |
| E8 (1111 bp) | chr13: 78 465 330-78 466 441 | H3K4me1 |
| E9(1322 bp) | chr13: 78 481 415-78 482 733 | DHS, H3K4me1 |
| E10 (1313 bp) | chr13: 78 493 407-78 494 720 | DHS |
| E11 (1125 bp) | chr13: 78 695 152-78 696 277 | DHS |
| Element (size) | hg19 coordinates | Epigenomic marks |
|---|---|---|
| E1 (1587 bp) | chr13: 78353421-78 355 008 | DHS, H3K4me1, H3K27ac |
| E2 (654 bp) | chr13: 78 394 428-78 395 082 | H3K4me1 |
| E3 (1853 bp) | chr13: 78 406 337-78 408 190 | DHS |
| E4 (777 bp) | chr13: 78 410 031-78 410 808 | H3K27ac |
| E5 (1217 bp) | chr13: 78 427 391-78 428 608 | DHS |
| E6 (1253 bp) | chr13: 78 434 192-78 435 445 | DHS |
| E7 (1054 bp) | chr13: 78439541-78 440 595 | DHS, H3K4me1, H3K27ac |
| E8 (1111 bp) | chr13: 78 465 330-78 466 441 | H3K4me1 |
| E9(1322 bp) | chr13: 78 481 415-78 482 733 | DHS, H3K4me1 |
| E10 (1313 bp) | chr13: 78 493 407-78 494 720 | DHS |
| E11 (1125 bp) | chr13: 78 695 152-78 696 277 | DHS |
These cis elements were assayed for enhancer activity in SK-N-SH cells; elements demonstrating significant activity over the basal promoter-only construct are marked in bold.
Eleven potential regulatory elements (hg19 coordinates) within a 341 kb EDNRB TAD marked by H3Kme1, H3K27ac and DHS marks in human fetal large intestine
| Element (size) | hg19 coordinates | Epigenomic marks |
|---|---|---|
| E1 (1587 bp) | chr13: 78353421-78 355 008 | DHS, H3K4me1, H3K27ac |
| E2 (654 bp) | chr13: 78 394 428-78 395 082 | H3K4me1 |
| E3 (1853 bp) | chr13: 78 406 337-78 408 190 | DHS |
| E4 (777 bp) | chr13: 78 410 031-78 410 808 | H3K27ac |
| E5 (1217 bp) | chr13: 78 427 391-78 428 608 | DHS |
| E6 (1253 bp) | chr13: 78 434 192-78 435 445 | DHS |
| E7 (1054 bp) | chr13: 78439541-78 440 595 | DHS, H3K4me1, H3K27ac |
| E8 (1111 bp) | chr13: 78 465 330-78 466 441 | H3K4me1 |
| E9(1322 bp) | chr13: 78 481 415-78 482 733 | DHS, H3K4me1 |
| E10 (1313 bp) | chr13: 78 493 407-78 494 720 | DHS |
| E11 (1125 bp) | chr13: 78 695 152-78 696 277 | DHS |
| Element (size) | hg19 coordinates | Epigenomic marks |
|---|---|---|
| E1 (1587 bp) | chr13: 78353421-78 355 008 | DHS, H3K4me1, H3K27ac |
| E2 (654 bp) | chr13: 78 394 428-78 395 082 | H3K4me1 |
| E3 (1853 bp) | chr13: 78 406 337-78 408 190 | DHS |
| E4 (777 bp) | chr13: 78 410 031-78 410 808 | H3K27ac |
| E5 (1217 bp) | chr13: 78 427 391-78 428 608 | DHS |
| E6 (1253 bp) | chr13: 78 434 192-78 435 445 | DHS |
| E7 (1054 bp) | chr13: 78439541-78 440 595 | DHS, H3K4me1, H3K27ac |
| E8 (1111 bp) | chr13: 78 465 330-78 466 441 | H3K4me1 |
| E9(1322 bp) | chr13: 78 481 415-78 482 733 | DHS, H3K4me1 |
| E10 (1313 bp) | chr13: 78 493 407-78 494 720 | DHS |
| E11 (1125 bp) | chr13: 78 695 152-78 696 277 | DHS |
These cis elements were assayed for enhancer activity in SK-N-SH cells; elements demonstrating significant activity over the basal promoter-only construct are marked in bold.
TFs regulating EDNRB
The specificity of enhancer function, as measured by its luciferase activity, of each of the 5 putative enhancers required identification of their cognate TFs. We first asked if SOX10, the major RET TF (7,16,17) regulates EDNRB as well. In the mouse, Sox10 directly regulates Ednrb in the developing ENS (18–21) by binding to a 540 bp element upstream of the gene (22, 23). First, we mapped this conserved (>70% base pair similarity between human and mouse) sequence to the human genome and discovered it mapped to E10, the potential EDNRB promoter (Fig. 2A). Second, we performed chromatin immunoprecipitation followed by quantitative polymerase chain reaction (ChIP-qPCR) in SK-N-SH for SOX10 to detect strong binding (19-fold enrichment; P = 6 × 10−3) to the SOX10 motif within this conserved sequence. Third, we demonstrated binding specificity by siRNA-mediated SOX10 knockdown and observing a 3-fold decrease (P = 2 × 10−3) in binding to the same sequence (Fig. 2A).
Identification of the cognate TFs bound to EDNRB enhancers. (A) Genomic map of the EDNRB locus with location for element E10 and ChIP-qPCR assays using a SOX10 antibody show enrichment of binding as compared to the background in SK-N-SH cells. The specificity of this binding is shown by siRNA knockdown of SOX10 with a concomitant reduction in the ChIP-qPCR signal. (B) Analogous assays for two putative GATA2 enhancers (E1 and E9) and ChIP-qPCR assays using a GATA2 antibody show enrichment of binding to E9 but not E1. (C) ChIP-qPCR assays on enhancer region E7 demonstrate specific binding of NKX2.5 but not NF-κB. Error bars represent SEs of the mean (*P < 0.01, **P < 0.001) for three independent biological replicates for each ChIP.
The common SOX10 regulator of RET and EDNRB prompted us to search for other shared TFs. We utilized published ChIP-seq data from the ENCODE project (24) to find that element E1 was bound by GATA2, RAD21 and CTCF, E4 by STAT3, E9 by GATA2 and E10 by CTCF and EZH2; we could not identify any known TF binding to element E7 in any of the ENCODE cell lines (24). Of these, GATA2 was of greatest interest because we have previously shown that this TF regulates RET and that a variant RET enhancer that fails to bind GATA2 is associated with HSCR (7).
For confirmation, we performed ChIP-qPCR for GATA2 in SK-N-SH cells for elements E1 and E9. To test for specificity, we showed that GATA2 strongly binds enhancer E9 with 26-fold enrichment over control IgG (P = 0.001) with a 2-fold decrease in binding with GATA2 knockdown using siRNA (P = 0.02) (Fig. 2B). We could not detect any binding of GATA2 at the E1 element as compared to the background in our assays, in contrast to ENCODE data (Fig. 2B). One potential reason for this discrepancy is that the ENCODE ChIP assays in SK-N-SH were performed using a different GATA2 antibody (http://antibodyregistry.org/search.php?q=AB_2616054).
Element E7 had the strongest reporter activity, was marked by all three epigenetic features and yet had no recognizable cognate TF in published data. We searched for binding sites within the E7 sequence using FIMO (25) and 890 validated TF motifs available in public motif databases, specifically the TRANSFAC motif matrix (26–29). We used the setting of “minimize false positives” to detect that only NKX2.5 (TCAAGTG; P = 1.3 × 10−6) and NF-κB (GGAAATTCCC; P = 0.008) had a matrix similarity and core similarity score of 1. The matrix similarity is a score that describes the quality of a match between a matrix and an arbitrary part of the input sequence. Analogously, the core similarity denotes the quality of a match between the core sequence of a matrix (i.e. the five most conserved positions within a matrix) and a part of the input sequence. No other TF had a matrix similarity score over 0.5, and the next best score was 0.45 for the TF AP1 within this sequence.
Of these, NKX2.5 was of primary interest because, although it is primarily a TF controlling early mesoderm differentiation and cardiac tissue morphogenesis in vertebrates, it also has a role in patterning in gut morphogenesis (30,31). In the mouse, Nkx2.5 along with Gata3 and Sox9 are expressed in undifferentiated cells in the pyloric mesenchyme, and loss of Nkx2.5 or Gata3 alters sphincter morphology, resulting in severe hypoplasia of a particular dorsal fascicle of longitudinal smooth muscles in the gut (32). Although NF-κB is also ubiquitously expressed in the developing gut and plays a significant role in the regulation of intestinal epithelial homeostasis and inflammation (33–35), it has no reported role in gut morphogenesis. We performed ChIP-qPCR for both NKX2.5 and NF-κB in SK-N-SH cells, independently, and detected significant binding (12-fold enrichment; P = 0.001) for NKX2.5 but no binding of NF-κB to E7 (0.28-fold enrichment; P = 0.8) (Fig. 2C). Additionally, we also demonstrated that NKX2.5 binding to E7 was reduced 3.6-fold (P = 8 × 10−3) under siRNA-mediated knockdown of NKX2.5 in SK-N-SH cells (Fig. 2C).
To prove that the identified TFs do indeed control EDNRB expression, we performed siRNA-mediated knockdown of each TF in SK-N-SH cells and measured both the TF transcript level and EDNRB gene expression, using Taqman-based qPCR assays. These data showed that EDNRB expression was indeed reduced by 38% (P = 4 × 10−4), 24% (P = 2.3 × 10−3) and 32% (P = 3.1 × 10−3) consequent to the knockdown of SOX10, GATA2 and NKX2.5 genes, respectively; as a control, knockdown of EDNRB by its specific siRNA reduced its expression 78% (P = 6.4 × 10−6) (Fig. 3A). Thus, each of these TFs control EDNRB gene expression (Fig. 3A). In contrast, knockdown of RARB, the remaining known RET TF (7), had no effect on EDNRB expression (Fig. 3A). For completeness, we checked whether RET expression was perturbed in NKX2.5 siRNA treated cells, as compared to control cells, and observed no change in gene expression (Supplementary Fig. 1). Thus NKX2.5 is an independent TF for EDNRB unlike SOX10 and GATA2.
TF-mediated in vitro control of gene expression. (A) siRNA-mediated knockdown of SOX10, GATA2 and NKX2.5 but not RARB in SK-N-SH cells downregulates EDNRB transcription. (B) siRNA-mediated knockdown of GATA2 affects activity of enhancer E9 but not E1, further confirming GATA2 binding to E9 only. Similar experiments, including using control siRNAs, show the specificity of binding of NKX2.5 and SOX10 to E7 and E10, respectively. (C) siRNA-mediated knockdown of GATA2, SOX10 and NKX2.5 demonstrates that loss of SOX10 also has an effect on activity of GATA2 regulated enhancer E9, highlighting feedback between the two TFs. The other enhancers are only affected by loss of their cognate TF. Pairwise comparisons are against a vector with basal promoter-only (black) for measuring enhancer activity and between untransfected and siRNA transfected cells for measuring TF specificity for both (B) and (C). Error bars represent SEs of the mean (*P < 0.01, **P < 0.001) for five independent biological replicates in all experiments.
We conducted additional tests to demonstrate that loss of each cognate TF leads to loss of enhancer activity. Specifically, we performed luciferase assays in SK-N-SH cells independently transfected with siRNAs against GATA2, SOX10, NKX2.5 or a negative control, to assess the activity of each putative enhancer. E1 did not show any loss of activity; E9 showed a 2-fold decrease in activity (P = 2.7 × 10−3), E7 a 1.9-fold drop in activity (P = 10−4) and E10 a 2.5-fold decrease (P = 3.8 × 10−4), when GATA2, NKX2.5 and SOX10 were knocked down, respectively (Fig. 3B). Thus, except for E1, the E7, E9 and E10 elements are enhancers that bind specific TFs and regulate EDNRB gene expression.
To test if there is a combinatorial activity and feedback between enhancers via the TFs, we also quantified the luciferase activity of each enhancer by knocking down all the TFs (GATA2, SOX10 and NKX2.5) independently with their cognate siRNAs. While enhancer E1 showed no effect, E7 and E10 had effects only when their cognate TFs, i.e. NKX2.5 and SOX10, respectively, were knocked down but showed no effects when other TFs were attenuated (Fig. 3C). Only GATA2 bound E9 enhancer showed a 1.5-fold decrease (P = 1.3 × 10−2) in activity when SOX10 was knocked down (Fig. 3C), highlighting a feedback, either direct or indirect between SOX10 and GATA2 in SK-N-SH cells.
In vivo effect of Sox10 in the developing mouse gut
Transcriptional studies in in vitro cell culture models are limited by the absence of multiple cells and tissue types, which plays a critical role in transcriptional regulation in vivo. To examine if our detected effect of SOX10 on EDNRB transcription in vitro is recapitulated in vivo during ENS development, we used a heterozygous Sox10 knockout mouse, which also develops aganglionosis (36,37). Multiple studies have previously demonstrated that Sox10 and Ednrb are epistatic in the developing mouse gut most likely through direct control of its transcription during ENCC migration and differentiation (18–22). We assayed for Ednrb, Sox10, Gata2, Nkx2.5 and Ret gene expression using Taqman-based qPCR in the developing mouse gut at E11.5 and E12.5, when the ENCCs cross the cecum and enter the colon. As expected, Sox10 is expressed at 50% levels, as compared to wild-type embryos at both stages, whereas Ednrb expression is also reduced by 30% at E11.5 (P = 3 × 10−5) and by 40% at E12.5 (P = 2.2 × 10−5). Additionally, we also detected 32% (P = 0.002) and 60% (P = 2.4 × 10−5) reduction in Gata2 expression at E11.5 and E12.5, respectively, recapitulating the Sox10-Gata2 feedback we had observed in SK-N-SH cells. There was also a 21% (P = 0.001) reduction in Nkx2.5 expression but only at E12.5. Expectedly, Ret gene expression is reduced by 40% at E11.5 (P = 3.6 × 10−5) and 50% (P = 2.5 × 10−6) at E12.5 in the Sox10 null heterozygote mouse gut, as we demonstrated previously (7) (Fig. 4). These results highlight that Sox10 not only controls the expression of Ednrb and Ret but also Gata2 and Nkx2.5 in a developmental stage-specific manner. Thus, a common regulatory network involving Ret and Ednrb, in the developing ENCC during enteric neurogenesis, is driven by Sox10 and Gata2 TFs.
RET–EDNRB GRN in the developing mouse ENS. Gene expression of Ednrb, Gata2, Nkx2.5 and Ret in the developing mouse gut at embryonic stage E11.5 and E12.5 in wild-type and Sox10 heterozygote embryos shows that Ednrb, Gata2 and Ret are transcriptionally affected at both developmental stages by loss of Sox10 expression. The effect on Nkx2.5 is only observed at E12.5. All pairwise comparisons are between wild-type and heterozygous embryonic guts within each developmental stage. Error bars represent SEs of the mean (*P < 0.01, **P < 0.001) for three independent embryos for each developmental stage and genotype.
Transcriptional feedback in the RET–EDNRB GRN
Every GRN is uniquely defined by the interactions between its constituent genes. Thus, we directly tested for feedback between EDNRB, RET and their TFs and tested whether this control was EDNRB dose-dependent. We varied EDNRB expression by varying the concentration of EDNRB siRNA and quantified the gene expression of SOX10, GATA2, NKX2.5, RARB, EDNRB and RET in SK-N-SH cells by Taqman qPCR (Fig. 5A). EDNRB siRNA concentrations up to 15 μM led to a maximum of 16% (P = 0.01) reduction, and concentrations up to 17 μm led to a maximum of 45% (P = 0.008) reduction in EDNRB expression but with no measurable expression changes in the other genes. However, at a concentration of 25 μM siRNA, EDNRB levels were reduced by 88% (P = 6.5 × 10−7), and expression of RET, SOX10 and GATA2 reduced by 32% (P = 0.006), 32% (P = 0.004) and 20% (P = 0.004), respectively. At 40 μM siRNA, the levels of EDNRB were reduced by ~90% (P = 2.7 × 10−8), while RET, SOX10 and GATA2 expression were unchanged as compared to the 25 μM siRNA dose; NKX2.5 and RARB levels remain unchanged throughout this concentration range having no effect from loss of EDNRB expression. It is interesting to note that, although there is reduced expression of SOX10, due to reduced EDNRB expression, this does not affect the levels of NKX2.5, even though we have observed Sox10 control of Nkx2.5 in the mouse gut at E12.5. We believe that there are two likely explanations: (1) there are additional regulatory controls on NKX2.5 in SK-N-SH cells as compared to the gut, and (2) the effect we see in the gut is developmental time-specific, which cannot be modeled in cultured cells. Thus, there is a positive dose-dependent feedback of EDNRB on the transcription of two of its TFs and RET, highlighting the bidirectionality of transcriptional control which lies at the heart of many regulatory networks.
Transcriptional feedback between RET, EDNRB and their TFs. (A) EDNRB gene expression in human SK-N-SH cells with increasing doses of EDNRB siRNA (12–40 μM). The transcriptional effect is seen in RET, SOX10 and GATA2 when EDNRB levels are significantly reduced to below 50% of wild-type levels. (B) Experiments as in (A) for RET loss of expression (12–25 μM) in SK-N-SH cells show a transcriptional feedback on EDNRB, SOX10 and GATA2 but only when RET expression falls below 50% of wild-type levels. NKX2.5 and RARB levels remain unchanged in both experiments. All pairwise comparisons are with transfections using control siRNAs. Error bars represent SEs of the mean (*P < 0.01, **P < 0.001) for five independent biological replicates.
Since EDNRB expression affects RET transcript levels as well, we investigated if the converse was true. We varied RET gene expression by varying the concentration of RET siRNA and quantified gene expression of the same genes as in the above (Fig. 5B). RET siRNA concentrations up to 15 μM led up to 50% reduction in RET gene expression (P = 2.35 × 10−5) with no significant loss of EDNRB expression. However, when RET siRNA concentration increased to 17 μM and 25 μM, RET expression further decreased to 25% (P = 4.1 × 10−8) and ~0% (P = 7.2 × 10−9), respectively, while EDNRB gene expression decreased to 80% at both levels (P = 0.003 at 17 μM and P = 0.0023 at 25 μM). Similarly, there was no change in SOX10 and GATA2 up to 15 μM of RET siRNA. SOX10 expression levels decreased more slowly to 80% (P = 0.05) and 20% (P = 2.0 × 10−9), respectively, while GATA2 decreased to 70% at 17 μM (P = 8 × 10−3) and 50% at 25 μM (P = 5.3 × 10−3). RARB and NKX2.5 levels were unchanged due to loss of RET expression.
Discussion
This study demonstrates that the two major genes controlling early ENS development and whose functional loss is associated with HSCR, RET and EDNRB are coregulated within one GRN. Specifically, although multiple enhancers regulate each gene (RET: 10 known, (7); EDNRB: 5 known, this study), the TFs GATA2 and SOX10 regulate both of them while RET and EDNRB demonstrate dose-dependent feedback on these same TFs, as well as on other genes within the GRN, at least in SK-N-SH cells. This bidirectional regulatory control common to both genes is the most parsimonious explanation for the strong epistasis between RET and EDNRB mutations observed in human and mouse aganglionosis (11,12) (Fig. 6).
The RET–EDNRB GRN. RET and EDNRB are the two major genes for ENS development and harbor multiple mutations leading to HSCR. These genes are coregulated within a larger GRN controlled by at least two common TFs with feedback and feed forward loops as significant features. The grey colored components of the GRN (GDNF, GFRA1 and CBL) were deciphered in our previous study (7).
Our study demonstrates that transcriptional feedback between EDNRB and its TFs is evident only when EDNRB gene expression levels fall below 50% (i.e. heterozygote levels). These results can therefore explain why, in the mouse, Ednrb hypomorphic and null mutations show aganglionosis penetrance only when its relative expression is <50% (38,39). Further, because Ednrb hypomorphic and null mutations show pigmentary anomalies arising from melanocytes but not aganglionosis in heterozygotes, the Ednrb gene dosage effect is probably absent in melanocytes. In other words, the effect of Ednrb loss of function mutations and its effect on the Ret–Ednrb GRN are cell-type specific in the mouse; hence, the same mutation will have different consequences in different cells. Our previous study of a EDNRB W276C hypomorphic mutation in HSCR in an Old Order Mennonite population, where its high frequency leads to both mutant heterozygotes and homozygotes, demonstrates dosage-dependent phenotypic defects analogous to the mouse (40). However, human EDNRB W276C heterozygotes, unlike the mouse, can also display aganglionosis, which likely results from a common SOX10-binding enhancer allele regulating both RET and EDNRB that segregates in the Old Order Mennonite population and predisposes to HSCR (40). Thus, parsimony dictates that the RET–EDNRB GRN is also cell-type specific in the human.
The effect of GATA2 and SOX10 on RET and EDNRB is minimally through the enhancers we have experimentally defined in SK-N-SH cells. However, the mechanism by which RET and EDNRB interact with GATA2 and SOX10 is unknown. This is unlikely to be through the canonical tyrosine kinase or G-protein coupled receptor signaling pathways since the genes in the RET signaling cascade are unaffected by siRNA knockdown in SK-N-SH cells or in the Ret null mutant homozygote mouse (7). We suspect that this is owing to specific posttranslational modifications of SOX10 and GATA2, in response to their lower expression, inhibiting their nuclear entry or making it inefficient. SOX10 is a nucleo-cytoplasmic shuttle protein and so provides some credence to this hypothesis (41).
We also hypothesize that sequence variants in at least the GATA2 and SOX10 binding enhancers, E9 and E10, would be HSCR-associated, but no such variants have been identified to date. Note that in the European ancestry population (42) a minimum of 146 polymorphisms (minor allele frequency [MAF], ≥1%) exist in the 11 putative enhancers reported here, 37 of which are in the 5 experimentally proven elements. The absence of associations refers to common (MAF, ≥10%) variants that have been tested (43,44). Rare variants could still be HSCR-associated, but identifying these will require more comprehensive studies.
Finally, for HSCR and gut development, gene expression of RET and EDNRB requires the activity of other TFs and enhancers involved, such as RARB and NKX2.5, and others yet to be identified. Thus, the RET–EDNRB GRN controlling ENS migration and differentiation, which currently comprises nine genes is undoubtedly much bigger. Note that EDNRB, GATA2 and SOX10 are expressed very early during neural crest cell specification and migration (45) and, subsequently, all the way through ENS development. Therefore, it is likely there would be both cell-type specific and cell-type invariant GRNs centered on these genes, as have been reported in melanocytes, migrating neural crest cells and HeLa cells (23). At the same time, loss of Ret expression in the developing ENS leads to hundreds of genes being affected (7). It is, therefore, not inconceivable that many genes jointly regulate RET and EDNRB and are commonly affected by their loss of function, and indeed, some of the observed effect on transcription could be mediated indirectly through other genes we have not uncovered yet. Thus, the GRN we have discovered is only the core of the network which involves multiple genes mutated in HSCR (Fig. 6) but is incomplete and clearly much bigger than what we have deciphered. Expanding the GRN would also shed light on other genes affected by and controlled by the specific TFs we have uncovered here and lead to a better understanding of how these interactions are wired. In other words, understanding the effect of sequence variation in ENS genes involved in HSCR and other neurocristopathies will require elucidating and understanding this dynamic GRN.
Materials and Methods
TAD definition and ChIP-seq peak calling
We utilized published HiC data (14) at the EDNRB locus from eight human cell lines and used Juicebox (46) to construct TADs by mapping the normalized data for each cell type at 5 kb resolution. The genomic coordinates (hg19) for TADs around EDNRB in each cell type are as follows: GM12878 (Chr13:78360000-78 700 000), HMEC (Chr13:78325000-79 865 000), HUVEC (chr13:78355000-78 700 000), IMR90 (chr13:78345000-78 700 000), K562 (Chr13:78705000-79 605 000), KBM7 (Chr13:78345000-78 705 000) and NHEK (Chr13:78350000-78 695 000). The common core TAD between these cells lines is a 341 kb region (Chr13:78359528-78 701 413).
Three epigenomic data sets from a 108-day-old human fetal large intestine, histone H3K27ac ChIP-seq (GSM1058765), histone H3K4me1 ChIP-seq (GSM1058775), and DNaseI-seq (GSM817188) were downloaded from the NIH Roadmap Epigenomics Project (http://www.roadmapepigenomics.org/data/tables/fetal). For each data set, MACS software v1.4 (47) with default settings was used to call “peaks” where sequence reads were significantly enriched. With the default peak-calling threshold (P < 10−5), 51 771, 61 689 and 66 930 genomic regions were identified in GSM1058765, GSM1058775, and GSM817188 respectively.
Regulatory element selection
Within the core TAD, we defined potential regulatory elements as segments overlapping peaks for any of the two histone marks or DHS; with multiple peaks, we used the overlap region. For region E9, which had a 3.9 kb overlap, we used the 1.3 kb region centered on the peak to keep all tested fragment sizes similar. For region E10, which overlapped a 5.7 kb DHS peak, we used the 1.3 kb fragment centered on the peak and overlapping the transcription start site of EDNRB.
Cell lines
The human neuroblastoma cell SK-N-SH, purchased from ATCC (no. HTB-11), was grown under standard conditions (DMEM +10% FBS and 1% penicillin streptomycin).
Reporter assays
Four hundred nanograms of firefly luciferase vector (pGL4.23, Promega Corporation) containing the DNA sequence of interest and 2 ng of Renilla luciferase vector (transfection control) were transiently transfected into cell line SK-N-SH (5–6 × 104 cells/well), using 6 ml of FuGENE HD transfection reagent (Roche Diagnostic, USA) in 100 ml of OPTIMEM medium (Invitrogen, USA). The cells were grown for 48 h, and luminescence was measured using a Dual Luciferase Reporter Assay System on a Tecan multidetection system luminometer, per the manufacturer's instructions. To check for the effect of loss of specific TFs on the reporter assay, we transfected cells with SOX10 siRNA (L-017192-00), GATA2 siRNA (L-009024-02), NKX2.5 siRNA (L-019795-00-0005) or a negative control non-targeting siRNA (D-001810-10) (Dharmacon) at a concentration of 25 μM for 24 h prior to transfecting the corresponding enhancer construct. All assays were performed in triplicate with independent readings (n = 9): the data presented are the means with their standard errors (SEs).
ChIP-qPCR assays
ChIP was performed thrice independently for each antibody using 1 × 106 SK-N-SH cells for each TF using the EZ-Magna ChIP kit (Millipore), as per the manufacturer's instruction, with the following modifications: the chromatin was sonicated for 30 s on and 30 s off for 10 cycles; sheared chromatin was preblocked with unconjugated beads for 4 h; and specific antibodies were separately conjugated to the beads for 4 h before immunoprecipitation was performed with the preblocked chromatin. The following antibodies were used: SOX10 (sc-17 342X, Santa Cruz), 10 mg; GATA2 (ab22849, Abcam), 8 mg; NKX2.5 (ab106923, Abcam), 5 mg; and NF-κB (17–10 060, Millipore Sigma), 5 mg. ChIP assays were also performed on cells 48 h after transfection with the following siRNAs at 25 μM to assess the specificity of TF binding: SOX10 (L-017192-00), GATA2 (L-009024-02) and NKX2.5 (L-019795-00-0005) (Dharmacon). qPCR assays were performed using SYBR green (Life Technologies) and specific primers against the E1 (E1_Fwd 5' CAAGAGGAAGTGCTGCCCTG 3'; E1_Rev 5' GCCTCTCTCAAGCACCTGAG 3'), E7 (E7_Fwd 5' CCCCAAGTATCAAGTGTAGCAGT 3′; E7_Rev 5' GATACATGATGTGTGCCCGA 3′), E9 (E9_Fwd 5' AGCTCTGCACTGAATCCAGG 3′; E9_Rev 5' CAGGCCTCCAGTGGTATGAA 3′) and E10 (E10_Fwd 5' TCCCTTCCCATCAATCACCG 3′; E10_Rev 5' ACTGATGCTGTCCAGGCATC 3′) regions. The data were normalized to input DNA, and enrichment was calculated by fold excess over ChIP performed with specific IgG as background signal. All assays were done in triplicate for each independent ChIP assay (n = 9).
Sox10 mouse
The generation of Sox10lacZ mice has been described previously (36). Heterozygous mice were intercrossed to generate all possible genotypes, and embryonic guts at E11.5 and E12.5 were dissected and genotyped; male embryonic guts were selected for further analysis. The tissues were washed in ice-cold phosphate buffered saline and snap frozen in liquid nitrogen for RNA extraction. Genotyping was performed on yolk sac DNA using primers specific to the Sox10 locus (Sox10Fwd: 5'-CAGGTGGGCGTTGGGCTCTT-3′ and Sox10Rev: 5'-CAGAGCTTGCCTAGTGTCTT-3′; 506 bp amplicon) for wild-type embryos and the lacZ transgene (LacZFwd: 5'-CTGGCAGGCGTTTCGTCAGTATCC-3′ and 5'-AGCGACATCCAGAGGCACTTCACCC-3′; 364 bp amplicon) for mutants. Genotyping for sex used primers mapping to the Kdm5c/d genes (Kdm5Fwd: 5'-CTGAAGCTTTTGGCTTTGAG-3′ and Kdm5Rev: 5'-CCGCTGCCAAATTCTTTGG-3′) resulting in two 331 bp bands specific for the X chromosome and an additional 301 bp Y chromosome band (48). All animal studies were reviewed and approved by the Institutional Animal Care and Use Committee of Johns Hopkins University (Protocol MO12M374) where these experiments were performed.
siRNA assays
EDNRB (L-003657-00-0005), RET (L-003170-00-0005), SOX10 (L-017192-00), GATA2 (L-009024-02), NKX2.5 (L-019795-00-0005) and RARB (L-003438-02) SMARTpool siRNAs (a combination of four individual siRNAs targeting each gene) along with ON-TARGET plus non-targeting siRNAs (D-001810-10, negative control) (Dharmacon, USA) were transfected at concentration ranges of 12 to 40 μM in SK-N-SH cells at a density of 104–105 cells using FuGene HD Transfection reagent (Promega Corporation, USA) per the manufacturer's instructions. Negative control siRNAs were always transfected at 25 μM concentration. To measure the efficacy of each gene specific siRNAs, we did a titration of each siRNA and measured gene expression by Taqman qPCR for its cognate gene. For SOX10, adding 15 μM, 17 μM and 25 μM siRNA led to transcript reduction by 28% (P = 3.2 × 10−3), 57% (P = 1.8 × 10−5) and 77% (P = 2.4 × 10−6), respectively. Corresponding assays for GATA2 led to reduction by 12% (P = 0.008), 49% (P = 2.2 × 10−4) and 79% (P = 1.4 × 10−6) and for NKX2.5, to reduction by 16% (P = 0.008), 40% (P = 1.2 × 10−4) and 72% (P = 1.3 × 10−5), respectively (Supplementary Fig. 2). For ChIP and luciferase assays in which the specific TF was knocked down, we used 25 μM of siRNA. The corresponding titration values for RET and EDNRB are detailed under results.
Total RNA was extracted from cells 48 h posttransfection and Taqman gene-specific assays conducted as described. Five independent transfections were used for each siRNA, and each Taqman assay was performed in triplicate (n = 15); P-values were calculated from pairwise two-tailed t tests, and the data were presented as means with their SEs.
Gene expression Taqman assays
Total RNA was extracted from SK-N-SH cells and individual mouse embryonic guts using TRIzol (Life Technologies, USA) and cleaned on RNeasy columns (QIAGEN, USA). Five hundred micrograms of total RNA was converted to cDNA using SuperScriptIII reverse transcriptase (Life Technologies, USA) using Oligo-dT primers. The diluted (1/5) total cDNA was subjected to Taqman gene expression (ThermoFisher Scientific) using transcript-specific probes and primers (Table S1). Human or mouse β-actin was used as an internal loading control, as appropriate, for normalization. Five independent biological samples for mouse fetal gut at each stage or five independent wells for SK-N-SH cells were used for RNA extraction and each assay performed in triplicate (n = 15). Relative fold change was calculated based on the 2DDCt (threshold cycle) method. For siRNA experiments, 2DDCt for negative control non-targeting control siRNA was set to unity; for measuring gene expression in mice guts, the 2DDCt value for E11.5 wild-type animals was set to unity. P-values were calculated from pairwise two-tailed t tests, and the data were presented as means with their SEs. Subsequently, P-values were calculated from pairwise two-tailed t tests, and the data were presented as the mean fold change with its SE.
Acknowledgements
The authors gratefully acknowledge the assistance of Dr. William J. Pavan (NHGRI, NIH) and Dr. Michael Wegner (Institut für Biochemie, Friedrich-Alexander-Universitat) for providing Sox10 mutant mice.
Conflict of Interest statement. None declared.
Funding
This work was supported by a National Institute of Health MERIT award [HD28088 to A.C.].
