BAC-E-GFP recapitulates Runx3 expression in TrkC neurons from earliest developmental stages
Our previous study analyzed the pattern of reporter gene expression in the presence of reporter-expressing bacterial artificial chromosomes (BACs) that span the entire Runx3 locus [3]. Among these, BAC-E (220 kb in length) that extends 5’ ~166 kb upstream of Runx3 P1 promoter, and 3’ reached the beginning of exon 6 (Fig. 1a), conferred a high level of GFP expression in all Runx3-expressing tissues at E14.5, including DRG TrkC neurons [3]. BAC-E includes the three Runx3 enhancers (R1, R2, and R3) shown to confer expression in TrkC neurons from the earliest neurogenic stages [3]. BAC-E-derived GFP is detected from early developmental stages of TrkC neurons at E11.5 (Fig. 1b). Therefore, isolating Runx3 expressing TrkC neurons using transgenic mice expressing BAC-E-GFP is possible and informative.
Accessible chromatin and H3K27ac-marked regions in E11.5 Runx3-expressing DRG neurons
We analyzed the chromatin state by ATAC-seq to get a first glance of accessibility in regulatory regions of TrkC neurons at this early (E11.5) differentiation stage. This analysis revealed ~ 64,000 regions of open chromatin of which ~ 30% were in gene promoters [(up to ± 1 kb from the transcription start site (TSS)] and ~ 70% in intron/intergenic regions (Fig. 2a). H3K27ac marked genomic regions are well established as signatures of putative active transcriptional enhancers [13]; however, H3K27ac also marks many promoter regions [14]. Using the CUT&RUN method on E11.5 TrkC DRG neurons, revealed 32,655 H3K27ac-marked genomic regions, ~ 40% of which resided in promoter regions (Fig. 2a). Analyzing the extent of overlap between all ATAC and H3K27ac-marked regions revealed that most (~ 64%) of ATAC peaks in promoters and only some (26%) in intron/intergenic regions were also marked with H3H27ac (Fig. 2a). This result could be explained, at least in part, by the temporal changes in chromatin accessibility that tend to precede acquisition of the H3K27ac marker [15].
Genome-wide association of ATAC peaks with genes, using GREAT, revealed 4,557 genes with ATAC peaks only in their promoters (P-only), 5,431 genes with ATAC peaks only in their intron/intergenic regions (I-only) and 8,971 genes with ATAC peaks in both their promoters and intron/intergenic regions (P-I). We then determined whether these three gene groups differed in enriched GO terms. Interestingly, the P-I gene group was the only one prominently enriched for multiple neuronal terms under all GO term categories (Supplementary Table 1). The P-only genes showed minimal enrichment for some neuronal GO terms under “Cell component” (GO:CC) and KEGG. The I-only group of genes was mainly enriched for terms associated with immune/inflammatory processes, although a few neuronal terms were also enriched (Supplementary Table 1). Taken together, the identity of the E11.5 Runx3- expressing DRG neurons was most prominently observed in the P-I group of genes, which harbors ATAC peaks both in their promoters and intron/intergenic regions that contain genomic enhancers. Moreover, only the P-I group revealed neuronal identity among the corresponding three groups of H3K27ac-marked genes, (Supplementary Table 1).
Runx3-bound genomic regions in E11.5 TrkC neurons
The genome-wide landscape of Runx3 binding revealed 13,443 consensus Runx3 peaks and 9,601 of them (71%) harbored at least one RBS (Runx3_RBS). The genomic distribution of Runx3_RBS peaks showed that ~ 38% resided in promoter and 62% in intronic/intergenic regions (Fig. 2b). Association of these Runx3_RBS peaks with genes, using GREAT, revealed 2,776 genes with P-only peaks, 4,248 with I-only peaks and 1,356 genes with both (P-I). GO annotation enrichment analysis revealed prominent enrichment for neuronal terms in the P-I and I-only gene groups, whereas the P-only gene group showed minor enrichment for neuronal terms and prominent enrichment for general terms such as metabolic process and transcription (Supplementary Table 1). Notably, the proportion of I-only genes among all genes with Runx3_RBS-bound regions (51%) was much larger than that of I-only genes with ATAC or K3K27ac peaks (29% and 18%, respectively). These results reveal that the Runx3_RBS peaks in intron/intergenic regions are mainly responsible for conferring the neuronal identity of the isolated cells. Of note, 2,638 out of the 4,248 genes included in the I-only Runx3_RBS group did not overlap with the I-only ATAC or H3K27ac gene groups. It could be suggested that these 2,638 I-only genes with Runx3_RBS peaks represent early neuronal genes, which are not yet fully activated by the acquisition of open chromatin and the H3K27ac mark in genomic regions harboring enhancers.
Next, to assess the chromatin accessibility and H3K27ac marker in Runx3-bound regions, we determined the overlap of genomic Runx3_RBS-bound peaks with ATAC and H3K27Ac peaks. Almost all (98%) Runx3_RBS-bound promoter regions overlapped with ATAC and 90% of them were also marked by H3K27ac. Among Runx3_RBS peaks residing in intron/intergenic regions, only 30% overlapped with ATAC-only, and 28% were marked by ATAC and H3K27ac (Fig. 2b). These results indicated that the proportion of promoter regions marked by ATAC and H3K27Ac and overlap with Runx3_RBS (94.6%, Fig. 2b) was larger than in all promoter regions (63.7%, Fig. 2a). In contrast, the proportion of ATAC and H3K27ac marked intron/intergenic regions overlapping Runx3_RBS was like that in all intron/intergenic regions. Of note, both the ATAC and H3K27ac signals in Runx3-bound promoters and intron/intergenic regions were significantly stronger than in ATAC and H3K27ac peaks in regions without Runx3 binding (p < 2.2e-16, Kolmogorov Smirnov test, determined on the non-logarithmic values; Supplementary Fig. 1). These results suggest that Runx3 could increase the degree of chromatin “openness”, as known to be caused by pioneer TFs, reviewed in Zaret et al [16].
Additional TFs collaborate with Runx3 in regulating gene expression
Homer de novo TF motif enrichment analysis in Runx3 peaks revealed that RUNX was the most enriched motif. Moreover, the motifs for TFs essential for DRG neurons development, POU4F1/Brn3a and ISL1, were also enriched in all Runx3 peaks and Runx3_RBS peaks as well as in Runx3_RBS peaks residing in intron/intergenic regions (Fig. 3). In contrast, promoter-bound Runx3_RBS peaks were enriched for RUNX and ELK1, a member of the ETS family, that is known to interact with RUNX family TFs [17], but were not enriched for POU4F1, ISL1 or other TFs that are known to be important for neuronal differentiation (Fig. 3).
These results correspond with the finding that genes associated with Runx3 peaks in P-only regions were not enriched for neuronal GO terms, unlike those genes associated with Runx3-bound intron/intergenic regions. The results also suggested that Runx3 might cooperate with Brn3a and Isl1 in regulating gene expression in E11.5 TrkC neurons by binding to enhancer regions. This possibility was supported by the RNA-seq results that showed high level Pou4f1/Brn3a and Isl1 expression in E11.5 TrkC neurons [3]. Furthermore, our and others’ studies have shown that Brn3a regulates expression of Runx3 and several other genes including Cartpt, Dcc, Gal, Pirt, Piezo2 and Scn7a at early developmental stages of DRG neurons [3, 18, 19]. Similarly, Isl1 was shown to regulate the expression of Cckar, Cbln2, Gal, Nova1 and Scn10a in DRG neurons [20], genes that are differentially expressed in Runx3-deficient neurons [3].
We conducted Brn3a and Isl1 CUT&RUN experiments and compared the genome-wide distribution of their bound regions to those of Runx3-bound regions to further explore the possibility that Brn3a and Isl1 cooperate with Runx3. The results revealed 32,823 Brn3a-bound and 16,036 Isl1-bound regions, of which 23,671 (72.1%) and 11,451 (71.4%) peaks, respectively, harbored their respective consensus motifs V$BRNF and V$LHXF (Genomatix sequence tools). However, while 38% of Runx3_RBS peaks resided in promoters, only 5% of Brn3a_BRNF and 4% of Isl1_LHXF peaks did so.
As expected, Homer de novo motif enrichment analysis (Supplementary Fig. 2) confirmed that the most enriched motifs in Brn3a and Isl1-bound regions were POU4F1 and ISL1, respectively, corresponding to the Genomatix V$BRNF and V$LHXF TF family motifs, respectively. Of note, whereas the top 15 enriched motifs in Brn3a-bound regions included ISL1 and RUNX motifs, neither RUNX nor POU4F1 motifs were included in the top 15 enriched motifs in Isl1-bound regions (Supplementary Fig. 2).
Analysis of the peaks of all three TFs revealed that they can be divided into seven peak categories with various combinations of Runx3, Brn3a, and Isl1 TFs (Fig. 4). Interestingly, overlap analysis between the three TF_TFBS peaks categories indicated that 5,999 out of the 9,488 Runx3_RBS-containg peaks (63%) and 15,400 out of the 23,331 Brn3a_BRNF-containing peaks (66%) did not overlap with the two other TF-bound peaks (Fig. 4a, marked as “1 Runx3_RBS-only” and “2 Brn3a_BRNF-only”). Of all 11,282 Isl1_LHXF-containing peaks, 4,080 peaks (36%), marked as “3 Isl1_LHXF-only”, did not overlap with the two other TF-bound peaks (Fig. 4a). All three TF_TFBS categories overlapped in only 1,620 peaks, marked as “7 Runx3_RBS & Brn3a_BRNF & Isl1_LHXF” (Fig. 4a). Furthermore, only 2,862 (30%) and 1,985 (21%) of the 9,488 Runx3_RBS-containing peaks overlapped with Brn3a_BRNF and Isl1_LHXF-containing peaks, respectively (Fig. 4a, marked as 3 + 7 and 4 + 7, respectively). Of the 23,331 Brn3a_BRNF-containing peaks, 6,620 overlapped with Isl1_LHXF-containing peaks (marked as 6 + 7). The results indicated that while Brn3a and Isl1 peaks prominently overlapped, Brn3a overlapped with Runx3 more than with Isl1.
Overlap analysis of the three TF_TFBS peak types with ATAC marked chromatin regions revealed that 40–50% of regions with overlapping Brn3a_BRNF and Isl1_LHXF or Runx3_RBS and Brn3a_BRNF and Isl1_LHXF peaks (Fig. 4b, marked as 6 or 7, respectively) were marked only with ATAC, while the rate was only 10–25% in other peak categories (Fig. 4b). Interestingly, while ATAC and H3K27ac marked 26–50% of Runx3_RBS-containing peaks (Fig. <link rid="fig4">4</link>c, 1,4,5 and 7), they marked only 8–14% of Runx3-lacking peaks (Fig. 4c, 2,3 and 6). These results suggested that Brn3a and Isl1 might act as more efficient pioneer TFs than Runx3 by binding to inaccessible chromatin regions and initiate later chromatin opening.
Brn3a and Isl1 TFs are expressed at the transition from neuronal precursors to post-mitotic TrkC neurons [21]. Brn3a acts to suppress the progenitor state and drive differentiation to the TrkC fate by activating genes required for TrkC neuronal functions and suppressing genes of other neuronal fates [19]. Based on the Brn3a-deficient phenotype, it was suggested that Brn3a precedes and is the principal regulator of Runx3 expression [19]. While Isl1 is also highly expressed in early developing Runx3-expressing neurons [3, 20] and Isl1-deficiet embryos show a delay in TrkC expression [20], it hardly affects Runx3 expression and newborn mice do not display proprioceptive defects [20]. We have previously shown that Runx3 expression is regulated by three highly conserved distant enhancers upstream of its P1 promoter (Fig. 1). Two of these enhancers (R1 and R3) harbored Brn3a binding sites. Mutation in these sites markedly affected the activity of these enhancers [3]. These three enhancers offer good examples of the combined action of Brn3a, Isl1 and Runx3. The genomic landscapes of Brn3a and Isl1 binding in TrkC neurons revealed that all three Runx3 enhancers (R1, R2 and R3) harbor overlapping Brn3a_BRNF and Isl1_LHXF peaks (Fig. 5a), confirming their function as mediators of Runx3 expression. Interestingly, the Brn3a_BRNF and Isl1_LHXF peaks in R3 overlapped with Runx3_RBS (Fig. 5a). The R1 and R2 enhancers also bind Runx3. Still, they lacked RBS (Fig. 5a), raising the possibility that Runx3 is recruited to these enhancers by forming a physical complex with the other two TFs. These results suggested a complex interaction and cooperation among the three TFs in regulating Runx3 expression through these three enhancers. Furthermore, the genomic binding distribution of these three TFs suggested that Brn3a and Isl1 might also cooperate with Runx3 in regulating the expression of certain TrkC neuron genes by binding to the same or separate genomic region in these genes.
Identification of high confidence Runx3 directly regulated genes
Gene expression analysis was carried out by RNA-seq on E11.5 isolated Runx3-expressing DRG neurons from heterozygous Runx3-P2GFP/+ and Runx3-deficient homozygous Runx3-P2GFP/GFP mice. Using strict parameters (fold change ≥ 2, FDR < 0.05), we identified 464 differentially expressed genes (DEGs) [3]. Under less stringent parameters for DEGs identification (fold change ≥ 1.8, padj < 0.05), 625 DEGs were detected in the Runx3-deficient neurons, 150 DEGs were down-regulated and 475 up-regulated.
We defined Runx3 HCT as DEGs that harbored at least one Runx3_RBS peak within their gene body (promoter, introns, 5’-UTR, or 3’-UTR) and/or in the upstream or downstream intergenic region that spans the distance between it and its neighboring RefSeq genes. We determined the genes corresponding to the 9601 Runx3_RBS peaks using Homer, which associates each peak with only one gene, and whose TSS is closest to the peak, and GREAT (default parameters). The latter assigns to each gene a regulatory domain consisting of a basal domain that extends 5 kb upstream and 1 kb downstream from its TSS and an extension up to the basal regulatory domain of the nearest upstream and downstream genes within 1 Mb. The Homer and GREAT analyses revealed 6,709 and 8,730 Runx3_RBS peak-associated genes. We then intersected these peak-associated genes with the 625 DEGs, which revealed 224 and 273 putative Runx3-regulated genes, respectively, but only 212 genes were common to both lists. Screening the genes not included among these 212 putative Runx3 targets indicated that many of their Runx3_RBS peaks resided within introns of neighboring non DEGs. Therefore, we employed the UCSC genome browser, to which the Runx3_RBS peak coordinates were uploaded, to screen manually all 625 DEGs for the presence of Runx3_RBS peaks according to the parameters mentioned at the beginning of this section and identified the genes under direct Runx3 regulation. This screening revealed 244 Runx3 HCT, 86 down-regulated, and 158 up-regulated (Supplementary Table 2). Interestingly, the proportion of down-regulated Runx3 HCT among the down-regulated DEGs (86 out of 150 genes, 57.3%) was significantly higher than the proportion of the up-regulated targets among the up-regulated DEGs (158 out of 475 genes, 33.3%; Chi squared statistic with Yate’s correction 22.0461; p < 0.00001, significant at p < 0.01; Fisher exact test statistic at p < 0.00001, significant at p < 0.01).
Recently, a single-cell RNA-seq study on E11.5 DRG identified several clusters of various cell types, two of which were identified as Runx3-expressing early and late proprioceptive neurons [22]. Therefore, it was interesting to determine which of our 86 down-regulated and 158 up-regulated Runx3 HCT were in these proprioceptive neuron clusters. This assessment revealed that 76 out of the 86 down-regulated (88%) and just 23 of the 158 up-regulated (15%) Runx3 HCT were in the E11.5 proprioceptive neuron clusters. Most up-regulated HCT were not included in these proprioceptive or in the mechanoceptive neuron clusters (Supplementary Table 3). These results confirmed that most detected down-regulated Runx3 HCT were indeed proprioceptive neuronal genes, whereas the majority of up-regulated Runx3 HCT were not expressed in these neurons, due to their repression by Runx3. In line with these observations, the 86 down-regulated Runx3 HCT were enriched in GO annotation for neuronal terms. In contrast, the 158 up-regulated Runx3 HCT were enriched for immune system, cell death, and only marginally for neuronal-related terms (Supplementary Fig. 3).
Next, we determined the presence of Runx3, Brn3a, or Isl1 peaks within promoters and intron/intergenic regions of the 244 Runx3 HCT and the distribution of ATAC and H3K27ac peaks that overlap with these peaks. As described earlier for the genome-wide distribution of ATAC and H3K27ac peaks, most target genes promoter-residing peaks were marked by ATAC and/or H3K27ac regardless of whether they harbored Runx3_RBS peaks. Seventy percent of the Runx3-bound intron/intergenic peaks were marked by ATAC and/or H3K27ac. Brn3a-only peaks were the least marked by ATAC and /or H3K27ac. The highest percentage of epigenetic marks in intron/intergenic peaks coincided with overlapping peaks containing Runx3, Brn3a, and Isl1 (Supplementary Table 2). These results suggested that Brn3a and Isl1 might collaborate with Runx3 in promoting chromatin accessibility in the target genes, thereby regulating a significant proportion of the Runx3 HCT, either by shared binding to the same genomic region and/or by binding at separate genomic regions within these Runx3 HCT (Supplementary Table 2).
The role of Runx3 HCT in regulating proprioceptive neuron functions
Runx3 is known to act both as a transcriptional activator and as a repressor [23]. This dual mode of action fits very well Runx3 requirements during TrkC neuron diversification, namely, to allow the induction of gene-specific neuronal lineage and to shut off the expression of genes specific to previous stages and/or other cell types [21]. Several DEGs in Runx3-deficient TrkC neurons [3] were classified as proprioceptive or non-proprioceptive genes [22]. Some of these DEGs are direct targets of Runx3 that are down-regulated or up-regulated in its absence (Fig. 6).
One up-regulated target gene in Runx3-deficient cells was Ntrk2, which encodes for TrkB, and is normally expressed in TrkB but not in TrkC neurons. Runx3 is expressed in the pre-bifurcation proprioceptive/mechanoceptive neurons during early developmental stages. These neurons co-express the TrkB and TrkC surface markers of mechanoceptors and proprioceptors, respectively. As differentiation proceeds and these two lineages become distinct, Runx3 remains expressed only in TrkC proprioceptors [24].
In Runx3-deficient mice, the arly hybrid TrkB/TrkC neurons fail to segregate into single TrkC neurons, indicating that Runx3 is crucial for repressing Ntrk2 expression in cells destined to become TrkC neurons [25].
The resulting up-regulation of Ntrk2 in Runx3-deficient cells reflects its function as a transcriptional repressor during TrkC neurons development. It was suggested that Runx3 binding to a putative regulatory element that contains a cluster of Runx-binding sites might be responsible for Ntrk2 suppression in TrkC neurons [26]. We have now confirmed the prediction that Ntrk2 was indeed a direct target of Runx3. Runx3-deficient E11.5 TrkC neurons express elevated levels of Ntrk2 [3]. Furthermore, we found that Runx3 binds to five intronic regions in the Ntrk2 gene (Fig. 5b). One of these regions, located in intron 5, matches the putative regulatory element mentioned above. Moreover, Ntrk2 expression induced by Brn3a [19], correlated with nine intronic Brn3a_BRNF peaks. Four of the Runx3_RBS peaks overlapped Brn3a peaks. It is of interest to determine whether these sites are involved in suppressing Ntrk2 in TrkC neurons, while additional Brn3a-only peaks alone or in combination with Isl1 are responsible for Ntrk2 activation. The most likely scenario is that during the early developmental stages of DRG TrkC neurons, Brn3a induces Runx3, and both TFs cooperate, possibly also with Isl1, to down-regulate Ntrk2 expression. Ntrk2 is not the only TrkB neuron-specific gene repressed by Runx3. Several other up-regulated Runx3 target genes, depicted in Fig. 6 (Gfra2, Nr5a2, Pou6f2, and Sgk1), were also identified as mechanoceptive genes [22]. Therefore, suppression of the TrkB specific genes mentioned above is part of the program by which Runx3 consolidates the TrkC neuronal lineage.
As indicated above (Supplementary Fig. 3), the up-regulated Runx3 target genes included genes with immune-related functions. Some of these genes (Gata2, Il10ra, Milr1, Itga2b, and Zfp36) are depicted in Fig. 6. This finding indicated that suppression of these genes is an important Runx3 function, as it prevents inappropriate expression of immune-related genes during TrkC neuron development.
Runx3 is instrumental in inducing the expression of additional genes essential for the proper development and function of sensory neurons; some are selectively expressed in proprioceptive neurons, and some are not (Fig. 6). One of the down-regulated genes in the absence of Runx3 was Piezo2, which encodes a mechanosensitive ion channel [27]. We have shown that Piezo2 expression in Runx3-deficient E11.5 TrkC DRG neurons was down-regulated compared to Runx3-expressing neurons [3]. Moreover, studies in mice have shown that a mouse line that lacked Piezo2 in its proprioceptive neurons (Pvalb-Cre:Piezo2f/f) had severely uncoordinated body movements and abnormal limb positions [28]. Deletion of Runx3 in neurons (Wnt1-Cre/Runx3f/f; [29]) recapitulated the ataxia phenotype that occurred in Runx3-deficient mice [2], and was associated with severe skeletal scoliosis. Deletion of Piezo2 in proprioceptive neurons (Pvalb-Cre:Piezo2f/f) but not in osteogenic (Col1a-Cre:Piezo2f/f) or chondrogenic (Col2a-Cre:Piezo2f/f) cells resulted in skeletal defects, including spine malalignment and hip dysplasia [30]. Furthermore, human PIEZO2 deficiency results in proprioception defects and hip dysplasia [31].
Piezo2 is also down-regulated in the DRG of Brn3a-deficient and Isl1-deficient embryos [19]. It was interesting to note in the present study that Piezo2 was a Runx3 HCT with one Runx3_RBS peak in intron 2 that overlaps with Brn3a_BRNF and Isl1_LHXF peaks and resides in a chromatin-accessible region with ATAC and H3K27ac markers (Fig. 5c). Two other Runx3_RBS peaks, found in an intergenic region upstream of its TSS, overlapped with Brn3a_BRNF and Isl1_LHXF peaks, but not with ATAC or H3K27ac (Fig. 5c). Piezo2 also harbored 3 additional Brn3a_BRNF and Isl1_LHXF overlapping peaks and three Brn3a_BRNF-only peaks, none of which overlapped with Runx3_RBS (Fig. 5c). These findings can explain, at least in part, the ataxia and skeletal abnormalities seen in Runx3 and Piezo2-deficient mice.
Kirrel3, which encodes a synaptic adhesion molecule, was previously shown to be expressed as early as E11.5 in DRG TrkC neurons and was suggested to play a role in axonal pathfinding, cell recognition, and synapse formation of DRG neurons on appropriate target cells, including the targeting of proprioceptive neurons on muscle spindles through the interaction with Nephrin [32]. We have now found that Kirrel3 harbored five Runx3_RBS peaks, four of which in introns and one in an intergenic region − 280 kb from its TSS (Fig. 5d). Four of the Runx3_RBS peaks in Kirrel3 also overlapped with Brn3a_BRNF and Isl1_LHXF peaks. Kirrel3 also harbored 13 separate Brn3a_BRNF and three Brn3a_BRNF_Isl1-LHXF overlapping peaks (Fig. 5d).
These results suggested that Runx3 might cooperate with the highly expressed Brn3a and Isl1, in negatively regulating Ntrk2 and positively regulating the expression of Piezo2 and Kirrel3 in DRG proprioceptive neurons already at the early stages of TrkC neuron development. However, while Kirrel3-deficient mice displayed certain sensory abnormalities, autism-like and hyperactive behavior phenotypes, no ataxia was reported [33]. Therefore, although Kirrel3 plays a role in proprioceptive neuron targeting on muscle spindles, even its complete loss is insufficient to induce the severe proprioceptive defects seen in Runx3-deficient mice.
Distinct properties of Runx3, Brn3a, and Isl1 peaks in down-regulated and up-regulated Runx3 HCT
To determine whether an association between TF peak categories in the HCT and GO annotation exists, we performed K-means clustering (using R software) of the up-regulated and down-regulated Runx3 HCT according to their peak categories and subjected each cluster to GO annotation. Analysis of the up-regulated Runx3 HCT revealed that Cluster 1 consisted of 77 genes, including 55 harboring just Runx3-only peaks, and was enriched for immune system process terms. Cluster 2 included 27 genes that harbored mainly Runx3-only and Brn3a-only peaks, did not reveal any enriched GO terms and contained similar numbers of genes with neuronal and immune system functions (Fig. 7a). Clusters 3 and 4 comprised 54 genes, all of which contained Runx3, Brn3a, and Isl1 peaks and were enriched for neuronal terms (Fig. 7a). The down-regulated target genes formed two clusters enriched for Brn3a-only peaks. Cluster 1 also contained Runx3-only peaks, while Cluster 2 did not (Fig. 7b). Genes in both clusters also contained Isl1 peaks and GO annotation analysis revealed enrichment for neuronal terms (Fig. 7b). This clustering analysis revealed that the clusters in which most genes contained Runx3, Brn3a and isl1 peaks were enriched for neuronal terms. Taken together, these analyses indicated that Runx3 cooperates with Brn3a and Isl1 in activating or suppressing neuronal genes. In contrast, non-neuronal immune-related gene suppression was mainly done by Runx3, with minimal cooperation with Brn3a and Isl1.
Additional Runx3 HCT harboring Runx3-bound regions lacking RBS
Our definition of Runx3 HCT considered only DEGs associated with Runx3-bound regions that harbored a consensus RBS. However, we realized that many of these Runx3 HCT also harbored Runx3 peaks but no RBS (Runx3_noRBS). Moreover, whole genome analysis of Runx3 peaks revealed that 92% of the 3,842 Runx3_noRBS peaks reside in chromatin accessible regions (ATAC-marked) and the genes associated with these Runx3_noRBS peaks were enriched for GO neuronal terms (Supplementary Table 4). Therefore, we searched for DEGs associated only with Runx3_noRBS peaks that could be regarded as Runx3 HCT in addition to the 244 Runx3 HCT genes mentioned in the previous subsection. To that end, we intercrossed the lists of genes with Runx3_noRBS peaks, the 244 Runx3 HCT genes and the 625 DEGs. This analysis revealed 37 such genes, 11 down-regulated and 26 up-regulated in Runx3-deficient neurons (Supplementary Table 4). Interestingly, of their 41 Runx3_noRBS peaks, 39 were marked by ATAC, and 25 overlapped with Brn3a and/or Isl1 peaks (Supplementary Table 4). All 11 down-regulated genes (Adamts17, Ctdspl, Dio3, Faim2, Ndrg1, P4ha2, Pam, Pcbp3, Pcp4, Rgs4, and Slc6a15) and five of the up-regulated genes (Atxn10, Cst3, Mrpl34, Nwd2, and Rassf7) were expressed in the proprioceptive DRG neuron clusters 7 and/or 9 and some also in the mechanoceptive clusters 10 and/or 11 [22].
Ctss, a Runx3_noRBS containing up-regulated gene (Supplementary Fig. 4), that was not expressed in the proprioceptive or mechanoceptive clusters, encodes cathepsin S, a papain-like cysteine protease known for its immune function in antigen presentation. It is now understood that cathepsin S also has a role in activating itch and pain (nociception) in DRG neurons [34]. The nociceptive activity results from cathepsin S functioning as a signaling molecule that activates protease-activated receptors 2 and 4 members of the G-protein-coupled receptor family. The results suggested that Runx3 might regulate these 37 genes by recruitment to the DNA through interaction with Brn3a and Isl1, bringing the total number of Runx3 HCT genes to 281.