PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Nat Neurosci. Author manuscript; available in PMC 2014 March 1.
Published in final edited form as:
PMCID: PMC3799941
NIHMSID: NIHMS505615

Saltatory remodeling of Hox chromatin in response to rostro-caudal patterning signals

Abstract

Hox genes controlling motor neuron subtype identity are expressed in rostro-caudal patterns that are spatially and temporally collinear with their chromosomal organization. Here we demonstrate that Hox chromatin is subdivided into discrete domains, controlled by rostro-caudal patterning signals that trigger rapid, domain-wide clearance of repressive H3K27me3 Polycomb modifications. Treatment of differentiating mouse neural progenitors with retinoic acid (RA) leads to activation and binding of RA receptors (RARs) to Hox1-5 chromatin domains, followed by a rapid domain-wide removal of H3K27me3 and acquisition of cervical spinal identity. Wnt and FGF signals induce expression of Cdx2 transcription factor that binds and clears H3K27me3 from Hox1-9 chromatin domains, leading to specification of brachial/thoracic spinal identity. We propose that rapid clearance of repressive modifications in response to transient patterning signals encodes global rostro-caudal neural identity and that maintenance of these chromatin domains ensures transmission of the positional identity to postmitotic motor neurons later in development.

Introduction

Development of a functional central nervous system (CNS) relies on mechanisms that assign precise positional identity onto dividing neural progenitors. Signaling factors released from localized sources within the embryo subdivide the nascent neural tube along two principal axes: dorso-ventral and rostro-caudal. Rostro-caudal patterning is established in early neural progenitors around the time of neural tube closure1 and in the developing hindbrain and spinal cord it is manifested by the differential expression of a phylogenetically conserved family of Hox transcription factors1-5. Although the activity of rostro-caudal patterning signals is transient, it has a lasting effect on the pattern of Hox gene expression that is maintained for the rest of embryonic development1. Maintenance of Hox gene expression is critical for the generation of segment-specific types of spinal nerve cells including columnar and pool specific motor neuron subtypes necessary for the proper wiring and functionality of spinal motor circuits2-4.

Mechanisms that control the establishment and maintenance of Hox gene expression have been under intense scrutiny. Mouse Hox genes are organized in four chromosomal clusters each harboring a subset of 13 paralogous Hox genes. Individual genes within these clusters are expressed in patterns that are spatially and temporally collinear with their physical chromosomal organization6,7. Mutations in the Polycomb group (PcG) complex of histone modifying enzymes lead to defects in the maintenance of Hox gene expression8, 9, indicating that epigenetic mechanisms and chromatin remodeling play an important role in the process of rostro-caudal patterning. Hox chromatin is found in a compacted and repressed state in pluripotent embryonic stem cells (ESCs)10 characterized by high density of histone H3 lysine 27 trimethyl modifications (K27me3) deposited by the PcG complex11, 12. The distribution of K27me3 modifications is altered during embryonic development resulting in a clearance of repressive modifications from chromatin domains harboring transcribed Hox genes13-15. Based on a correlation between chromatin modifications and Hox gene expression, it has been proposed that progressive removal of the repressive modifications from Hox chromatin controls the temporally progressive collinear pattern of Hox gene expression14.

Rostro-caudal patterning of the embryo is controlled by diffusible patterning signals. Specification of cervical, brachial and thoracic spinal cord identity depends on opposing gradients of retinoic acid (RA) and fibroblast growth factor (FGF) signals5, 16-18. RA secreted by paraxial mesoderm activates retinoid receptors (RARs) bound to genomic sites and initiates the recruitment of additional RARs within Hox chromatin domains harboring Hox1-Hox5 paralog genes19. FGF signaling activates expression of more caudal Hox6-9 paralog genes controlling the brachial and thoracic spinal identity5, 16-18. Caudal type homeobox protein Cdx2 is a candidate transcription factor linking FGF signaling to the regulation of caudal brachial and thoracic identity18. Wnt and FGF signals jointly induce the expression of Cdx220-22 and Cdx2 binding sites have been identified in a cis-regulatory element located between Hoxc8 and Hoxc9 genes23.

To study dynamic interactions between extrinsic signals and changes in chromatin architecture during rostro-caudal patterning, we developed an in vitro differentiation system that faithfully recapitulates normal development of cervical and brachio-thoracic spinal motor neurons. We provide evidence that patterning signals specifying cervical and brachio-thoracic identity activate RARs and Cdx2 transcription factors that are recruited to distinct Hox chromatin domains. We demonstrate that this pattern of repressive H3K27me3 modifications is rapidly altered following the exposure to rostro-caudal patterning signals. The repressive modifications are cleared from Hox chromatin domains occupied by RARs or Cdx2 transcription factors in a domain-wide saltatory process, instead of the proposed temporally progressive pattern of a shifting boundary between repressed and accessible chromatin. We propose that changes in H3K27me3 chromatin modifications are controlled by transient patterning signals and that stable maintenance of repressed and accessible Hox chromatin domains from progenitors to postmitotic motor neurons encodes the positional identity of differentiating cells and ensures proper specification of motor neuron subtype identity. Based on our findings we propose a comprehensive model of rostro-caudal patterning that integrates effects of extrinsic patterning signals, activation of developmentally regulated transcription factors, and changes in Hox chromatin modifications during neural development.

Results

Domain-wide changes in Hox chromatin modifications

To characterize mechanisms that underlie establishment and maintenance of the rostro-caudal pattern of Hox gene expression, we examined molecular mechanisms leading to specification of spinal motor neurons with distinct rostro-caudal positional identities. We took advantage of an in vitro mouse embryonic stem cell (ESC) to spinal motor neuron differentiation system that recapitulates normal motor neuron development in a highly homogeneous and purely neural context24. This system eliminates potentially confounding influences of mesodermal and endodermal lineages in the developing embryo. The high yield and efficiency of motor neuron differentiation makes this system amenable to analysis of temporal changes in chromatin modifications by chromatin immunoprecipitation (ChIP) analysis25. ESCs were converted to cervical spinal motor neurons by joint exposure of cells on Day 2 of differentiation to retinoic acid (RA) and ventralizing signal sonic hedgehog agonist (Hh)24. RA treatment triggers the recruitment of retinoic acid receptor (RAR) to genomic sites localized within the 3′ end of Hox clusters 19, leading to the specification of rostral cervical motor neurons expressing Hox4 and Hox5 paralog transcription factors by Day 6 of differentiation16, 24 (Fig. 1A).

