Global Reprogramming of Apoptosis-Related Genes during Brain Development

To enable long-term survival, mammalian adult neurons exhibit unique apoptosis competence. Questions remain as to whether and how neurons globally reprogram the expression of apoptotic genes during development. We systematically examined the in vivo expression of 1923 apoptosis-related genes and associated histone modifications at eight developmental ages of mouse brains. Most apoptotic genes displayed consistent temporal patterns across the forebrain, midbrain, and hindbrain, suggesting ubiquitous robust developmental reprogramming. Although both anti- and pro-apoptotic genes can be up- or downregulated, half the regulatory events in the classical apoptosis pathway are downregulation of pro-apoptotic genes. Reduced expression in initiator caspases, apoptosome, and pro-apoptotic Bcl-2 family members restrains effector caspase activation and attenuates neuronal apoptosis. The developmental downregulation of apoptotic genes is attributed to decreasing histone-3-lysine-4-trimethylation (H3K4me3) signals at promoters, where histone-3-lysine-27-trimethylation (H3K27me3) rarely changes. By contrast, repressive H3K27me3 marks are lost in the upregulated gene groups, for which developmental H3K4me3 changes are not predictive. Hence, developing brains remove epigenetic H3K4me3 and H3K27me3 marks on different apoptotic gene groups, contributing to their downregulation and upregulation, respectively. As such, neurons drastically alter global apoptotic gene expression during development to transform apoptosis controls. Research into neuronal cell death should consider maturation stages as a biological variable.


Introduction
Apoptosis is a ubiquitous regulated cell death pathway in multi-cellular organisms, essential for development, tissue homeostasis, and immune system functions [1]. Various physiological or pathological stimuli can induce apoptosis. Depending on the origin of a triggering stimulus, the signaling cascades leading to apoptosis are commonly attributed to the intrinsic (mitochondrial) pathway and the extrinsic (death receptor) pathway, though the two pathways can influence one another [2][3][4][5]. Both pathways converge on activation of executioner caspases that cause cellular structure breakdowns and DNA fragmentation [6][7][8].
The significance of apoptosis in the mammalian nervous system is well documented [9,10]. Dysregulated apoptosis is implicated in a wide spectrum of neurodevelopmental disorders and neurodegenerative diseases. During development, apoptosis plays an essential role in establishing appropriate neuronal cell numbers and forming robust neural circuits [11,12]. For example, peripheral neurons depend on target-derived trophic factors for survival and without them undergo apoptosis [13,14]. As such, the formation of neural circuits is a selection process through apoptotic control.
Despite the importance of apoptosis for neuronal development, neurons reduce in apoptotic competence as they mature [15,16]. In the peripheral nervous system, withdrawal Datasets of mRNA expression (mRNA-seq) and histone ChIP-seq were downloaded from the mouse Encyclopedia of DNA Elements (ENCODE) project (http://encodeproject.org) (accessed on 1 October 2021). These corresponded to three brain tissues (forebrain, midbrain, and hindbrain) of the Mus musculus C57BL/6 strain, and eight developing ages (embryonic days 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, and 16.5, as well as to postnatal day 0). The histone modification marks included H3K4me3 and H3K27me3 for the same brain regions and developing ages. Data from the two biological replicates of each experimental configuration were averaged.

RNA-Seq Data Analysis
Processed gene quantification data of two biological replicates from the forebrain, midbrain, and hindbrain tissues were downloaded from ENCODE. Expressed apoptosisrelated genes were defined using the following cutoffs: transcripts per million (TPM) ≥ 1 in at least one-third of data points and maximum TPM ≥ 5 among all data points. The Z-score normalization on log 2 (TPM + 1) was performed before the hierarchical clustering analysis using the UPGMA agglomeration method. Two major clusters were labeled as the upregulation cluster and the downregulation cluster based on their overall chronological trends. The same analyses were performed for data from the three regions separately.

ChIP-Seq Data Analysis
For the ChIP-seq data, "bed narrowPeak" files containing the peaks of signal enrichment and "bigwig" files containing signal p-values were obtained from ENCODE. For each sample, the file containing either replicated peaks or pseudoreplicated peaks was downloaded based on the default recommendation from ENCODE. Defined by ENCODE, replicated peaks are peak calls from pooled replicates and pseudoreplicated peaks are calls from two analysis partitions. Significant peaks were annotated with the ChIPpeakAnno package [40] using mm10 Ensembl genes as references. Peaks residing from −2 kb to 2 kb around a transcription start site (TSS) were defined as histone modification signals over the promoter region. Peak intensity data were obtained from the bed narrowPeak files. For a promoter overlapping with multiple peaks, the peak with maximum intensity was taken for further analysis. The ChIP-seq tracks were visualized using the UCSC genome browser.

Regression Fit between RNA-Seq (or ChIP-Seq) Data and Ages
Linear regression models were fit between expression levels (log 2 (TPM + 1)) of each gene in the upregulation or downregulation clusters and the considered age time points (P0 was coded as day 22). For the H3K4me3 or H3K27me3 ChIP-seq data, linear regression models were fit between the promoter peak intensity and ages. The age coefficient for each gene was denoted as β Ex , β K4 , and β K27 for expression, H3K4me3, and H3K27me3 signals respectively. The coefficient represents the change of the expression or the histone modification when the age advances one day.

KEGG Apoptosis Pathway
The classical apoptosis pathway (entry: hsa04210) was acquired from KEGG (Kyoto Encyclopedia of Genes and Genomes, https://www.genome.jp/kegg/, accessed on 1 October 2021). and then replotted for visualization purposes. Genes from the upregulation cluster in forebrain, midbrain, or hindbrain regions were labeled as "upregulation cluster" and genes from the downregulation cluster in any of the three brain regions were labeled as "downregulation cluster". Genes were color labeled if they positively or nega- To characterize expression dynamics of apoptosis-related genes in developing mouse brains, we extracted the expression levels (TPM) of 1923 apoptotic genes in the forebrain region across seven embryonic age points (E10.5, E11.5, E12.5, E13.5, E14.5, E15.5, and E16.5) and one postnatal (P0) age point. Many of these genes are well known for their roles in the classic apoptosis pathway, including caspases, Bcl-2 family proteins, death ligands and receptors, etc. Others are upstream modulators of the classical apoptotic genes, including signaling molecules and transcription factors, etc. The expression values from the ENCODE RNA-seq datasets were log-transformed (log 2 (TPM + 1)) before further Z-score normalization across developmental ages.
A total of 1306 apoptotic genes passed expression filtration (details in Materials and Methods) and were subject to hierarchical clustering, which unbiasedly defined five clusters ( Figure 1a). Clusters 1 and 4 were the two major gene groups whose gene expression, respectively, decreased and increased chronologically during brain development ( Figure 1b). Cluster 1 contained 603 downregulated apoptotic genes. The median fold change between the last (P0) and the first (E10.5) age points was 0.45 (P0 vs. E10.5). On the contrary, cluster 4 contained 434 apoptotic genes increasing their expression levels with a median 2.45-fold change.
The remaining genes belonged to three smaller clusters. Cluster 2 contained 77 genes, the expression of which were reduced during the early embryonic stage of development but appeared to be upregulated around birth. By contrast, cluster 5 contained 189 genes, modestly increasing expression in the early phase of forebrain development but decreasing expression from E16.5 to P0. Cluster 3 contained only three genes without a clear pattern of temporal regulation. In summary, a majority of the apoptotic genes (1037 out of 1306) falling within cluster 1 or 4 exhibited nearly unidirectional expression changes during mouse forebrain development. We called clusters 1 and 4 the "downregulation cluster" and the "upregulation cluster", respectively, for further analyses.

Temporal Expression Patterns of Apoptosis-Related Genes Are Consistent across Different Brain Regions during Development
To determine whether the dynamic regulation of these apoptosis-related genes was unique to forebrain development or ubiquitous to general neuronal development, we examined mouse midbrain and hindbrain transcriptomic RNA-seq data from the ENCODE project. Specifically, we determined the temporal expression patterns of the five forebrain clusters during midbrain and hindbrain development. Overall, the temporal patterns for each of these clusters followed a similar trend among the three brain regions ( Figure S1). In particular, forebrain cluster 1 displayed a downregulation trend and cluster 4 consistently exhibited an upregulation pattern in both the midbrain and the hindbrain regions.
We also conducted expression clustering ab initio based on the time-series RNAseq data from the midbrain and hindbrain regions. As shown in Figure 2a, a majority of apoptotic genes belonged to either a downregulation cluster (615 genes, or 47% of 1317 total expressed genes in the midbrains; 594 genes, or 44% of 1340 total expressed genes in the hindbrains) or an upregulation cluster (518 and 545 genes for midbrains and hindbrains, respectively; or 39% and 41%, respectively, of their total expressed genes).
More importantly, genes in the downregulation clusters largely overlapped between the three brain regions (Figure 2b). A total of 434 downregulated apoptotic genes (56.4% of the cluster union) were identified in all three brain regions. The genes in the upregulation clusters, likewise, shared 344 (50.7%) of apoptotic genes among the three regions. Genes in the upregulation and downregulation clusters of these three brain regions are listed in Table S2. This shows that the regulation of apoptotic genes during neurogenesis is largely

Temporal Expression Patterns of Apoptosis-Related Genes Are Consistent across Different Brain Regions during Development
To determine whether the dynamic regulation of these apoptosis-related genes was unique to forebrain development or ubiquitous to general neuronal development, we examined mouse midbrain and hindbrain transcriptomic RNA-seq data from the ENCODE project. Specifically, we determined the temporal expression patterns of the five forebrain clusters during midbrain and hindbrain development. Overall, the temporal patterns for each of these clusters followed a similar trend among the three brain regions ( Figure S1). In particular, forebrain cluster 1 displayed a downregulation trend and cluster 4 consistently exhibited an upregulation pattern in both the midbrain and the hindbrain regions.
We also conducted expression clustering ab initio based on the time-series RNA-seq data from the midbrain and hindbrain regions. As shown in Figure 2a, a majority of apoptotic genes belonged to either a downregulation cluster (615 genes, or 47% of 1317 total expressed genes in the midbrains; 594 genes, or 44% of 1340 total expressed genes in the

The Regulation of the Classical Apoptosis Pathway during Mouse Brain Development
We focused on genes in the classic apoptosis pathway, defined by the KEGG database, for their expression dynamics during brain development. These genes were most extensively reported for their roles in apoptosis, whereas other GO-annotated apoptotic genes were studied in a more limited context, for example, in a specific cell type or under a specific treatment condition. As shown in Figure 3, out of the 79 classic apoptosis genes, 39 genes displayed a decreasing expression trend during development in any of the brain regions, while 22 genes increased their expression chronologically. The large proportion of genes changing expression (61/79 or 77%) demonstrates, again, that differentiating neurons drastically alter their regulation of apoptosis during development. hindbrains) or an upregulation cluster (518 and 545 genes for midbrains and hindbrains, respectively; or 39% and 41%, respectively, of their total expressed genes). More importantly, genes in the downregulation clusters largely overlapped between the three brain regions (Figure 2b). A total of 434 downregulated apoptotic genes (56.4% of the cluster union) were identified in all three brain regions. The genes in the upregulation clusters, likewise, shared 344 (50.7%) of apoptotic genes among the three regions. Genes in the upregulation and downregulation clusters of these three brain regions are listed in Table S2. This shows that the regulation of apoptotic genes during neurogenesis is largely consistent across brain regions, indicating their robust transcriptional reprogramming is intrinsic to neuronal differentiation and ubiquitous to neuronal types.
In summary, both anti-and pro-apoptotic genes can be up-or downregulated, but nearly half of the developmentally regulatory events in the classical apoptosis pathway are the downregulation of pro-apoptotic genes, which makes up the largest regulatory gene group (28/(28 + 13 + 8 + 9) = 48%). Interestingly, in the downregulation cluster, classical pro-apoptotic genes were more downregulated than anti-apoptotic genes (foldchange medians: 0.38 vs. 0.67; p-value of the one-sided Wilcoxon signed-rank test: 0.04). In summary, both anti-and pro-apoptotic genes can be up-or downregulated, but nearly half of the developmentally regulatory events in the classical apoptosis pathway are the downregulation of pro-apoptotic genes, which makes up the largest regulatory gene group (28/(28 + 13 + 8 + 9) = 48%). Interestingly, in the downregulation cluster, classical pro-apoptotic genes were more downregulated than anti-apoptotic genes (foldchange medians: 0.38 vs. 0.67; p-value of the one-sided Wilcoxon signed-rank test: 0.04). Therefore, by number or magnitude, the pro-apoptotic genes are more downregulated than the anti-apoptotic genes.

Distinct Histone Modifications Underly the Developmental Regulation of Apoptotic Genes
To further quantify the temporal expression trends of apoptotic genes, we conducted regression fitting between the gene expression level of each gene and time points (see the Materials and Methods section for details) to derive its regression coefficient of gene expression (β Ex ) across development. The sign (+ or −) of a coefficient indicates the direction of the developmental change, whereas the absolute value of the coefficients indicate the "speed", or the "magnitude of the developmental change per time unit". As expected, the regression coefficients β Ex were negative for the downregulation clusters and positive for the upregulation clusters in all three brain regions (Figure 4a,b).  Regression coefficients (βEx) between mRNA expression and ages for the downregulation clusters (a) and the upregulation clusters (b) in the forebrain, midbrain, and hindbrain. Regression coefficients (βK4) between H3K4me3 signals and ages for genes in the downregulation clusters (c) and the upregulation clusters (d) in the forebrain, midbrain, and hindbrain. Regression coefficients (βK27) between H3K27me3 signals and ages for genes in the downregulation clusters (e) and the upregulation clusters (f) in the forebrain, midbrain, and hindbrain. All three brain regions show consistent patterns.

Developmental Remodeling of H3K4me3 and H3K27me3 Marks Underscores Expression Changes of Pidd1, Casp6, Mapk8, and Tradd
To illustrate the different epigenetic regulations as driving factors of developmental gene expression changes, we examined representative genes from the KEGG apoptotic pathway (Figure 3). For example, Pidd1 (P53-induced death domain protein 1) and Casp6 are two genes in the downregulation clusters of all three brain regions. PIDD1, caspase-2, and adaptor protein CRADD form a multiprotein complex, PIDDosome, that activates caspase-2 to initiate apoptosis [41,42]. CASP6 belongs to the cysteine-aspartic acid protease (caspase) gene family and is important for the execution phase of cell death [6][7][8]. As shown in Figure 5a,b, the H3K4me3 peaks at the promoter regions of Pidd1 and Casp6 diminished over time. By contrast, the H3K27me3 marks exhibited no pattern of regulation. As a result, Pidd1 and Casp6 mRNA expressions were downregulated 9-and 5-fold from E10.5 to P0 in the brain, respectively (Figure 5a,b). Notably, the H3K4me3 peak of Regression coefficients (β Ex ) between mRNA expression and ages for the downregulation clusters (a) and the upregulation clusters (b) in the forebrain, midbrain, and hindbrain. Regression coefficients (β K4 ) between H3K4me3 signals and ages for genes in the downregulation clusters (c) and the upregulation clusters (d) in the forebrain, midbrain, and hindbrain. Regression coefficients (β K27 ) between H3K27me3 signals and ages for genes in the downregulation clusters (e) and the upregulation clusters (f) in the forebrain, midbrain, and hindbrain. All three brain regions show consistent patterns.
To understand the mechanisms regulating the developmental expression dynamics of apoptotic genes, we investigated the temporal changes of histone mark H3K4me3, which is known to activate transcription and is highly enriched in gene promoter regions near transcriptional start sites [37]. We also examined repressive histone mark H3K27me3, which also resides around promoter areas but causes downregulation of nearby genes. We utilized the ENCODE ChIP-seq data corresponding to the same developmental ages and the same brain regions for precise inference. For each gene, we obtained the histone mark intensities at its promoter and quantified their changes along the brain development period through a regression fitting between histone mark intensities and time points (see Materials and Methods for details). The regression coefficients (β K4 , β K27 ) resulting from this analysis represented the change of histone modification intensities over a time unit.
For the downregulation clusters, we found the H3K4me3 coefficients β K4 mostly agreed with the developmental gene expression trends. Over three-quarters of genes displayed a negative β K4 , consistent with the notion that H3K4me3 is positively correlated with transcription (Figure 4c). For these genes, the loss of H3K4me3 at promoters during development implicates reduced transcription or expression levels over time. Interestingly, the repressive H3K27me3 mark showed insignificant contributions to the expression changes, as the median β K27 hovered around 0 and most values were relatively small (Figure 4d). In other words, the H3K27me3 patterns do not change much for most genes in the downregulation clusters.
The upregulation clusters exhibited notably different patterns of β K4 and β K27 than the downregulation clusters. If H3K4me3 alteration was the major driving force of expression changes, the upregulation cluster was expected to contain mostly positive β K4 . However, the H3K4me3 mark was not accordant with the expression changes, as median β K4 was about 0 for the forebrain and midbrain regions (Figure 4e). The hindbrain even atypically exhibited negative β K4 for most upregulated genes (median: −0.19). Therefore, H3K4me3 is not predictive of the developmental increase in expression of most apoptotic genes in the upregulation clusters.
By contrast, the loss of repressive H3K27me3 paralleled the increase in apoptotic gene expression. The β K27 of the upregulation cluster clearly displayed asymmetric regulatory directions, with over three-quarters of them in a negative range and the others in a negligibly positive range (Figure 4f). This meant, for most genes, peak intensities of H3K27me3 dropped over the developmental period. Taken together, the developmental upregulation of apoptotic genes can be attributed to the consistent removal of H3K27me3 marks from the promoters and, for only a selection of genes, to the addition of H3K4me3 marks.

Developmental Remodeling of H3K4me3 and H3K27me3 Marks Underscores Expression Changes of Pidd1, Casp6, Mapk8, and Tradd
To illustrate the different epigenetic regulations as driving factors of developmental gene expression changes, we examined representative genes from the KEGG apoptotic pathway (Figure 3). For example, Pidd1 (P53-induced death domain protein 1) and Casp6 are two genes in the downregulation clusters of all three brain regions. PIDD1, caspase-2, and adaptor protein CRADD form a multiprotein complex, PIDDosome, that activates caspase-2 to initiate apoptosis [41,42]. CASP6 belongs to the cysteine-aspartic acid protease (caspase) gene family and is important for the execution phase of cell death [6][7][8].
As shown in Figure 5a,b, the H3K4me3 peaks at the promoter regions of Pidd1 and Casp6 diminished over time. By contrast, the H3K27me3 marks exhibited no pattern of regulation. As a result, Pidd1 and Casp6 mRNA expressions were downregulated 9-and 5-fold from E10.5 to P0 in the brain, respectively (Figure 5a,b). Notably, the H3K4me3 peak of Pidd1 showed the largest drop between E10.5 and E11.5, consistent with the steepest decline in mRNA expression at the same time. Therefore, changes in H3K4me3 rather than H3K27me3 are the main contributor to developmental downregulation of Pidd1 and Casp6.
The upregulation cluster included Mapk8 and Tradd (Figure 5c,d). Mapk8 encodes mitogen-activated protein kinase 8, belonging to the MAP kinase family [43]. MAPK8 plays an important role in TNF alpha-mediated cell death and UV-induced apoptosis [44][45][46]. TRADD encodes TNFR1-associated death domain protein, which functions in the TNFR1 signaling pathway to regulate TNFα-induced extrinsic apoptosis and proteasomal stressinduced apoptosis [47][48][49][50][51]. For both genes, there were no clear trends of developmental changes in H3K4me3. However, their H3K27me3 signals markedly declined over time (Figure 5c,d). This led to an overall increase in mRNA expression of these two genes during development.
Pidd1 showed the largest drop between E10.5 and E11.5, consistent with the steepest decline in mRNA expression at the same time. Therefore, changes in H3K4me3 rather than H3K27me3 are the main contributor to developmental downregulation of Pidd1 and Casp6. The upregulation cluster included Mapk8 and Tradd (Figure 5c,d). Mapk8 encodes mitogen-activated protein kinase 8, belonging to the MAP kinase family [43]. MAPK8 plays an important role in TNF alpha-mediated cell death and UV-induced apoptosis [44][45][46]. TRADD encodes TNFR1-associated death domain protein, which functions in the TNFR1 signaling pathway to regulate TNFα-induced extrinsic apoptosis and proteasomal stressinduced apoptosis [47][48][49][50][51]. For both genes, there were no clear trends of developmental changes in H3K4me3. However, their H3K27me3 signals markedly declined over time (Figure 5c,d). This led to an overall increase in mRNA expression of these two genes during development.

Discussion
Neurons contain a unique regulatory program for apoptosis, as cell death needs to be attenuated to allow continuous neuronal survival. The regulatory program, as regards to apoptotic gene expression, is presumably determined during neuronal development. We found most apoptosis-related genes, including those in the classic apoptosis pathway,

Discussion
Neurons contain a unique regulatory program for apoptosis, as cell death needs to be attenuated to allow continuous neuronal survival. The regulatory program, as regards to apoptotic gene expression, is presumably determined during neuronal development. We found most apoptosis-related genes, including those in the classic apoptosis pathway, are subject to developmental regulation, consistent with the idea that neurons transform how they regulate apoptosis during differentiation. This developmental reprogramming is conserved, sharing most of the upregulated and downregulated genes between different brain regions (forebrain, midbrain, and hindbrain), and, therefore, is intrinsic to the process of neural differentiation and ubiquitous to various neuronal types. Interestingly, neither the up-nor downregulatory trend is restricted to one category of anti-or pro-apoptotic genes (determined by GO terms). This is probably not surprising given the vast number of apoptosis-related genes considered in the study.
However, the core genes in the classical apoptosis pathway (defined by KEGG) exhibit skewed regulation. Specifically, among the pro-apoptotic genes, more are downregulated than upregulated. Among the downregulated genes, pro-apoptotic genes reduce expression more than anti-apoptotic genes do. Notably, while anti-apoptotic Bcl-2 family proteins (Bcl-2 and Bcl-xl) are upregulated, most pro-apoptotic Bcl-2 family proteins (Bax, Bak, Noxa, Bim) and caspases (Casp2, Casp6, Casp8, Casp9) are downregulated. This would reduce the overall "fuel" for apoptosis, agreeing with the notion that the balance and relative ratio of pro-and anti-apoptotic genes determine cell death outcomes and that neurons attenuate apoptosis during development [15,16].
Casp3 and Casp7 are the exceptions, as they increase expression during embryonic development. They are best known as the executioner caspases of apoptosis for destruction of cellular proteins and structures, irrespective of upstream death signals. In neurons, caspase-3 is also important for dendritic trimming and synaptic pruning, which is necessary during the early postnatal period to establish robust and responsive neural circuits [52][53][54][55][56]. Furthermore, caspase-3 has non-apoptotic functions for synaptic plasticity in adult neurons [57,58]. How this localized function of caspase-3 is activated and restricted without inducing a full-blown cell death program is an interesting subject. The essential role of caspase-3 for neural circuit formation likely restrains its developmental downregulation, or even requires its upregulation so as to be distributed broadly in distant dendritic processes.
Despite upregulation of effector caspases Casp3 and Casp7, neurons uniformly turn down initiator caspases (Casp2, Casp8, and Casp9). This constitutes an effective strategy to attenuate apoptosis, because, in healthy cells, caspases are zymogens with weak or little protease activity and a hierarchy of caspase activation is needed to amplify the destructive signals leading to apoptosis [7,[59][60][61]. For example, cytochrome c, released from mitochondria, interacts with Apaf1 to form apoptosome, a complex that recruits caspase-9, causing its self-cleavage and activation, which in turn cleaves and activates caspase-3. Cytochrome c, Apaf1, and Casp9 all decrease mRNA expression during development (Figure 3), effectively reducing the chance of caspase-3 activation.
Another regulatory checkpoint restricting apoptosis is the Bcl-2 family of proteins. Many apoptotic signals converge on BAX and BAK1 in mitochondria by inducing their oligomerization or formation of mitochondrial transition pores, leading to mitochondrial membrane permeabilization and release of cytochrome c and DIABLO [62,63]. BCL-xL and BCL-2 can interact with BAX and BAK1 to inhibit mitochondrial membrane permeabilization [64,65]. Intracellular ratios between anti-and pro-apoptotic BCL-2 family proteins determine cells' responsiveness to death stimuli [66,67]. This ratio is magnified by the upregulation of Bcl-xl and Bcl-2 and simultaneous downregulation of Bax and Bak1, thereby significantly diminishing neurons' vulnerability towards apoptosis activation. Notably, Bak1 is subject to further post-transcriptional regulation. Increased splicing of Bak1 exon 5 during development generates an alternative isoform, degraded by the nonsense-mediated mRNA decay pathway without productive translation, such that the BAK1 protein is nearly absent in neurons [24]. The abatement of apoptosome components and the increasing ratio of anti-to pro-apoptotic Bcl-2 family members present two critical "choke" points to attenuate apoptosis in neurons.
The developmental regulation of apoptotic genes coincides with widespread changes in histone modification marks at their promoter regions. This developmental rearrangement of H3K4me3 and H3K27me3 for apoptotic genes may constitute a larger epigenome program. The H3K4me3 landscape continuously changes during brain development and has been implicated in directing neuronal lineage differentiation [38,68,69]. H3K27me3 regulation has also been shown to orchestrate several aspects of neurogenesis, including the balance of self-renewal and differentiation of neural progenitors, as well as the switch from neurogenesis to gliogenesis [70][71][72][73].
We found the upregulation and downregulation clusters of apoptotic genes to undergo distinct dynamic histone modification mark transitions. The predominant driving force for the developmental downregulation of apoptotic genes is reducing H3K4me3 signals on their promoters where H3K7me3 does not significantly change. By contrast, downregulation of H3K27me3 is responsible for developmental upregulation of apoptotic genes, while changes in H3K4me3 are not predictive of the increased gene expression. The differential re-modeling of H3K4me3 and H3K27me3 during development for the upregulation vs. downregulation clusters suggests divergent epigenetic regulatory mechanisms of apoptotic genes, which is an interesting topic for future investigation.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10.339 0/cells10112901/s1, Figure S1: Expression of apoptotic genes in midbrain and hindbrain regions. Table S1: GO terms related to apoptosis. Table S2: Apoptotic genes in the downregulation and upregulation cluster in the three brain regions. Table S3: Accession IDs of RNA-seq and CHIP-Seq data analyzed in this study.
Funding: This research was funded by NIH grant number R01NS104041 and R01MH116220 to S.Z. and R01GM137428 to L.C.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
All datasets used in this study can be accessed at the ENCODE portal (http://encodeproject.org, accessed on 1 October 2021), using the accession IDs provided in Table S3.