Figure 1
Domain-wide clearance of repressive histone modifications and changes in Hox gene expression in response to RA

To understand how RA signaling is translated into a developmentally stable pattern of Hox gene expression, we examined changes in the distribution of repressive PcG mediated histone H3 modifications (H3K27me3 referred to as “K27me3”) that were previously implicated in the maintenance of Hox gene expression patterns. ChIP analysis revealed that the density and distribution of K27me3 modifications changed dramatically during the differentiation of ESCs into motor neurons. While the Hox clusters are covered by a high density of K27me3 modifications at ESC stage, the repressive mark was reduced to a background level in a 3′ Hox chromatin domain spanning Hox 1-5 genes in Day 7 motor neurons (Fig. 1B and Supplementary Fig. 1). The K27me3-depleted domain corresponds to the domain bound by RARs after the addition of RA and harbors all Hox genes expressed in cervical motor neurons (Fig. 1B)19. Together, these data indicate that large-scale changes in histone modification patterns accompany differentiation of ESCs into spinal motor neurons of defined rostro-caudal positional identity, recapitulating embryonic development.

Chromatin state analysis of tail buds at different developmental stages suggested a directional and temporally progressive removal of K27me3 modifications from Hox clusters 14. To determine the dynamics of K27me3 removal during neuronal development, we performed a time-series analysis of histone modifications during ESC differentiation into spinal motor neurons. The density of K27me3 modifications over Hox chromatin remained high until the addition of patterning signals on Day 2 of differentiation (Fig. 1C and Supplementary Fig. 2). After one day of RA/Hh treatment (Day 3 of differentiation), we observed a domain wide decrease in the density of K27me3 modifications spanning the Hox1-5 genes (Fig. 1C and Supplementary Fig. 2). By Day 4 (motor neuron progenitor stage) the pattern of K27me3 modifications over Hox genes was indistinguishable from the one observed in postmitotic motor neurons on Day 7 of differentiation (Fig. 1C and Supplementary Fig. 2). The rapid removal of repressive modifications is driven by RA treatment. Day 2 embryoid bodies that were treated with RA for 8 hours exhibited significant loss of repressive modifications in Hox1-5 domain compared to control untreated cells (Supplementary Fig. 2). Analysis of Ring1b and Suz12, components of Polycomb repressive complexes PRC1 (involved in chromatin compaction and repression of transcription 26 and PRC2 (methyltansferase complex required for the maintenance of H3K27me3), revealed that both complexes are rapidly displaced from activated Hox chromatin domain following RA treatment (Supplementary Fig. 3). Jointly, these results indicate that RA-mediated activation and recruitment of RARs to locations within Hox1-5 chromatin initiates saltatory (rapid, domain-wide) removal of Polycomb repressors and K27me3 histone modifications during motor neuron differentiation.

To directly compare the kinetics of K27me3 removal with transcription of Hox genes employing the same platform, we profiled histone H3 lysine 79 dimethyl modifications (“K79me2”). K79me2 modifications are correlated with the presence of elongating Pol2 complex and thus can be used as a proxy for transcriptional activity27. In contrast to the domain-wide changes in K27me3 density, a time-series analysis revealed temporally and spatially progressive accumulation of K79me2 modifications within the K27me3 depleted domain during spinal motor neuron differentiation (Fig. 1C). This is consistent with observed temporally progressive and collinear activation of Hox gene expression during motor neuron differentiation (Supplementary Fig. 2b) and reminiscent of the progressive collinear activation of Hox gene transcription following RA treatment of embryonal carcinoma cell lines 6. The dichotomy between the dynamics of K27me3 and K79me2 modifications is best illustrated by a comparison of Hoxa1 and Hoxa5 genes. While the rapid clearance of K27me3 from Hoxa1 coincides with K79me2 accumulation, the similarly rapid clearance of K27me3 modifications from the Hoxa5 gene is followed by a delayed, gradual accumulation of K79me2 modifications (Fig. 1D). Based on these observations, we conclude that Hox chromatin is partitioned into transcriptionally accessible (K27me3 low) and repressed (K27me3 high) domains in response to retinoid signaling. The rapid kinetics of K27me3 clearance cannot explain the temporally progressive pattern of Hox gene transcription indicating that the level of PcG repressive chromatin modifications is permissive rather than instructive for activation of Hox gene expression.

Suz12 activity is required to maintain Hox gene expression

To study the role of K27me3 modifications we asked whether ESCs harboring mutations in key PRC2 components (Eed, Ezh2, and Suz12) can be directed to differentiate into motor neurons in vitro. Genetic inactivation of PRC2 complex leads to the loss of K27me3 modifications and early embryonic lethality28. Likewise, ESC lines in which PRC2 function is disrupted (Eed−/−, Ezh2−/−)29, 30 failed to differentiate along neural lineage. However, we determined that previously reported β-gal gene-trap mutation of Suz12 31 is a hypomorph and ESCs homozygous for this allele maintain a low level of K27me3 (Fig. 2A and Supplementary Fig. 4). The hypomorphic line was able to differentiate into motor neurons, albeit with a lower efficiency (Fig. 2B). Analysis of Hox gene expression in Suz12βgal/βgal ESCs differentiated into cervical spinal motor neurons revealed a significant increase in the expression of brachial spinal cord Hox genes Hoxc6 and Hoxa7 (Fig. 2B and 2C) as compared to wild type cells. Consistent with this, a subset of the differentiated Suz12β-gal/β-gal cells acquired the identity of limb innervating motor neurons (Foxp1+, Raldh2+, Lhx3) found in brachial but not cervical spinal cord (Fig. 2B and 2D). Importantly, Suz12β-gal/β-gal mutant cells retained collinear expression of Hox genes during motor neuron differentiation (Fig. 2E). Together these findings provide compelling evidence that the PRC2 complex is critical for maintained repression of Hox genes located in K27me3-high chromatin territory and for proper specification of motor neuron subtype identity.

Figure 2
Patterns of Hox gene expression in Suz12βgal/βgal cells differentiated to motor neurons

Wnt3A and FGF2 signals caudalize motor neurons

Our observation that RA signaling leads to the clearance of repressive chromatin modifications only from a defined set of rostral Hox genes suggested that other patterning signals might control the removal of K27me3 modifications from more caudal Hox chromatin domains. To test this hypothesis, we examined mechanisms underlying the specification of more caudal spinal motor neurons. Wnt and FGF patterning signals are necessary for the specification of brachial and thoracic motor neurons expressing Hox6-9 transcription factors 16 and are sufficient to induce caudal spinal Hox gene expression in chick neural tube explants 17. We have systematically optimized caudalization of differentiating ESCs by Wnt3a and FGF2 treatment and identified conditions (150 ng/ml Wnt3a; 100 ng/ml FGF2) that gave rise to a mixture of motor neurons expressing the thoracic marker Hoxc9 (28.9 ± 3.1% of all motor neurons) and brachial markers Hoxc6 (20.8 ± 2.4% of all motor neurons) and Hoxc8 (50.4 ± 7.1 % of all motor neurons), as well as the limb innervating motor neuron marker FoxP1 (28.6 ± 3.9 % of all motor neurons) 32, 33 (Fig. 3A-C).

Figure 3
Wnt3A and FGF2 induce Cdx2 and caudalize differentiating motor neurons

To identify transcription factors that mediate induction of Hox6-9 genes by extrinsic Wnt and FGF signals, we focused on the evolutionarily conserved family of caudal homeobox (Cdx) transcription factors that are necessary for the development of caudal structures18, 34. Of the three mammalian Cdx homologues, Cdx2 exhibited the greatest degree of up-regulation (336 ± 22 fold) at 24 hours following Hh/RA/Wnt/FGF treatment (Fig. 3D). Activation of Cdx2 expression was dependent on a synergistic action of Wnt and FGF, as treatment of cells with Wnt3a alone16 or FGF2 alone failed to induce high levels of Cdx2 expression (Fig. 3D and Supplementary Fig. 5). Cdx2 protein expression was transient, dropping to basal levels by 48 hours of treatment (Fig. 3E). These results raise the possibility that transiently expressed Cdx2 might mediate the synergistic caudalizing activity of Wnt and FGF signals and regulate Hox gene expression.

Cdx2 induces caudal Hox expression in motor neurons

To determine whether Cdx2 activity is sufficient to activate caudal Hox gene expression independently of Wnt and FGF signaling, we generated a doxycycline (Dox) inducible V5-tagged Cdx2 ESC line (iCdx2) 25. Functionality of the epitope tagged Cdx2 is demonstrated by its ability to suppress Oct4 expression following its induction in ESCs 35 (Fig. 4A). Next, we asked if the expression of Cdx2 is sufficient to induce caudal Hox gene expression in the context of RA/Hh differentiated cervical motor neurons. Combined doxycycline and FGF2 treatment of RA/Hh treated embryoid bodies on Day 3 resulted in an efficient induction of Hoxc6, Hoxc8 and Hoxd9 expression (Fig. 4B-D). Expression of Cdx2 had no effect on Hox gene expression in the presence of FGF receptor inhibitor PD173074 (Supplementary Fig. 5B) and similarly, FGF2 treatment alone was not sufficient to induce caudal Hox gene expression (Fig. 4B and 4C). These results indicate that FGF signaling plays a dual role during neural rostro-caudal patterning - first it synergizes with Wnt to induce Cdx2 expression, and second it synergizes with Cdx2 to activate caudal Hox gene expression.

Figure 4
Cdx2 induces caudal Hox gene expression during motor neuron differentiation

Cdx2 binds and controls Hox chromatin modifications

If Cdx2 directly regulates Hox gene expression, it should bind to genomic locations within Hox clusters. We performed Cdx2 ChIP-seq experiments 48 hours after doxycycline treatment and identified 39,493 Cdx2 binding sites (p<0.001), including 118 binding sites within the Hox clusters (Fig. 5A). The primary DNA sequence motif overrepresented under Cdx2 ChIP-seq peaks corresponds to the previously described Cdx2 binding motif GYMATAAA 36 (Fig. 5B). Cdx2 selectively binds within the Hox1 to Hox9 paralog gene chromatin domain (Fig. 5C and Supplementary Fig. 6) and the 5′ extent of Cdx2 association with Hox chromatin is in register with the last 5′ Hox gene (Hox9 paralogs) induced in caudalized motor neurons (Fig 5C and Supplementary Fig. 6). Taken together, these results revealed that Cdx2 associates with regulatory regions proximal to transcriptionally active Hox genes, suggesting its involvement in the direct regulation of Hox gene expression during the specification of brachial and thoracic motor neurons.

Figure 5
Cdx2 directly controls Hox gene expression and chromatin modifications

To determine whether transient expression of Cdx2 caused changes in PcG mediated chromatin modifications, we measured the status of K27me3 by ChIP-seq analysis 36 hours after the Dox and FGF treatment (equivalent to 24 hours of Cdx2 expression). We observed a domain-wide removal of the repressive histone modifications up to a boundary between Hox9 and Hox10 paralogs (Fig. 5C and Supplementary Fig. 6). The K27me3 depleted domain spans the genomic territory occupied by Cdx2 transcription factor and harbors Hox genes expressed by Cdx2/FGF2 caudalized motor neurons (Fig. 5C and Supplementary Fig. 6). These results suggest that, analogous to RA signaling, transient expression of Cdx2 in neural progenitors results in a global change in the Hox chromatin landscape, making a subset of brachial and thoracic Hox genes available for transcription in postmitotic spinal motor neurons.

Discussion

We have determined a sequence of molecular events leading from the reception of transient patterning signals to the establishment of lasting rostro-caudal identity in neural tissue. Based on the combined analysis of extracellular signals, transcription factors, and chromatin states, we conclude that activity of RA and Wnt/FGF patterning signals is transmitted by RAR and Cdx2 transcription factors that bind to discrete Hox chromatin domains. Transcription factor binding is accompanied by a saltatory clearance of repressive K27me3 histone modifications, partitioning Hox clusters into activated and repressed domains. The K27me3 pattern is maintained through progenitor cell divisions until postmitotic motor neurons, defining a rostro-caudally appropriate subset of Hox genes available for expression and specification of motor neuron subtype identity. Consistent with the notion that patterning signals controlling postmitotic motor neuron Hox gene expression act early in development, we have recently demonstrated that Hox genes are not properly induced during transcriptional programming of motor neurons from ESCs that bypasses neural progenitor stages37.

Analyses of chromatin states in developing embryos and differentiating ESCs revealed that the temporally progressive activation of Hox gene expression is paralleled by a shift of the boundary between K27me3 low and high Hox chromatin domains and by a visible chromatin decondensation 14, 38. These results admit two interpretations: (1) the removal of K27me3 progresses one gene at a time from 3′ to 5′ ends of each Hox cluster during the period of rostro-caudal patterning, or (2) the different observed chromatin states are the product of discrete signaling events acting on distinct Hox chromatin domains. The results presented in this work favor the latter interpretation.

Time-course analysis of changes in Hox chromatin modifications revealed removal of repressive modifications that is rapid (within 24 hours following RA treatment) and that does not progress in a directional 3′ to 5′ manner at the timescale examined (see the 8 hour time point [Day2+ RA] in supplementary Fig S2). Importantly, kinetics of K27me3 removal does not correlate with progressive upregulation of Hox1-5 gene expression over a four day period. We propose that what appears like progressive directional removal of repressive modifications from Hox chromatin in vivo, is an effect of sequentially acting patterning signals, each initiating clearance of K27me3 modifications from a different subdomain of Hox chromatin. Early in development retinoids activate rostral Hox1-5 domain expressed in the hindbrain and cervical spinal cord. Later specified brachial and thoracic segments are exposed to Wnt and FGF signals emanating from the node region17, leading to Cdx2 induction and activation of more posterior Hox6-9 domains. We speculate that other combinations of signals including FGFs and Gdf115, 39 will be involved in activation of Hox10-Hox13 chromatin domain.

The rapid removal of K27me3 modifications is likely accompanied with relaxation of Hox chromatin and large scale chromatin remodeling 40. We observe that PRC1 protein Ring1b (Supplementary Fig. 3), necessary for the maintenance of condensed and repressed Hox chromatin26, is rapidly displaced from the activated domain. Thus our data are consistent with the reported formation of differential chromatin boundaries at different developmental times and rostrocaudal spinal cord positions3, 14, 41 and with observed large scale remodeling of Hox chromatin observed by FISH10, 26, 38.

Maintenance of K27me3 high domains is necessary for stable repression of Hox genes. While complete clearance of K27me3 modifications in Eed and Ezh2 null ESCs is incompatible with neural differentiation, erosion of K27me3 modifications in a Suz12 hypomorph cell line results in a de-repression of brachial Hox genes in cervical motor neurons. Ectopic expression of Hoxc6 and Hoxc7 leads to a change in a subtype identity of resulting motor neurons and to generation of ectopic limb innervating motor neurons expressing Foxp1. These findings are consistent with the observed de-repression of caudal Hox genes in Polycomb complex mutants8, 9, 11, including a recent study demonstrating that downregulation of PRC1 gene Bmi1 in the developing spinal cord in vivo leads to a de-repression of Hoxc9 gene in brachial spinal cord and a conversion of brachial motor neurons to thoracic ones41.

The final expression pattern of Hox genes in the nervous system is highly complex and cannot be explained only by binary division of Hox chromatin into permissive and repressed territories. Expression of individual Hox genes within the permissive chromatin domain is fine-tuned by secondary transcriptional interactions. For example, brachial Hox4-8 transcription factors located in the permissive chromatin domain engage in cross-repressive interactions leading to the diversification of motor neurons into distinct motor pool subtypes2. Similarly, Hoxc9 transcription factor is critical for direct repression of cervical and brachial Hox genes (Hox1-7) in permissive chromatin territory, ensuring generation of motor neurons with proper thoracic identity3. Interestingly, these secondary repressive events do not lead to reacquisition of K27me3 repressive Polycomb-catalyzed modifications3 (Supplementary Fig. 7).

Expression of Hox genes in developing limb mesenchyme is controlled by distal regulatory elements that are brought into physical proximity of Hox genes by chromatin looping42. However, the transcription factors binding to the distal regulatory elements have not been identified and the mechanisms through which these elements are recruited to the proximity of selected Hox genes remain unknown. We speculate that RAR and Cdx2 transcription factors bound within Hox genes might coordinate the long distance chromatin looping and provide the necessary local anchor points demarcating the precise chromatin domains to be activated during rostro-caudal patterning. Our observation that RAR and Cdx2 transcription factors bind in the vicinity of previously identified distal regulatory elements42 (Fig. 6B and data not shown), raises the possibility that the long distance chromatin looping is based on homotypic RAR-RAR and Cdx2-Cdx2 interactions, analogous to the previously described long distance interactions among chromatin sites occupied by estrogen receptors43.

Figure 6
The sequential activity of RAR and Cdx2 establishes two distinct chromatin states

Based on our findings and previous reports, we propose a four-step model of Hox gene regulation during rostro-caudal patterning of the developing neural tube (Supplementary Fig. 8). First, patterning signals activate relevant transcription factors that bind to cis regulatory elements within Hox clusters and to distal enhancers resulting in chromatin looping. Second, the transcription factors recruit transactivators that initiate rapid, domain-wide clearance of repressive H3K27me3 modifications from the bound chromatin domains leading to chromatin decondensation and transcriptional activation of Hox genes. Third, the established permissive/repressive chromatin domains are maintained until the postmitotic stage. Fourth, secondary transcriptional regulations of individual Hox genes within activated domains lead to fine-grain patterns of Hox gene expression observed in postmitotic motor neurons.

Together, this study provides insight into the precisely orchestrated interplay between signaling molecules, transcription factors, and chromatin modifications that jointly contribute to neural tube patterning. These findings open new approaches for directed differentiation of ESCs to the caudal spinal cord nerve cells by extrinsic signals or by forced expression of Cdx2 transcription factor. The study highlights the utility of ESCs as an optimal system to study mechanisms governing embryonic development, to test the inferred principles, and eventually to employ the developmental logic for the production of cell types that will advance our understanding of the nervous system in health and disease.

Methods

Cell culture and motor neuron differentiation

Hb9-GFP (HBG3) Suz12β-gal/β-gal and Eed mutant ES cells were differentiated as previously described 11, 24, 31. Briefly, ES cells were trypsinized and seeded at 5×105 cells/ml in ANDFK medium (Advanced DMEM/F12:Neurobasal (1:1) Medium, 10% Knockout-SR, Pen/Strep, 2 mM L-Glutamine, and 0.1 mM 2-mercaptoethanol) to initiate formation of embryoid bodies (Day 0). Medium was exchanged on Day 1, Day 2 and Day 5 of differentiation. Patterning of embryoid bodies was induced by supplementing media on Day 2 with 1 μM all-trans-Retinoic acid (RA, Sigma) and 0.5 μM agonist of hedgehog signaling (SAG, Calbiochem). For Wnt and FGF differentiation of ES cells, 150 ng/ml Wnt3A (R&D Systems) and 100 ng/ml bFGF (PeproTech), were added along with RA and Hh on day 2. To generate caudal brachial and thoracic ES motor neurons using inducible Cdx2 and FGF, imCdx2-V5 ES cells were induced to differentiate on day 2 with RA and Hh. On day 3, cultures were treated with Doxycycline (Dox) (2 μg/ml) for 6 to 8 hours. Cultures were collected, washed with PBS, and treated with 100 nM RA, 100 ng/ml FGF, and 100 nM Hh agonist and medium was changed on day 5. For ChIP experiments, the same conditions were used but scaled to seed 1×107 cells on Day 0. All differentiation was performed multiple times, analyzing at least three healthy cultures (n=1). Data distribution was assumed to be normal but this was not formally tested. No statistical methods were used to pre-determine sample sizes but our sample sizes are similar to those generally employed in the field.

Immunocytochemistry

Embryoid bodies were fixed with 4% paraformaldehyde in PBS, embedded in OCT (Tissue-Tek) and sectioned for staining: 24 hours at 4C for primary antibodies and 4 hours at RT for secondary antibodies. After staining, samples were mounted with Aqua Poly Mount (Polyscience). Images were acquired with a LSM 510 Carl Zeiss confocal microscope. Antibodies used in this study include: goat anti-Oct3/4 (ab27985, Abcam. 1:500), mouse monoclonal anti-Cdx2 (CDX2-88, BioGenex, 1:500), rabbit anti-Sox1 (kindly provided by Sara I. Wilson. 1:1,000), mouse anti-Isl1 and mouse anti-Hb9 (Developmental Studies Hybridoma Bank, 1:100), goat anti-Hoxc6 (Santa Cruz); rabbit anti-Hoxc4 (1:1,000), guinea pig anti-Hoxa5 (1:1,000), mouse anti-Hoxc8 mouse (1:100), guinea pig anti-Hb9 (1:5,000), rabbit and guinea pig anti-Foxp1 (1:1,000) are gifts from Tom Jessell. Alexa 488–(A11001, A11008,, FITC-(715-095-150, 706-095-148), Cy3-(715-165-150, 711-165-152, 706-165-148) and Cy5-(715-175-150, 711-175-152) conjugated secondary antibodies were used 1:2,000.

Quantifications of ES motor neuron cultures

To determine motor neuron differentiation efficiency and the fraction of Hoxc8 expressing ES motor neurons, EBs were dissociated on day 6 of differentiation. On day 7, we examined the number of Hb9+ cells and the total number of DAPI+ cells using immunocytochemistry or the expression of Hoxc8 in post-mitotic Hb9+ ES motor neurons.

Expression analysis

Total RNA was extracted from ES cells or embryoid bodies using Qiagen RNAeasy kit. For quantitative PCR (qPCR) analysis, cDNA was synthesized using SuperScript III (Invitrogen) and amplified using SYBR green brilliant PCR amplification kit (Stratagene) and Mx3000 thermocycler (Stratagene). For GeneChip expression analysis, RNA was amplified using Ovation amplification and labeling kit (NuGen) and hybridized to Affymetrix Mouse Genome 430 2.0 microarrays. Expression experiments were performed in biological triplicate for each analyzed time-point. Arrays were scanned using the GeneChip Scanner 3000. Data analysis was carried out using the affylmGUI BioConductor package 44. GCRMA normalization 45 was performed across all arrays, followed by linear model fitting using Limma 46. Differentially expressed genes were defined by ranking all probesets by the moderated t-statistic-derived p-value (adjusted for multiple testing using Benjamini & Hochberg’s method47) and setting thresholds of p<0.001 and a fold-change of at least 2. All arrays were submitted to the NIH GEO database under accession numbers GSE19372 (RA/Hh-dependent motor neuron differentiation timeseries) and GSE39422 (response to Cdx2 expression in motor neuron progenitors).

qPCR primers

Cdx1: acagccggtacatcactatcc, cttgtttactttgcgctccttg. Cdx2: tagtcgatacatcaccatcagg, tgattttcctctccttggctct. Cdx4: gcaatagatacatcaccatcagg, actttgcacggaacctccag. Hoxa5: tgtacgtggaagtgttcctgtc, gtcacagttttcgtcacagagc. Hoxc6: tagttctgagcagggcagga, cgagttaggtagcggttgaagt. Hoxc8: gtaaatcctccgccaacactaa, cgctttctggtcaaataaggat. Hoxc9: gcaagcacaaagaggagaagg, cgtctggtacttggtgtaggg. HPRT: agcaggtgttctagtcctgtgg, acgcagcaactgacatttctaa. B-actin: tgagagggaaatcgtgcgtgacat, accgctcgttgccaatagtgatga.

ChIP-chip protocols

ChIP protocols were adapted from http://jura.wi.mit.edu/young_public/hESregulation/ChIP.html. Briefly, approximately 3×10e7 cells taken from each developmental time point were cross-linked using formaldehyde and snap-frozen in liquid nitrogen. Cells were thawed on ice, resuspended in 5ml lysis buffer 1 (50 mM Hepes-KOH, pH 7.5, 140 mM NaCl, 1 mM EDTA, 10% glycerol, 0.5% NP-40, 0.25% Triton X-100) and mixed on a rotating platform at 4°C for 5 minutes. Samples were spun down for 3 minutes at 3000rpm, resuspended in 5ml lysis buffer 2 (10 mM Tris-HCl, pH 8.0, 200 mM NaCl, 1 mM EDTA, 0.5 mM EGTA), and mixed on a rotating platform for 5 minutes at room temperature. Samples were spun down once more, resuspended in lysis buffer 3 (10 mM Tris-HCl, pH 8.0, 100 mM NaCl, 1 mM EDTA, 0.5 mM EGTA, 0.1% Na-Deoxycholate, 0.5% N-lauroylsarcosine) and sonicated using a Misonix 3000 model sonicator to sheer cross-linked DNA to an average fragment size of approximately 500bp. Triton X-100 was added to the lysate after sonication to final concentrations of 1% and the lysate spun down to pellet cell debris. The resulting whole-cell extract supernatant was incubated on a rotating mixer overnight at 4°C with 100 μL of Dynal Protein G magnetic beads that had been preincubated for 24 hours with 10 μg of the appropriate antibody in a PBS/BSA solution. The following antibodies were used for ChIP experiments: antibodies against total H3 (Abcam, AB1791), H3K4me3 (Abcam, ab8580), H3K27me3 (Abcam, ab6002), H3K79me2 (Abcam, ab3594), and pan-RAR (Santa Cruz Biotechnology, sc- 773). After approximately 16 hours of bead-lysate incubation, beads were collected with a Dynal magnet. ChIP samples probing for transcription factor binding were washed with the following regimen, mixing on a rotating mixer at 4°C for 5 minutes per buffer: low-salt buffer (20 mM Tris at pH 8.1, 150 mM NaCl, 2 mM EDTA, 1% Triton X-100, 0.1% SDS), high-salt buffer (20 mM Tris at pH 8.1, 500 mM NaCl, 2 mM EDTA, 1% Triton X-100, 0.1% SDS), LiCl buffer (10 mM Tris at pH 8.1, 250 mM LiCl, 1 mM EDTA, 1% deoxycholate, 1% NP-40), and TE containing 50 mM NaCl. ChIP samples probing for histone and chromatin marks were washed 4 times with RIPA buffer (50 mM Hepes-KOH, pH 7.6, 500 mM LiCl, 1 mM EDTA, 1% NP-40, 0.7% Na-Deoxycholate) and then once with TE containing 50 mM NaCl, again mixing on a rotating mixer at 4°C for 5 minutes per buffer. After the final bead wash, samples were spun down to collect and discard excess wash solution, and bound antibody-protein-DNA fragment complexes were eluted from the beads by incubation in elution buffer at 65°C with occasional vortexing. Cross-links were reversed by overnight incubation at 65°C. Samples were digested with RNase A and Proteinase K to remove proteins and contaminating nucleic acids, and the DNA fragments precipitated with cold EtOH. Resulting purified DNA fragments were amplified via ligation-mediated (LM) PCR, labeled with a BioPrime CGH Genomic Labeling System (Invitrogen, 18095-011), and the labeled DNA co-hybridized to custom Agilent DNA microarrays.

Array design & hybridization

Labeled DNA samples were co-hybridized to custom Agilent DNA microarrays using the Oligo aCGH/ChIP-on-chip Hybridization Kit (Agilent, 5188-5220) at 65°C for approximately 16 hours. Arrays were washed according to previously published protocols (http://jura.wi.mit.edu/young_public/hESregulation/ChIP.html).

The Agilent 244K probe arrays were designed to tile at least 400Kbp surrounding each of the four Hox clusters with an average probe spacing of 110bp. The remaining probes were used to tile +/− 30Kbp around the transcription start site of selected genes exhibiting differential patterns of expression during motor neuron differentiation, with probe spacing between 100bp and 150bp. The array design and all ChIP-chip array data were submitted to the NIH GEO database under accession number GSE19447.

ChIP-chip data analysis

Arrays were scanned at dual PMT intensities (10% and 100%) at a 5um resolution using an Agilent microarray scanner. Feature extraction was performed with Agilent Feature Extraction software. Background subtracted values were normalized with (1) median normalization, (2) line fitting normalization, (3) quantile normalization; all normalization code was implemented in SQL as part of our in-house microarray database. Median normalization accounts for differences in the amount of dye hybridized between the two channels. This normalization multiplied IP intensities by median(control)/median(IP) such that the median intensity of the two channels is the same. Line fitting normalization is similar to Loess normalization but with a simpler model; it assumes that the bulk of the data should fall along the line y=x. Line fitting performed linear regression on the IP values as a function of the control and then rotated the datapoints such that the resulting line had slope one and intercept zero; this transformation is performed in log space. Finally, in order to allow comparison of probe intensities across time-points, arrays were quantile normalized 48. Heatmap style figures (Figures: 1C, S2) were generated in MATLAB by plotting the time-series of smoothed array probe intensities, where the smoothing calculated average intensities across a sliding window of 500bp and an offset of 250bp.

Histone modification quantification

A 2-state Hidden Markov Model was used to find the boundaries of H3K27me3-enriched domains at each time-point. Initial probabilities, transition probabilities, and Gaussian distributions for each state were estimated using the scaled Baum-Welch learner in the JaHMM Java library (http://www.run.montefiore.ulg.ac.be/~francois/software/jahmm).

To compare the dynamics of H3K27me3 and H3K79me2 at the promoter regions of Hox genes, we calculated the mean intensity of the probes in 1kb regions centered around the transcription start site of each gene at each time point. To determine the significance of differences in observed H3K27me3 and H3K79me2 levels at individual Hox genes, we infer the experimental variance for a given enrichment level by modeling mean-variance relationships across experimental replicates for all 1Kbp regions represented on the tiling arrays.

In order to determine whether a linear relationship exists between the dynamics of repressive histone modifications and Hox gene transcription after the addition of RA, we summarized the amount of H3K27me3 present near Hoxa1-5 by taking the average enrichment of the modification in a 4 kb window centered around the transcription start site (TSS) for each gene. The H3K27me3 enrichment and expression values (obtained by microarray analysis) were each normalized to the maximal value obtained for each gene over the course of the time series. A standard linear additive regression model was fitted to the H3K27me3/expression pairs from the entire time series using the least-squares method. The coefficient of determination (R2) for the fit of this model to these data is 0.7307. This may be interpreted as the proportion of variability in the data that is captured by the model. For the entire time series, a linear additive model captures a good portion of the variability in the data as noted by the observation that the onset of gene transcription occurs simultaneously with the initial domain wide removal of the repressive mark. We then limited our analysis to only the data obtained after the exposure to RA. R2 for the fit of the model to these data is equal to 0.0874. A much lower portion of the variability in these data is captured by a linear additive model. Thus, there is very little evidence of a linear relationship between H3K27me3 and expression after the initial onset of transcription that follows exposure to RA.

ChIP-seq protocols

Samples analyzed by ChIP-seq were prepared similarly to ChIP-chip samples except that purified DNA fragments were not amplified via LM-PCR but were instead processed according to a modified version of the Illumina/Solexa sequencing protocol (Illumina, http://www.illumina.com/pages.ilmn?ID=252). ChIP-sequencing reads for RAR were aligned to the mouse genome (version mm8) as described previously 19 using Bowtie49version 0.9.9.2 with options -k 2 --best. ChIP-sequencing reads for Cdx2 and H3K27me3 were aligned to the mouse genome (version mm8) using Bowtie version 0.12.5 with options -q --best --strata -m 1 -p 4 --chunkmbs 1024. Only uniquely mapping reads were analyzed further. Multiple hits aligning to an individual nucleotide position are discarded above the level expected at a 10−7 probability from a per-base Poisson model of the uniquely mappable portion of the mouse genome. RAR binding event analysis was performed as described previously 19, and Cdx2 binding event analysis performed using GPS 50. Motifs occurring at RAR and Cdx2 peaks were discovered using GimmeMotifs (with settings “-w 200 -a large -g mm8 -f 0.5 -l 500”), which combines results from the motif-finders MDmodule, MEME, GADEM, MotifSampler, trawler, Improbizer, MoAn, and BioProspector. Scaling ratios to allow comparison of H3K27me3 levels across experiments were determined using a regression analysis of read counts occurring in 10Kbp windows along the mouse genome. Raw sequencing data (FASTQ format) was submitted to the NIH GEO/SRA database under accession numbers GSE19409 (RAR), GSE39433 (Cdx2 and H3K27me3 in Cdx2 induced cells), and GSE47485 (H3K27me3 for PcG cell lines).

Supplementary Material

Acknowledgment

We would like to thank Vladimir Korinek (IMG, Prague) for providing us with purified Wnt3A protein, George Daley (Harvard Medical School) for providing Cdx2 construct, Michael Kyba (University of Minnesota) for sharing reagents for construction of the inducible Cdx2 ESC line, and Tom Jessell and Jeremy Dasen (NYU) for sharing antibodies and constructive discussion of our findings. EOM was the David and Sylvia Lieb Fellow of the Damon Runyon Cancer Research Foundation (DRG-1937-07) and the work was supported by Helmsley foundation grant, NIH grant P01 NS055923 (DKG, HW), R01 NS058502 (HW) and The Richard and Susan Smith Family Foundation, Chestnut Hill, MA (L.A.B.).

Footnotes

Author Contributions:

E.O.M., S.M., D.K.G. and H.W. conceived the experiments, analyzed the data and wrote the manuscript. E.O.M. generated and validated inducible cell lines, performed cell differentiations and expression analyses; S.M. and C.R. performed all computational and statistical analyses of genomic, expression and sequencing data; M.P. and T.P. performed and optimized caudalization experiments; S.R.T. and L.A.B. performed analysis of Prc2 null and hypomorph cell lines; S.McC. and R.A.Y. performed ChIP-chip experiments.

References

1. Ensini M, Tsuchida TN, Belting HG, Jessell TM. The control of rostrocaudal pattern in the developing spinal cord: specification of motor neuron subtype identity is initiated by signals from paraxial mesoderm. Development. 1998;125:969–982. [PubMed]
2. Dasen JS, Tice BC, Brenner-Morton S, Jessell TM. A Hox regulatory network establishes motor neuron pool identity and target-muscle connectivity. Cell. 2005;123:477–491. [PubMed]
3. Jung H, et al. Global control of motor neuron topography mediated by the repressive actions of a single hox gene. Neuron. 2010;67:781–796. [PMC free article] [PubMed]
4. Wu Y, Wang G, Scott SA, Capecchi MR. Hoxc10 and Hoxd10 regulate mouse columnar, divisional and motor pool identity of lumbar motoneurons. Development. 2008;135:171–182. [PubMed]
5. Liu JP, Laufer E, Jessell TM. Assigning the positional identity of spinal motor neurons: rostrocaudal patterning of Hox-c expression by FGFs, Gdf11, and retinoids. Neuron. 2001;32:997–1012. [PubMed]
6. Simeone A, et al. Sequential activation of HOX2 homeobox genes by retinoic acid in human embryonal carcinoma cells. Nature. 1990;346:763–766. [PubMed]
7. Papalopulu N, Lovell-Badge R, Krumlauf R. The expression of murine Hox-2 genes is dependent on the differentiation pathway and displays a collinear sensitivity to retinoic acid in F9 cells and Xenopus embryos. Nucleic Acids Res. 1991;19:5497–5506. [PMC free article] [PubMed]
8. Akasaka T, et al. Mice doubly deficient for the Polycomb Group genes Mel18 and Bmi1 reveal synergy and requirement for maintenance but not initiation of Hox gene expression. Development. 2001;128:1587–1597. [PubMed]
9. Beuchle D, Struhl G, Muller J. Polycomb group proteins and heritable silencing of Drosophila Hox genes. Development. 2001;128:993–1004. [PubMed]
10. Chambeyron S, Bickmore WA. Chromatin decondensation and nuclear reorganization of the HoxB locus upon induction of transcription. Genes Dev. 2004;18:1119–1130. [PubMed]
11. Boyer LA, et al. Polycomb complexes repress developmental regulators in murine embryonic stem cells. Nature. 2006;441:349–353. [PubMed]
12. Bernstein BE, et al. A bivalent chromatin structure marks key developmental genes in embryonic stem cells. Cell. 2006;125:315–326. [PubMed]
13. Rinn JL, et al. Functional demarcation of active and silent chromatin domains in human HOX loci by noncoding RNAs. Cell. 2007;129:1311–1323. [PMC free article] [PubMed]
14. Soshnikova N, Duboule D. Epigenetic temporal control of mouse Hox genes in vivo. Science. 2009;324:1320–1323. [PubMed]
15. Kashyap V, Gudas LJ. Epigenetic regulatory mechanisms distinguish retinoic acid-mediated transcriptional responses in stem cells and fibroblasts. J Biol Chem. 2010;285:14534–14548. [PubMed]
16. Peljto M, Dasen JS, Mazzoni EO, Jessell TM, Wichterle H. Functional diversity of ESC-derived motor neuron subtypes revealed through intraspinal transplantation. Cell Stem Cell. 2010;7:355–366. [PMC free article] [PubMed]
17. Nordstrom U, Maier E, Jessell TM, Edlund T. An early role for WNT signaling in specifying neural patterns of Cdx and Hox gene expression and motor neuron subtype identity. PLoS Biol. 2006;4:e252. [PMC free article] [PubMed]
18. Bel-Vialar S, Itasaki N, Krumlauf R. Initiating Hox gene expression: in the early chick neural tube differential sensitivity to FGF and RA signaling subdivides the HoxB genes in two distinct groups. Development. 2002;129:5103–5115. [PubMed]
19. Mahony S, et al. Ligand-dependent dynamics of retinoic acid receptor binding during early neurogenesis. Genome Biol. 2011;12:R2. [PMC free article] [PubMed]
20. Shimizu T, Bae YK, Hibi M. Cdx-Hox code controls competence for responding to Fgfs and retinoic acid in zebrafish neural tissue. Development. 2006;133:4709–4719. [PubMed]
21. Skromne I, Thorsen D, Hale M, Prince VE, Ho RK. Repression of the hindbrain developmental program by Cdx factors is required for the specification of the vertebrate spinal cord. Development. 2007;134:2147–2158. [PMC free article] [PubMed]
22. Chawengsaksophak K, de Graaff W, Rossant J, Deschamps J, Beck F. Cdx2 is essential for axial elongation in mouse development. Proc Natl Acad Sci U S A. 2004;101:7641–7645. [PubMed]
23. Taylor JK, Levy T, Suh ER, Traber PG. Activation of enhancer elements by the homeobox gene Cdx2 is cell line specific. Nucleic Acids Res. 1997;25:2293–2300. [PMC free article] [PubMed]
24. Wichterle H, Lieberam I, Porter JA, Jessell TM. Directed differentiation of embryonic stem cells into motor neurons. Cell. 2002;110:385–397. [PubMed]
25. Mazzoni EO, et al. Embryonic stem cell-based mapping of developmental transcriptional programs. Nat Methods. 2011;8:1056–1058. [PMC free article] [PubMed]
26. Eskeland R, et al. Ring1B compacts chromatin structure and represses gene expression independent of histone ubiquitination. Mol Cell. 2010;38:452–464. [PMC free article] [PubMed]
27. Krogan NJ, et al. The Paf1 complex is required for histone H3 methylation by COMPASS and Dot1p: linking transcriptional elongation to histone methylation. Mol Cell. 2003;11:721–729. [PubMed]
28. Surface LE, Thornton SR, Boyer LA. Polycomb group proteins set the stage for early lineage commitment. Cell Stem Cell. 2010;7:288–298. [PubMed]
29. Montgomery ND, et al. The murine polycomb group protein Eed is required for global histone H3 lysine-27 methylation. Curr Biol. 2005;15:942–947. [PubMed]
30. Shen X, et al. EZH1 mediates methylation on histone H3 lysine 27 and complements EZH2 in maintaining stem cell identity and executing pluripotency. Mol Cell. 2008;32:491–502. [PMC free article] [PubMed]
31. Pasini D, Bracken AP, Hansen JB, Capillo M, Helin K. The polycomb group protein Suz12 is required for embryonic stem cell differentiation. Mol Cell Biol. 2007;27:3769–3779. [PMC free article] [PubMed]
32. Rousso DL, Gaber ZB, Wellik D, Morrisey EE, Novitch BG. Coordinated actions of the forkhead protein Foxp1 and Hox proteins in the columnar organization of spinal motor neurons. Neuron. 2008;59:226–240. [PMC free article] [PubMed]
33. Dasen JS, De Camilli A, Wang B, Tucker PW, Jessell TM. Hox repertoires for motor neuron diversity and connectivity gated by a single accessory factor, FoxP1. Cell. 2008;134:304–316. [PubMed]
34. Moreno E, Morata G. Caudal is the Hox gene that specifies the most posterior Drosophile segment. Nature. 1999;400:873–877. [PubMed]
35. Niwa H, et al. Interaction between Oct3/4 and Cdx2 determines trophectoderm differentiation. Cell. 2005;123:917–929. [PubMed]
36. Berger MF, et al. Variation in homeodomain DNA binding revealed by high-resolution analysis of sequence preferences. Cell. 2008;133:1266–1276. [PMC free article] [PubMed]
37. Mazzoni EO, et al. Synergistic recruitment of transcription factors to cell type specific enhancers programs motor neuron identity. Nature Neuroscience. In the press. [PMC free article] [PubMed]
38. Chambeyron S, Da Silva NR, Lawson KA, Bickmore WA. Nuclear reorganisation of the Hoxb complex during mouse embryonic development. Development. 2005;132:2215–2223. [PubMed]
39. Liu JP. The function of growth/differentiation factor 11 (Gdf11) in rostrocaudal patterning of the developing spinal cord. Development. 2006;133:2865–2874. [PubMed]
40. Bickmore WA, Mahy NL, Chambeyron S. Do higher-order chromatin structure and nuclear reorganization play a role in regulating Hox gene expression during development? Cold Spring Harb Symp Quant Biol. 2004;69:251–257. [PubMed]
41. Golden MG, Dasen JS. Polycomb repressive complex 1 activities determine the columnar organization of motor neurons. Genes Dev. 2012;26:2236–2250. [PubMed]
42. Montavon T, et al. A regulatory archipelago controls Hox genes transcription in digits. Cell. 2011;147:1132–1145. [PubMed]
43. Fullwood MJ, et al. An oestrogen-receptor-alpha-bound human chromatin interactome. Nature. 2009;462:58–64. [PMC free article] [PubMed]
44. Wettenhall JM, Simpson KM, Satterley K, Smyth GK. affylmGUI: a graphical user interface for linear modeling of single channel microarray data. Bioinformatics. 2006;22:897–899. [PubMed]
45. Wu ZJ, Irizarry RA, Gentleman R, Martinez-Murillo F, Spencer F. A model based background adjustment for oligonucleotide expression arrays. J. Am. Stat. Assoc. 2004;99:909–917.
46. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3 Article3. [PubMed]
47. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Statist, Soc. B. 1995;57:289–300.
48. Bolstad BM, Irizarry RA, Astrand M, Speed TP. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 2003;19:185–193. [PubMed]
49. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25. [PMC free article] [PubMed]
50. Guo Y, et al. Discovering homotypic binding events at high spatial resolution. Bioinformatics. 2010;26:3028–3034. [PMC free article] [PubMed]