Cell-type and fetal-sex-specific targets of prenatal alcohol exposure in developing mouse cerebral cortex.

Prenatal alcohol exposure (PAE) results in cerebral cortical dysgenesis. Single-cell RNA sequencing was performed on murine fetal cerebral cortical cells from six timed pregnancies, to decipher persistent cell- and sex-specific effects of an episode of PAE during early neurogenesis. We found, in an analysis of 38 distinct neural subpopulations across 8 lineage subtypes, that PAE altered neural maturation and cell cycle and disrupted gene co-expression networks. Whereas most differentially regulated genes were inhibited, particularly in females, PAE also induced sex-independent neural expression of fetal hemoglobin, a presumptive epigenetic stress adaptation. PAE inhibited Bcl11a, Htt, Ctnnb1, and other upstream regulators of differentially expressed genes and inhibited several autism-linked genes, suggesting that neurodevelopmental disorders share underlying mechanisms. PAE females exhibited neural loss of X-inactivation, with correlated activation of autosomal genes and evidence for spliceosome dysfunction. Thus, episodic PAE persistently alters the developing neural transcriptome, contributing to sex- and cell-type-specific teratology.

Introduction

Congenital neurodevelopmental disabilities are increasingly common. A recent National Health Interview Survey (NHIS) found that 17.8% of children surveyed in 2015–2017 were reported to have at least one developmental disability, including attention-deficit/hyperactivity disorder (ADHD, 9.5%), autism spectrum disorders (ASD, 2.5%), and intellectual disabilities (ID, 1.2%). Moreover, those rates were higher than those reported in 2009–2011 and 2012–2014 surveys (Zablotsky et al., 2019). Perturbations in the maternal-fetal environment are likely to be a significant, but unresolved, contributor to the increasing prevalence of neurodevelopmental disorders (Modabbernia et al., 2017). Recent state-wide biomarker-based studies show that third trimester prenatal alcohol exposure (PAE) in the US occurs in ∼8% of births (Bakhireva et al., 2017; Umer et al., 2020), a much higher estimate than seen with self-report in the National Survey on Drug Use and Health (England et al., 2020). Therefore, alcohol is a common maternal-fetal perturbagen, and likely a leading cause of congenital anomalies and neurodevelopmental disabilities. PAE is known to result in growth and neurodevelopment deficits, collectively termed fetal alcohol spectrum disorders (FASDs). FASDs are common, with a recent case ascertainment study conservatively estimating a prevalence of 1.1%–9.8% of school-aged children in the US (May et al., 2018). Intriguingly, individuals with diagnoses along the FASD continuum are also documented to be co-morbid for ASD, ADHD (Bishop et al., 2007; Infante et al., 2015; Lange et al., 2017; Mukherjee et al., 2011; Riikonen et al., 2005; Stevens et al., 2013), and a number of other diagnostically distinct neuropsychiatric conditions (Popova et al., 2016). The prevalence of these comorbidities in FASD populations is significantly higher than in the general population (Popova et al., 2016), suggesting that many congenital neurodevelopmental disorders have a shared etiology, with origins in adverse maternal-fetal environments.

Heavy PAE is a well-established cause of brain growth deficits in human populations (Norman et al., 2009), and these outcomes have been replicated in mouse models of PAE (e.g. (Marquardt and Brigman, 2016; O'Leary-Moore et al., 2010)). We hypothesize that the critical period of telencephalic neurogenesis, which extends from the mid-first trimester through the second trimester in humans (Workman et al., 2013), is a particularly important window of vulnerability for PAE-induced brain growth deficits. A number of birth-dating studies (Takahashi et al., 1995, 1999) have shown that the pseudostratified epithelium of the murine fetal ventricular zone (VZ) generates most neurons of the cortical plate between gestational days 11 and 17, with founder neural stem cell (NSC) populations completing approximately 50% of neurogenic cycles by GD 14.5. Mouse models of PAE between gestational days 12.5 and 14.5, corresponding to the first half of the cortical neurogenic period, and the generation of layer V and VI neurons (Takahashi et al., 1995), resulted in ventricular dilation and decreased thickness of the developing cortical plate (Sudheendran et al., 2013; Tsai et al., 2014). We also observed in ex vivo models of the murine cortical neuroepithelium that although ethanol exposure did not kill neural stem cells, it did result in increased asymmetric cell proliferation, loss of cells expressing markers of stem cell identity, and increased indices of premature maturation (Camarillo and Miranda, 2008; Prock and Miranda, 2007; Santillano et al., 2005; Tingling et al., 2013). These data suggested that PAE during the period of neurogenesis results in diminished brain growth due to premature stem cell maturation and, consequently, loss of NSCs.

In the current study, we used a single-cell RNA sequencing (scRNAseq) approach to assess the persistent effects of a single episode of binge-like maternal ethanol exposure on the developmental trajectory of fetal cerebral cortical cells. ScRNAseq enables the characterization different cell types based on individual cell transcriptional profiles (Loo et al., 2019; Macosko et al., 2015) and has been used to successfully assess cell-type specific mediators of neurodevelopmental disorders, including ASD (Anderson et al., 2020; Co et al., 2019; Fazel Darbandi et al., 2020) and Rett syndrome, including allele-specific X chromosome activation status, in humans and in mouse models (Renthal et al., 2018). We also assessed whether the effects of PAE on the cortical development transcriptome were dependent on fetal sex, as a previous scRNAseq study has shown sex differences in adult cortical NSC differentiation trajectories (Mizrak et al., 2019) and a limited amount of research in rodent models has indicated that sex differences occur in PAE-associated behaviors that appear later in life (Raineki et al., 2016; Varlinskaya and Mooney, 2014). Therefore, these behavioral outcomes may be mediated by early molecular adaptations during neurodevelopment.

Our data, focused on the proliferative cells of the VZ, the zone of secondary proliferation, i.e., sub-ventricular zone (SVZ), and layer V/VI neurons, indicate that each of these subpopulations are molecularly and developmentally heterogeneous, with up to 22 different cell types identified within the VZ and SVZ and 8 different neuronal cell types assigned to layers V/VI. PAE significantly increased cell cycle, mainly in SVZ-type cells, and altered developmental trajectory in VZ, SVZ, and neuronal populations, and inhibited a network of genes associated with protein synthesis, while broadly inducing hemoglobin gene expression across neural sub-populations. Moreover, we identified significant sex differences in the response to PAE, which in females was linked to loss of X-inactivation and differential expression of neurodevelopmental-disorder-associated genes.

Results

Sample characteristics

Pregnant dams, at GD 12.5, were exposed to ethanol vapor for 45 min, during which they exhibited incomplete sedation and impaired mobility with loss of righting reflex. This exposure resulted in a blood ethanol content (BEC) between 210 and 260 mg/dL at 30 min following the end of the exposure period (Table 1). All dams regained their righting reflex and mobility by 90 min following the exposure period.

Table 1

Dam and sample IDs and maternal blood ethanol levels

Dam ID	Maternal blood ethanol concentration (BEC)a	Sample IDs
9	260 mg/dL	PAE_male_1PAE_female_1
10	0 mg/dL	control_male_1control_female_1
12	250 mg/dL	PAE_male_2PAE_female_2
2	0 mg/dL	control_male_3control_female_3
1	210 mg/dL	PAE_male_3PAE_female_3
8	0 mg/dL	control_male_4control_female_4

Assessed 30 min following ethanol vapor exposure.

At GD 14.5, we micro-dissected fetal mouse dorsal telencephalon from the litters of 6 separate dams. For each dam we combined two female pups to form one female sample and two male pups to form one male sample (see Figure S1 for schematic of experimental protocol). From these collected neuroepithelium tissues, we sequenced 139,195 single cells. The number of single cells sequenced ranged from 7,475 to 18,612 cells/sample. The mean number of reads per cell ranged from 21,662 to 60,579, and the number of genes detected per sample averaged 19,933 in all samples (complete sequencing metrics presented in Table S1). The percent of reads mapping to mitochondrial genes in all cells averaged 6.74%, within expected parameters for assessing live cells ((10xGenomics, 2019), Figure S2A).

Single-cell RNA sequencing captured the heterogeneity of the murine fetal developing cortex

A graph-based clustering algorithm grouped the 139,195 individual cells into 33 separate clusters (Figure 1A). We used previously published markers of cell identity (Loo et al., 2019) and differential gene expression to identify and annotate each of the 33 clusters. The clusters were identified as one of nine cell types: stem/progenitor cells, classified as VZ, VZ/SVZ, SVZ, transient progenitor cells; excitatory neurons, classified as lateral layer V-VI, rostral layer V-VI, layer I; interneurons; and non-neural cells, including endothelia, microglia, and erythrocytes. Five clusters expressed high levels of ventricular zone (VZ) markers Prom1 (Holmberg Olausson et al., 2014) and Pax6 (Englund et al., 2005). Overlapping Pax6, NeuroG2 (Johnson et al., 2015), and Eomes (Englund et al., 2005) expression marked the clusters at the junction between SVZ/VZ zones. As Sox2 expression was shown to decrease from VZ to SVZ to cortical plate (Vinci et al., 2016), we ordered the clusters within VZ, SVZ/VZ, and SVZ in order of decreasing mean Sox2 expression. Radial glia markers were expressed in VZ, VZ/SVZ, and some of SVZ clusters. We used markers summarized by Loo and colleagues (Loo et al., 2019) to annotate clusters belonging to layer I, interneurons, lateral and rostral layer V-VI, microglia, and endothelial cells (Figure 1B, percent of each cell type and are shown in Table 2, for Figure 1D, annotation of clusters based on enrichment of cell lineage markers are presented in Table S2). We used the combined expression of Slc4a1 (Reithmeier et al., 2016), Alas2 (Chiabrando et al., 2014), and hemoglobin genes (Kingsley et al., 2006) to identify erythrocytes.

Figure 1

scRNAseq reveals the heterogeneity of GD 14.5 murine developing cortex

(A) t-SNE plots of the sequenced cells showing the breakdown of the cells to 33 clusters and the mean log2 expression levels of markers of different cell types. Red and blue ellipses correspond to magnified microglia and endothelial cell populations, respectively.

(B) Heatmap showing expression of markers of mature and immature cell types in various cell types; black lines separate clusters, with the major cell type for each cluster denoted above; red lines separate markers of each major subtype. Each vertical line (yellow shades) represents the Z score of the individual cell expression level for a random subset of 800 cells within the cluster. ∗ denotes cluster 1, which was further subclustered.

(C) t-SNE plot showing the breakdown of cluster 1 into 8 subclusters, heatmap showing Z score normalized expression of subcluster gene markers (complete heatmap presented in Figure S2), and bar chart showing the percent of cells in each subcluster in the three phases of cell cycle.

(D) Bar chart shows the percent of all sequenced cells in each cluster (top panel) and the percent of cells in each of the three cell-cycle phases in each cluster (bottom panel). VZ: ventricular zone, SVZ: subventricular zone, TPC: transient progenitor cells. Clusters 26 and 29 are erythrocytes, and 30 and 31 are microglia and endothelial cells, respectively.

Table 2

Percents of cells in each cell lineage

Cell type	Percent of cells (%)
VZ	24.91
VZ/SVZ	8.08
SVZ	22.1
Transient progenitor cells	8.57
Lateral layer V–VI	17.87
Rostral layer V–VI	9.66
Layer I	1.29
Interneurons	5.73
Non-neural	1.79

For clusters that were not directly determinable, i.e., clusters #4 and #18, we utilized the FindMarkers feature in the Seurat package, which identifies differentially expressed genes between two groups of cells using a Wilcoxon Rank-Sum test. We identified upregulated genes in clusters #4 and #18, compared with their UMAP-algorithm-assigned neighboring clusters #13 and #20 (see Figure S2C). Gene ontology analysis of the upregulated genes showed that those differentially upregulated genes were associated with axonal and neuronal projection guidance.

Cluster #1 cells (6.5% of the total cell population) expressed markers of both immature and mature cell lineages. To further define this large, potentially heterogeneous group of cells, cluster #1 (Figure 1B, asterisk) was further subclustered into eight distinct sub-clusters (Figure 1C, sub-clusters 1.0–1.7). We used the significantly enriched genes in each cluster, identified by the FindMarkers tool in the Seurat pipeline, to identify functions of distinct sub-clusters in cluster #1 (Table S3, representation of enriched genes in each subcluster, Figures 1C and S2B). The sub-clusters within cluster #1 contain cells that exhibit identities of path-finding cells (sub-clusters 1.2 and 1.7), migrating interneuron-like cells (sub-cluster 1.5), the capacity to form neuronal projections and expressing markers of glutamatergic synapses (sub-cluster 1.4), and markers of cell division (sub-clusters 1.1 and 1.7) (Figures 1C, Table 3).

Table 3

Gene ontology of biological processes for sub-cluster-specific enriched genes

Subcluster	Key gene ontology biological processes for subcluster marker genes	Percent of cluster 1 cells
1.0	Positive regulation of cell morphogenesis involved in differentiation; synapse organization	19.16
1.1	Mitotic spindle assembly; ribonucleoprotein complex subunit organization	17.89
1.2	Regulation of axon extension involved in axon guidance; chemotaxis	15.98
1.3	Microtubule cytoskeleton organization; apoptosis	13.09
1.4	Regulation of synaptic transmission, glutamatergic; positive regulation of excitatory postsynaptic potential; positive regulation of filopodium assembly	12.36
1.5	Regulation of membrane potential; cerebral cortex GABAergic interneuron differentiation; synapse organization	10.68
1.6	Positive regulation of cytokinesis; regulation of mitotic cell cycle	8.52
1.7	Regulation of axon extension/axon guidance; regulation of migration/localization	2.33

Cell-cycle composition (proportion of cells in G1, S or G2/M phases) and pseudotime score distribution were used to further order each cell cluster within its major annotated cell-type. In the case of cluster #1, we also quantified the percent of cells in each cell-cycle phase in each of the sub-clusters. All VZ-annotated clusters were predominantly in S and G2/M phases, and 3 out of 6 SVZ clusters were also predominantly in S and G2/M phases. Two SVZ clusters, clusters #13 and #20, which shared a profile of cell identity markers, differed in their cell-cycle phase composition, with cluster #13 cells predominantly in G2/M phase, whereas cluster #20 was predominantly in S phase. In contrast, more mature, neuronal-lineage-committed cells, including those within clusters defined as layer V-VI neurons, or as interneurons, were predominantly in G1 phase, as expected for post-mitotic cells. Non-neural cell types, erythrocytes, microglia, and endothelial cells were composed of cells in all three phases of cell cycle (Figure 1D).

We also identified four clusters that shared the interneuron markers, Dlx1, Dlx2, and Sst. Two of the interneuron clusters, including the largest of the four clusters, expressed markers of medial-ganglionic-eminence-derived cortical interneurons in addition to higher levels of genes involved in the synthesis of gamma-aminobutyric acid, GAD1, GAD2. We used the t-test-based FindMarkers feature in the Seurat package to identify genes enriched in each of the interneuron clusters compared with the total interneuron population (Table S4), resulting in the confirmation that these are four separate interneuron populations present in the dorsal telencephalon at this time point.

We then utilized the UMAP algorithm and ordered the cells based on their progression across a UMAP-initialized trajectory. The UMAP algorithm resulted in one major contiguous cluster, for which a trajectory was constructed and pseudotime was computed, and two additional minor clusters that were not further assessed in this analysis. We then generated a pseudotime score for all the cells in the contiguous cluster based on their progression across the trajectory (Figure 2A), and the pseudotime score distribution in each graph-based cluster was used to order that cluster within each of the cell types contributing to the contiguous UMAP cluster (Figure 2B). The pattern of expression of Pax6, Neurog2, and Eomes and Tbr1, markers for VZ, SVZ, and neurons respectively (Englund et al., 2005; Loo et al., 2019), along the pseuodtime continuum indicated that progression across pseudotime reflects a maturation process. Gas1, a marker of proliferative radial glia (Telley et al., 2016), shows high expression in VZ cells that decreases in expression as SVZ cells progress through pseudotime. Plxna4, a gene crucial for neuronal radial migration (Chen et al., 2008), shows increased expression at the pseudotime score corresponding to SVZ-type cells that decreases with pseudotime progression toward neurons of lateral layer V-VI (Figure 2C). The two UMAP-identified minor, non-contiguous clusters were composed of two interneuron clusters, graph-based clusters #33 and #27, two transient progenitor cell clusters, graph-based clusters #4 and #18, and two SVZ clusters, graph-based clusters #13 and #20, suggesting sub-population heterogeneity with respect to developmental stage.

Figure 2

Pseudotime analysis orders the cells based on sequence of gene expression changes through neuronal maturation

(A) UMAP-initialized trajectories of sequenced cells, color coded by identified cell types. Gray line represents the structure of the trajectory graph, and red circle denotes the root of the trajectory.

(B) Ridge plot of pseudotime score distribution in each of the t-SNE identified clusters. Clusters in minor UMAP trajectories do not map to a pseudotime score.

(C) Expression level of Pax6 (VZ marker), Gas1 (radial glia marker), Neurog2 (VZ/SVZ marker), Eomes (SVZ marker), Tbr1 (neuronal marker), and Plxna4 (axonal marker) across the pseudotime score continuum. VZ: ventricular zone, SVZ: subventricular zone. Clusters 26 and 29 are erythrocytes, and 30 and 31 are microglia and endothelial cells, respectively.

PAE did not alter the profile of cell clusters derived from fetal telencephalic cells but did result in sex- and cluster-specific alterations in developmental trajectory and cell cycle

To assess whether an episode of PAE resulted in altered cell identity or developmental trajectory, we calculated the percent of cells from each sample in each of the 33 clusters. Neither PAE nor fetal sex affected the percent of cells comprising the clusters (two-way ANOVA dependent variable: number of cells in each sample, independent variables: exposure, sex; ANOVA summary tables for pseudotime scores of 33 clusters are in Table S5). We next asked if PAE or fetal sex altered the maturation trajectory, as represented by the distribution of pseudotime scores across the cells within each cluster. We found that eight clusters exhibited a significant interaction effect between fetal sex and PAE, four clusters exhibited a main effect of PAE, and two clusters exhibited a main effect of sex (Figure 3A; ANOVA summary tables for the pseudotime scores of 33 clusters are in Table S5). Of the eight clusters that showed an interaction effect, two VZ clusters, 1 SVZ cluster, and 2 lateral layer V-VI clusters exhibited PAE-induced leftward shifts in pseudotime score, i.e., toward diminished maturity, in males, but rightward shift, indicative of increased maturity, in female clusters. PAE exposure resulted in higher pseudotime scores in two VZ clusters, but lower pseudotime scores, in one lateral layer V-VI cluster (Figures 3B and S3).

Figure 3

Effects of PAE on cell trajectory and cell-cycle phase

(A) Significant main (Exp: exposure; sex) or interaction (Int) effects of sex and alcohol exposure (∗ p values <0.05, two-way ANOVA) on pseudotime scores.

(B) Density plots showing the distribution of pseudotime scores by cell number in select cell clusters comparing male (blue) and female (pink), control (solid line), and alcohol-exposed (dotted) conditions.

(C) Significant main or interaction effects of sex and ethanol exposure (∗ p values <0.05, two-way ANOVA) on cell-cycle phases across cell clusters.

(D) Percent of cells in each of the three cycle phases in select clusters, with males (blue) and females (red), control (light color), and ethanol (dark color) values shown. (€: significant main effect of exposure; $: significant main effect of sex, two-way ANOVA); data are represented as mean ± SEM. VZ: ventricular zone, SVZ: subventricular zone, TPC: transient progenitor cells. Clusters 26 and 29 are erythrocytes, and 30 and 31 are microglia and endothelial cells, respectively.

We next tested if PAE or fetal sex affected the percent of cells in each of the cell-cycle phases in each cluster (Figure 3A). Cluster #10 (lateral layer V-VI) exhibited a main effect of PAE on the percent of cells in G1 phase and main effects of both sex and PAE on the percent of cells in G2M phase. There was a main effect of PAE on the percent of cells in both G1 and S phase in cluster #19 (SVZ cells) and an interaction effect of sex and alcohol exposure on the percent of cells in G1 phase. Microglia (cluster #30) exhibited a main effect of PAE on the percent of cells in S phase (Figures 3C and S4A).

PAE effects on cell death potential in the developing cortex

We next investigated whether PAE resulted in differential expression of caspases, as an index of cell-death susceptibility (Thornberry and Lazebnik, 1998). Our data show that at GD 14.5, developing cortical cells in all clusters expressed Caspases 2, 3, and 6 (expression levels in various clusters are shown in Figure S4B), although there was within-cluster variation, with some cells expressing low-to-undetectable transcript numbers. Other caspases such as caspase 1 were expressed in virtually no cells. Out of 33 clusters, Caspase 2 was downregulated by PAE in one female cluster, and Caspase 3 was downregulated by PAE in two male clusters and three female clusters, whereas Caspase 6 was upregulated by PAE in one female cluster. These data collectively show that PAE does not have a significant and uniform effect on the potential for execution of apoptosis.

PAE results in distinct sex-specific patterns of gene expression dysregulation primarily in migrating cell populations

To investigate the effects of PAE on transcription profiles in males and females and to identify the baseline transcriptional expression differences between male and female developing cortices, we utilized the DEsingle package to identify differentially expressed genes. Ethanol exposure resulted in more differentially expressed genes in female clusters compared with male clusters. An SVZ cluster and lateral layer V-VI neuron cluster exhibited the highest numbers of alcohol-dysregulated genes in females (Figure 4A), with majority of genes downregulated with general differential expression (DEg-type, Figure S5A). Interestingly, hemoglobin genes, and Hbb-y, the embryonic variant of the beta hemoglobin-subunit, in particular, were among the few genes that were consistently upregulated in both males and females and at all stages of neural maturation, from VZ and SVZ cells, to mature neurons. Moreover, male-derived fetal neural cells also exhibited a developmental-stage-independent upregulation of additional hemoglobin genes, including Hba-a1 and Hba-a2 and Hbb-bs (Figures 4B and S7G–S7J). Overall, the highest number of dysregulated genes in male clusters were observed in rostral layer V-VI neurons. Few to no differentially expressed genes were observed in non-neural cell types in both males and females (box plots of magnitude of fold change of significant genes in each cluster are presented in Figure S7, complete results from DEsingle analyses are presented in Tables S6 and S7).

Figure 4

Prenatal alcohol exposure results in gene dysregulation in murine developing cortex in lineage and sex specific manner

(A) Bar charts depicting the number genes with dysregulated expression following ethanol exposure in each female (top panel) and male (bottom panel) cluster and the breakdown of the dysregulated genes into each of the three dysregulation subtypes.

(B) Bar charts depicting the number of female (left) and male (right) clusters (x axis) in which each of the genes (on y axis) was downregulated (top panels) or upregulated (bottom panels). Bars are color coded by major cell type (Gm21887 is part of an X-Y syntenic region containing Gm47283 (Erdr1) and considered an Erdr1 paralog (Merkenschlager et al., 2019)).

(C) Heatmap of IPA for biological pathways with significant enrichment (Benjamini-Hochberg corrected p value < 0.05) for PAE-dysregulated genes in female (top panel) and male (bottom panel) clusters. Significant pathways are marked by boxes; gray color indicates no known direction of activation; color scale indicates the Z score of activation. VZ: ventricular zone, SVZ: subventricular zone, TPC: transient progenitor cells. Clusters 26 and 29 are erythrocytes, and 30 and 31 are microglia and endothelial cells, respectively.

We next investigated whether alcohol exposure dysregulates the same genes in males and females. We calculated the percent of overlap between differentially expressed genes in each cluster PAE relative to controls, in males and females. Most clusters exhibited minimal overlap between DE genes in males and females (ranging from 0.0%–31.1%, average = 6.9% overlap). However, cluster #2, with a mixed VZ/radial glial identity, exhibited the highest percentage (31.1%) of overlap between males and females (Figure S5C).

We identified frequently dysregulated genes, which are downregulated or upregulated in >50% of the clusters (Figure 4B), and found that genes expressed in both males and females were altered in different directions by PAE. For example, ubiquitin-b (Ubb) was upregulated in 17 female clusters but downregulated in 21 male clusters (Figure S7B). Intriguingly, the female-specific X-chromosome-inactivation lncRNA, Xist, was significantly downregulated by PAE in 19 out of 33 female clusters (see below, Figure S7A).

To further understand the potential biological consequences for DE genes, we performed pathway-overrepresentation analysis (IPA) on DE genes in each cluster separately for male and female cells. We then compared biological pathways significantly enriched for DE gene expression to identify common dysregulated pathways across lineages and pathways that were dysregulated in a cell-type- or sex-specific manner. Our analyses showed that PAE resulted in a wide disruption in translation control (EIF2, eI4) and stress/nutrient (mTOR) pathways in both male and female cortices across multiple stages of neural maturation. Other pathways like 14-3-3 mediated signaling, epithelial adherens junction signaling, and phagosome maturation were identified as significant in similar clusters in males and females but were more dysregulated in male clusters (Figure 4C). Sex-specific responses to PAE were observed at discrete stages of maturation. For instance, although pathways associated with mitochondrial dysfunction and disruption in oxidative phosphorylation were identified in male rostral layer V-VI clusters, these same pathways were identified as disrupted in females in lateral layer V-VI clusters (Z scores and p values for all the pathways for each cluster are presented in Table S9).

Overall, the transcriptome of interneuron clusters that were present in this region of the developing cerebral cortex were among the least affected by alcohol exposure, in contrast to significant effects on interneurons observed in other PAE models (Marguet et al., 2020; Skorput et al., 2015). However, one interneuron cluster, cluster #11, was specifically affected by PAE in both males and females, but in opposite directions. IPA showed sexual dimorphism in the pathways enriched for DE genes, such that although EIF2, synaptogenesis, estrogen receptor signaling pathways, and oxidative phosphorylation were predicted to be activated in female cells in cluster #11, they were predicted to be inhibited in equivalent male cells (Figure 4C).

Using the upstream regulator analysis tool in IPA, we further identified candidate genes upstream of DE genes in each of the male and female clusters and genes identified as “significant upstream regulators” across multiple clusters (Table 4, Table S9). A subset of those upstream regulators (e.g., Bcl11a, Ctnnb1, Htt, see Figures S7C–S7E) was also dysregulated by PAE in some clusters. For instance, Bcl11a, a member of the BAF Swi/Snf chromatin-remodeling complex (Kadoch et al., 2013), which is itself an ethanol target in neural progenitors (Burrowes et al., 2017), was identified as a significant up-stream regulator of up to four dysregulated genes in 19 female and 21 male clusters and was itself downregulated in 10 female clusters and 3 male clusters, including both male and female cortical layer V-VI neurons. Htt, the Huntington gene transcript that facilitates radial cortical neuronal migration and survival development (Reiner et al., 2003; Tong et al., 2011), was a predicted upstream regulator in 21 male and female clusters with a predicted influence on up to 134 dysregulated genes in these clusters and was itself downregulated in 3 female clusters. Finally, β-Catenin, Ctnnb1, a canonical component of Wnt-signaling (MacDonald et al., 2009), was significantly downregulated in migrating SVZ and SVZ/radial glia clusters in both males and females. It is predicted to be a significant upstream regulator of up to 166 significantly dysregulated genes (mostly downregulated genes) by alcohol exposure in those clusters and was itself dysregulated in 12 male (upregulated in 3, downregulated in 9) and 8 female (upregulated and downregulated in 4 each) clusters.

Table 4

Significant upstream regulator of DE genes across multiple clusters

Female		Male
Transcriptional regulator	# of clusters	Transcriptional regulator	# of clusters
YAP1	23	BCL11A	21
CTNNB1	22	HTT	21
MYC	22	MYC	20
HTT	21	MYCN	20
NFKBIA	21	YAP1	20
YY1	21	CTNNB1	19
MLXIPL	20	MECP2	19
MYCN	20	PML	19
TP53	20	TP53	19
BCL11A	19

Next, we investigated whether DE genes were implicated in other neurodevelopmental disorders and may underlie common mechanisms of neurodevelopmental delay (NDD). We utilized the gene-scoring module in SFARI-gene (https://gene.sfari.org/database/human-gene/) and the dataset of genes implicated in autism spectrum disorders (ASD), ranked based on the strength of the evidence implicating that gene in the etiology of ASD (Abrahams et al., 2013; Larsen et al., 2016), and determined which genes were dysregulated by PAE. Intriguingly, although there is a large sex bias in presentation of ASD toward males (Werling, 2016; Werling and Geschwind, 2013), we found that a larger number of ASD-implicated genes were altered by PAE in female neural clusters compared with male, with up to 26% of ASD-identified NDD genes downregulated, mainly in VZ clusters and lateral layer V-VI neurons. Fewer NDD-related genes (up to 14%) were dysregulated in male clusters, in one VZ cluster, and the rostral layer V-VI neuron cluster (Figure S9).

Gene co-expression networks are dysregulated by alcohol exposure

Genes are typically not expressed individually, but instead in subsets, or “modules,” as has been previously documented in transcriptomic and functional studies (Kimbrough et al., 2020). Disruption in gene expression modules has been shown in neurodevelopmental and neurodegenerative disorders (Mariani et al., 2015; Parikshak et al., 2015; Tebbenkamp et al., 2014). Therefore, to investigate if PAE resulted in dysregulation of gene module expression, we constructed gene co-expression networks using weighted gene co-expression network analysis (WGCNA) and identified modules of highly interconnected genes as well as “hub genes” for each module and have previously been used in scRNAseq analyses (Cha and Lee, 2020; Luo et al., 2015; Lv et al., 2019; Mao et al., 2019). Co-expression networks were first constructed independently for the top most variant 20% of genes in each control male and female cluster and identified hub genes for each module. We next quantified the number of genes dysregulated by PAE in each of the modules and ranked the modules according to the proportion of dysregulated genes (Figure 5A, Table S10).

Figure 5

WGCNA shows that PAE disrupts gene networks

(A) Bar charts showing the top ethanol-dysregulated gene modules (y axis) in male (top panel) and female (bottom panel) clusters, color coded by cluster number, which is also indicated by _# in the module name. Hub-gene and number of genes constituting each module are shown, with arrows indicating direction of ethanol-induced expression change for each significantly altered hub-gene in the respective module. Ratio of module dysregulated genes relative to total number of genes expressed in the module is shown on x axis.

(B) Venn diagram of the overlap between Eef1a1 module genes in female and male cluster #2, female cluster #12, and male cluster #10.

(C) Histogram of the number of resampling-with-replacement iterations of male cluster #3, in which each gene was identified as an eigen (hub) gene. Dotted lines indicate 50% and 80% thresholds; hub genes identified in the original analysis are underlined in blue.

(D) Heatmap of Eef1a1 module gene (individual horizontal lines) regulation by ethanol exposure in female (top panel) and male (bottom panel) clusters. Color legend indicates the fold change in expression. VZ: ventricular zone, SVZ: subventricular zone, TPC: transient progenitor cells. Clusters 26 and 29 are erythrocytes, and 30 and 31 are microglia and endothelial cells, respectively.

Eef1a1, a critical gene for translational control, is a hub-gene dysregulated by PAE

Our data showed that for both male and female cells, the hub genes in four out of the top five dysregulated modules, were also dysregulated by PAE. For example, Eef1a1, the fetal neural-isoform of the translation elongation factor (Cho et al., 2012; Pan et al., 2004), was both downregulated by PAE (Figure S7F) and identified as a hub gene for two highly dysregulated modules in both males and females (female cluster #10 blue module, female cluster #12 blue module, male cluster #3 brown module, male cluster #10 blue module). These “Eef1a1 hub modules” share common gene members, but also include cluster- and module-specific genes (Figure 5B), and are collectively associated with ribonucleoprotein complex assembly and translation (Table 5, Table S11). Eef1a1 was shown to be a stable male cluster #3 hub gene in modules identified by bootstrap-WGCNA analyses (Figure 5C); in addition, the same analyses validated the stability of modules identified in male cluster #3 as shown by the high fraction of iterations with high similarity indices between the bootstrap-identified modules and each of the modules identified in the original analysis (Figure S8). The gene members of Eef1a1-containing modules exhibited a similar PAE-related pattern of downregulation in 2 VZ clusters (clusters #6 and #12), 1 SVZ cluster (cluster #3), a lateral layer V-VI (cluster #10), and a rostral layer V-VI (cluster #17) clusters in both males and females (Figure 5D).

Table 5

IPA function annotation for Eef1a1 combined modules genes (586 genes)

Function annotation	p value	# Molecules
Metabolism of protein	3.02 X 10−54	121
Synthesis of protein	2.25 X 10−68	98
Expression of RNA	170 X 10−11	95
Translation	3.10 X 10−82	94
Apoptosis	1.35 X 10−07	94

In contrast, several clusters exhibited sex-specific responses to PAE. For example, genes in Eef1a1-containing modules were generally downregulated in one male interneuron cluster (cluster #11), but few of those genes were downregulated in the same female cluster. In females, Eef1a1-associated genes were downregulated in response to PAE in 3 neuronal clusters belonging to lateral layer V-VI (clusters #23, #10, and #8, Figure 5D).

Female-specific responses to PAE mediated by Xist downregulation

Female clusters showed a higher number of dysregulated genes compared with male clusters. We sought to investigate possible female-specific mechanisms that explained this increased dysregulation. As expected, the X-inactivation long non-coding RNA (lncRNA) Xist was only expressed in female cells (Figure 6A). Xist was significantly downregulated by ethanol exposure in 19 out of 33 clusters (Figure 4B). Interestingly the Xist antagonist, Tsix, was upregulated by PAE in 7 clusters (Figure S7A). WGCNA also identified Xist as a hub for modules in one VZ cluster, one VZ/SVZ cluster, and one Layer I cluster (yellow_6 (Figure 6B), yellow_14, brown_25, respectively). Those modules exhibited 63, 13, and 8 downregulated genes out of 91, 54, and 207 genes, respectively. Bootstrap-WGCNA analyses validated the stability of Xist as a hub gene, and WGCNA identified modules in cluster #6, one of the PAE highly dysregulated clusters (Figures 6C and S8). Genes from the above three Xist-centered modules exhibited little overlap, compared with Eef1a1-associated modules (Figure 6D), suggesting that Xist-associated modules serve divergent functions in different cell types. We asked if Xist-associated genes were specifically downregulated in this cluster in females or if they were regulated by other mechanisms in the other female clusters or in male clusters. Xist modules, collectively, contained 352 genes, 146 of which were downregulated in female VZ cluster #6, 186 genes downregulated in female lateral layer V-VI cluster #10, 154 genes downregulated in female VZ cluster #12, and 107 genes downregulated in female rostral layer V-VI cluster #17. However, this pattern of dysregulation was not evident in any of the male clusters (Figure 6E).

Figure 6

Female-specific responses to PAE are mediated by ethanol-induced disruption of female-specific Xist gene expression

(A) t-SNE plot of sequenced female (left) and male (right) cells color coded by the log2 expression level of female-specific Xist (top) and male-specific Ddx3y (bottom).

(B) Topological overlap matrix of gene members of blue, brown, green, and yellow modules of cluster 6 (left), and zoomed yellow module genes dissimilarity matrix, Xist hub gene, and arrows showing key, highly connected genes. Light colors represent low topological overlap and progressively dark blue represents higher topological overlap.

(C) Histogram of the number of resampling-with-replacement iterations of female cluster #6, in which each gene was identified as hub gene. Dotted lines indicate 50% and 80% thresholds; hub genes identified in the original analysis are underlined in blue.

(D) Venn diagram shows the overlap between gene members of Xist module genes in female clusters 6 and 14 yellow modules and cluster 25 brown module.

(E) Heatmap showing regulation of Xist module genes (individual horizontal lines) by ethanol exposure in female (top panel) and male (bottom panel) clusters. Color legend indicates the fold change of expression.

(F) Location of X chromosome expressed genes in GD 14.5 murine developing cortices (first vertical panel, black horizontal lines); known X-inactivated genes (orange) and genes that escape X-inactivation (blue) (data derived from (Berletch et al., 2010; Carrel and Willard, 2005)) (second vertical panel); and X chromosome genes upregulated by ethanol exposure in female clusters (third and fourth vertical panels, false-discovery rate corrected p value threshold of 0.05 and 0.2, respectively, calculated with the DE Single analysis pipeline). Width of the line indicates the number of clusters in which the respective gene is upregulated. Genes upregulated in more than five clusters are labeled by gene name and number of clusters that show upregulation. VZ: ventricular zone, SVZ: subventricular zone, TPC: transient progenitor cells. Clusters 26 and 29 are erythrocytes, and 30 and 31 are microglia and endothelial cells, respectively.

Due to the high percentage of dysregulated genes in co-expression module yellow_6, within a VZ cluster, we asked if genes belonging to this specific Xist module in cluster #6 were also dysregulated in other clusters. A subset of the Xist module genes were downregulated in 29 clusters in females, whereas they were downregulated in only 7 male clusters (Figure S10A). Gene members of this Xist module, including Celf2 (downregulated in 18 female clusters), Luc7l2 and Tia2 (dysregulated in 17 female clusters), Ddx17 (downregulated in 12 female clusters), and Rbm25 (downregulated in 12 female clusters), are all regulators of spliceosome assembly and alternative mRNA splicing activity (Figure S10C).

RNA splicing was a significant and common functional annotation for genes that were enriched in the Xist modules (Table 6). We therefore further investigated the effect of PAE on the expression of spliceosome-activity-associated transcripts in male and female fetal cortical cells. We extracted the genes falling into each of the following gene ontologies: positive regulation of mRNA splicing via spliceosome, negative regulation of mRNA splicing via spliceosome, and spliceosomal complex assembly. Our results show that VZ-type cells (clusters #3, #6, and #12) and lateral layer V-VI-type cells (cluster #10) exhibited similar dysregulation of spliceosomal function genes in female cells, but this dysregulation was not observed in males (Figure S10A). In addition, a group of mRNA splice site selection genes, Celf2, Luc7l2, and Luc7l3, were downregulated in more than 10 female clusters but not in any male clusters. Luc7l3 specifically was upregulated in 12 male clusters (Figure S10B).

Table 6

IPA function annotation for Xist combined modules genes (352 genes)

Function annotation	p value	# Molecules
Transcription	2.65 × 10−19	101
Expression of RNA	687 × 10−19	104
Transcription of RNA	397 × 10−17	87
Transcription of DNA	188 × 10−16	77
Splicing of RNA	209 × 10−16	29
Organismal death	170 × 10−15	101
Activation of DNA endogenous promoter	135 × 10−14	64
Processing of RNA	273 × 10−14	33
Processing of mRNA	104 × 10−11	24
Perinatal death	140 × 10−11	36
Splicing of mRNA	604 × 10−11	21

Because X chromosome inactivation in females is maintained by the expression of Xist lncRNA (Marahrens et al., 1998), we expected that the downregulation of Xist lncRNA would result in the partial loss of X-inactivation and the upregulation of genes harbored within X chromosome inactivation regions. We surveyed the effects of ethanol exposure on the expression of the X chromosome genes in females and males, and found, as expected, a pattern of X-chromosome-associated genes upregulated in the female clusters compared with male clusters (Figure S11). As further evidence that this increase in X chromosome gene expression is due to decreased X-inactivation, we determined if the upregulated genes were located within known chromosomal regions of X-inactivation (Berletch et al., 2010; Carrel and Willard, 2005) and found that upregulated expression of X chromosome genes (specifically genes frequently upregulated across clusters) were indeed located within known X-inactivation regions (Figure 6F).

Cis-regulatory effects of Xist lncRNA on X chromosome inactivation are well established, although to our knowledge, an association between loss of Xist and transcript expression from autosome-encoded genes has not been documented. We therefore investigated the association between Xist expression and autosomal gene expression by assessing autosomal genes whose induction was significantly correlated with the loss of Xist in each of the 33 cell-type clusters. In each cluster we identified the chromosomal location of the genes showing significant negative correlation with Xist and quantified the number of significantly negatively correlated genes on each chromosome. Our results show a high number of autosomal genes (both raw number and normalized to the number of genes expressed from each respective chromosome) significantly negatively correlated with Xist in clusters #1 and #14 and to a lesser degree, in cluster #6. All three clusters were characterized by high percentage of cells in S and G2/M phases. Our analyses identified chromosomes 2, 4, 5, 7, and 11 as having the highest frequency of genes whose expression was inversely correlated with Xist expression (Figure S12). These data do not by themselves support a causal, trans-regulatory link between loss of Xist and upregulation of autosome-encoded genes but do point to an associative linkage that merits further investigation.

Discussion

In this study, we used scRNAseq transcriptomics analysis, which facilitates characterizing different cell types based on transcriptional profiles of individual cells (Loo et al., 2019; Macosko et al., 2015), to assess the impact of a single, temporally restricted episode of PAE at gestational day 12.5. In mice, this exposure period corresponds to the emergence of the subventricular zone (SVZ) and cerebral cortical pre-plate and the birthdate of the early neurons of the cortical plate, layer V/VI corticothalamic and corticospinal neurons (Caviness, 1982; Takahashi et al., 1995, 1999; Workman et al., 2013). We assessed cerebral cortical development after 48 h, i.e., at GD 14.5, which encompasses the time period for the completion of ∼50% of the founder-population cell cycles required for cortical neurogenesis (Takahashi et al., 1995) and for the migration of layer V/VI neurons to their appropriate position within the cortical plate (Workman et al., 2013). Our analysis showed that there is a high degree of cellular heterogeneity associated with each assessed cortical laminae and that this cell-population diversity was replicated across six separate pregnancies and in both male and female fetal cortices. For example, we documented 15 sub-populations of immature neural progenitors with identities related to the VZ and SVZ. The major VZ cluster, cluster 1, could itself be sub-divided into 8 additional sub-clusters, suggesting that at least 22 different neural precursor subtypes contribute to the neuroepithelium of the VZ and SVZ at this developmental stage. We also identified 10 distinct neuronal sub-populations with identities related to layers V, VI, and I (pre-plate origin), as well as four sub-populations of cortical interneurons.

Pseudotime analysis ordered most of the identified neural cell clusters along a temporal continuum, suggesting that these molecularly heterogeneous clusters represented stages along a developmental continuum. Consistent with this finding, developmental-stage-selective mRNA transcripts (Englund et al., 2005; Telley et al., 2016) exhibited discrete expression windows, with radial glial markers Pax6 and Gas1 preferentially expressed in cells with earlier pseudotime scores, neurogenic markers Neurog2 and Eomes expressed in cells with mid-timeline pseudotime scores, and post-mitotic neuronal markers Tbr1 and Plxna4 expressed in cells with later pseudotime scores. At GD 14.5, there were few cells identified as having endothelial or microglial identities and none with mature astrocyte or oligodendrocyte identity, an expected outcome, because the period of cortical neurogenesis precedes gliogenesis (Morrow et al., 2001). Cell-proliferative capacity was mainly relegated to the more immature cell types, although residual proliferative capacity was observed among two interneuron populations. VZ clusters differed in numbers of cells assigned to the stage of mitosis. For example, clusters 1.6, 21, and 6 contained cells predominantly in mitosis, whereas clusters 2 and 12 contained cells predominantly in S-phase, suggesting that these cell clusters were also spatially segregated by interkinetic motion into apical and basal compartments of the VZ and exposed to different microenvironments (McConnell, 1995; Taverna and Huttner, 2010) that may modulate the acute response to PAE. Interestingly, cluster 14, which was identified as an intermediate stage between VZ and SVZ, had an equivalent number of mitotic-type cells as VZ clusters 21 and 6 and is consistent with the presence of an asymmetrically dividing VZ fraction (Noctor et al., 2001), destined to exit the VZ toward the SVZ. Importantly, we did not observe any significant differences in cell population distribution due to fetal sex. Both male and female fetal cortices exhibited similar patterns of cellular heterogeneity.

PAE has been well documented to result in microcephaly in both human populations (Archibald et al., 2001; Norman et al., 2009; Roussotte et al., 2012) and in animal studies (Maier et al., 1999b; O'Leary-Moore et al., 2010). Consistent with our previous observations in ex vivo model systems (Prock and Miranda, 2007; Santillano et al., 2005), we did not find evidence for extensive cell death due to ethanol exposure. For instance, there was no loss within identified clusters, as the numbers of cells in each cluster were not altered due to PAE, nor was there a loss of cellular subtypes, as all identified clusters were present in all groups. Surprisingly, no cells expressed caspase-7, and some cell clusters with SVZ/transient progenitor identity also expressed low to undetectable levels of Caspases 3 and 6, which have all been defined as identified as and final common effector caspases (Thornberry and Lazebnik, 1998), suggesting overall developmental-stage-related resistance to apoptosis induction. Moreover, there were no changes in mitochondrial gene reads and either no change or decreased expression of caspase gene transcripts due to PAE, suggesting that apoptosis and the potential for cell death remained unaffected by exposure. However, we found evidence that PAE resulted in sex-independent increases in the proportion of cells in S-phase in one SVZ subpopulation, consistent with previously published reports (Miller, 1989; Miller and Nowakowski, 1991; Santillano et al., 2005). PAE also resulted in shifts in pseudotime scores across multiple immature and mature cell types, suggesting an impact on maturation. However, the effect of PAE on pseudotime shifts was also sex dependent. In SVZ cluster 13 and lateral layer V/VI neurons (clusters 22 and 23), for example, fetal male cells experienced a pseudotime delay, whereas females exhibited a pseudotime advance, suggesting that ethanol inhibited maturation in males but facilitated maturation in identical female cells. This conclusion is based on a significant statistical interaction effect between sex and exposure, which acknowledges the contribution of baseline differences between males and females. Moreover, an analysis of differentially expressed genes showed that PAE dysregulated more genes in female cells, compared with male cells, across developmental stages, and in general, a majority of genes in both sexes were downregulated rather than upregulated by PAE. In this study, the transcriptome of microglia, endothelial cells, and erythrocytes were relatively un-responsive to ethanol exposure, although microglia did exhibit a significant increase in numbers of cells assigned to S-phase. One possible explanation is that the smaller number of non-neural cell types obtained in our analysis may have resulted in decreased statistical power to discover small effect-size changes due to an episode of PAE. Although no genes in microglia and endothelial cells exhibited dysregulation that exceeded the false-discovery rate, DESingle analyses did identify candidate dysregulated genes in those cell types using more relaxed criteria (p < 0.05 and 0.005 respectively, Figure S5B), suggesting at least a minimal level of responsiveness to this single episode of PAE. In addition, the published literature documents that more prolonged ethanol exposures, particularly at later stages of neural development, do alter the behavior of non-neural cells (Boschen et al., 2016; Drew et al., 2015; Girault et al., 2017; Jégou et al., 2012), suggesting that the timing and pattern of exposure is also a critical determinant of non-neural cell adaptations to PAE.

Despite the discrepancy in numbers of dysregulated genes between males and females, pathway overrepresentation analysis showed that male and female cells exhibited some common vulnerable pathways related to inhibition of translation initiation, suggesting that sex differences in gene responses to PAE could nevertheless result in some common, sex-independent functional outcomes. WGCNA analysis with bootstrap resampling further emphasized the importance of inhibition of translation as a key response to PAE, because translation initiation genes were hub genes in a number of gene modules in both males and females. Our data identified Eef1a1, the eukaryotic translation elongation factor, as hub gene in highly dysregulated modules within four cell clusters, and expression of Eef1a1 itself was downregulated by ethanol exposure. Moreover, those Eef1a1-associated modules contained genes associated with ribosomal biogenesis, translational processes, and apoptosis. Eef1a1-associated genes exhibited a similar pattern of dysregulation in migrating cells and rostral layer V/VI neurons in both males and females but sexually dimorphic pattern of dysregulation in SVZ/VZ progenitors, lateral layer V/VI neurons, and interneurons. In addition, genes differentially expressed due to PAE in both male and female cortices were significantly enriched for pathways related to eukaryotic initiation factor 2 (eIF2) and mTor (mammalian target of rapamycin, an inhibitor of translation (Fournier et al., 2013)) in most of the clusters. These data are consistent with previous reports showing that that ribosomal biogenesis and disruption in translational processes are targets of PAE (Berres et al., 2017). Perturbations in global protein synthesis and dysregulation of mTor have also been implicated in ASD and other neurodevelopmental disorders (reviewed in (Mohammad et al., 2019)) and in adult alcohol use disorders (Morisot et al., 2019), suggesting that perturbation of translation control is not only an important component of a number of disorders of development but that mTOR may also be commonly altered by alcohol exposure.

WGCNA by focusing analysis on highly intercorrelated genes rather than on single genes supports the discovery of gene networks and their eigen (hub) genes, and therefore, potential mechanistic contexts for the effects of a developmental perturbagen like ethanol. These analyses, when applied to scRNAseq data (Cha and Lee, 2020; Luo et al., 2015; Lv et al., 2019; Mao et al., 2019), have the potential to identify cell-type specific mechanisms of vulnerability, although there are certainly challenges with implementing network modeling approaches for single-cell data, related primarily to the quality of the input data (Cha and Lee, 2020), that warrant caution. It is important to note that in this study each male/female fetal sample pair was obtained from a separate pregnancy, and in the case of the ethanol exposure paradigm, a range of maternal blood alcohol content was obtained despite stringently controlling exposure conditions, suggesting inter-animal differences in ethanol metabolism. Moreover, PAE effects are likely to be dependent on the developmental state at exposure and assessment time points, which may be moderated not only by fetal sex (DiPietro and Voegtline, 2017) but also by other intrauterine sources of between-fetus variation (Ryan and Vandenbergh, 2002). All these factors are expected to maximize between-sample variance. However, despite these sources of variance, bootstrap resampling analyses suggest that gene networks identified by WCGNA were generally stable, suggesting that the networks represent true cell-type-specific transcription states and not between-sample variance (Figures 5C, 6C, and S8).

A key finding in our exploration of the basis of sex differences in the PAE response was our identification that the X chromosome inactivation factor Xist was a hub gene in a number of gene co-expression modules and was itself downregulated in a number of immature and mature cortical cell types. Co-incident with loss of Xist, we found that the antisense lncRNA Tsix, which serves as an inhibitor of Xist (Ohhata et al., 2008; Sado et al., 2005), was upregulated along with a number of X chromosome genes within the X-inactivation region, indicating decreased Xist leads to decreased X-inactivation. These data suggest that in females, loss of X chromosome inactivation represents an important and sex-specific response to PAE. We found that gene members of female, Xist hub-containing modules, which were also downregulated by PAE, were preferentially related to RNA splicing. Cortical development, lamination, and cell fate are influenced by alternative splicing events, and important neural transcription factors exhibit maturation-stage specific isoforms (reviewed in (Porter et al., 2018)). Spliceosome deficits have also been identified in neurodevelopmental (Sanders et al., 2020; Smith and Sadee, 2011) and craniofacial disorders (Lehalle et al., 2015). These outcomes are consistent with some previous transcriptomic studies in both murine (Alberry et al., 2020) and avian (Berres et al., 2017; Garic et al., 2014) models of PAE that have also implicated spliceosome function as an important developmental target of ethanol. The Xist locus appears to be particularly vulnerable to developmental pathology. Other environmental exposures during development, for example exposure to bisphenol-A, have been shown to result in loss of X chromosome inactivation, which may mediate the neurodevelopmental effects of these chemical teratogens (Kumamoto and Oshio, 2013). X chromosome skewness (non-random X chromosome inactivation) in the absence of chromosomal mutations has been reported in females with autism spectrum disorders (Talebizadeh et al., 2005), X-linked mental retardation (Plenge et al., 2002), and Rett syndrome (Knudsen et al., 2006). Therefore, loss of X chromosome inactivation may contribute to etiology of a number of neurodevelopmental disorders. Moreover, although Xist is known to act in cis to inactivate the X chromosome, our data show that loss of Xist is also associated with the induction of autosomal-encoded genes across a number of chromosomes, suggesting a link between Xist and autosome activity. However, these data document associations and do not by themselves suggest causality. Further studies will be needed to identify any potential mechanistic linkage between loss of X-inactivation and the induction of autosomal genes. However, our data do indicate that this associative linkage is remarkably cell type specific and appears to be restricted to VZ and SVZ progenitors, suggesting that these developing cell stages may be particularly vulnerable to widespread dysregulation of gene expression across multiple chromosomes, along with loss of Xist.

Another key finding was that hemoglobin genes, embryonic beta-like subunit (Hbb-y) in females and multiple alpha (Hba-a1, Hba-a2), beta (Hbb-bs), and embryonic beta-like (Hbb-y) subunits in males, were one of the few gene classes that were consistently and significantly upregulated by PAE at multiple stages of neural maturation. Hemoglobin gene expression has previously been documented in neurons (Richter et al., 2009) and may serve as an adaptive neural response to stress. Importantly, beta-hemoglobin subunits have been shown to be present within the nucleus of neural cells, associated with nuclear histones and to increase histone H3 trimethylation (H3K4me3 (Brown et al., 2016)), an epigenetic activation mark, that, along with histone methyltransferases, has also been documented to be increased following PAE in rodent models (Chater-Diehl et al., 2017; Schaffner et al., 2020). This evidence raises the intriguing possibility that neural expression of hemoglobin genes, and beta globin genes in particular, may mediate some of the epigenetic effects of PAE. Moreover, Bcl11a, a member of the Swi/Snf chromatin remodeling complex (Kadoch et al., 2013) and inhibitor of embryonic and fetal hemoglobin genes (Guda et al., 2015; Sankaran et al., 2008; Xu et al., 2013), was itself inhibited by PAE, in both males and females, in layer V/VI neurons as well as in other neural clusters. Bcl11a haploinsufficiency in mice results in microcephaly, impaired cognitive, and social behaviors (Dias et al., 2016), and in humans, microdeletions within the BCL11A gene locus result in both retention of fetal hemoglobin and neurological symptoms such as developmental delay and autism (Basak et al., 2015; Funnell et al., 2015). Interestingly, Bcl11a was identified as a member of a class of upstream genes in dysregulated gene networks, including Htt and Ctnb1, which are also linked to developmental disorders. For instance, Htt gene variants have been implicated in Lopes-Maciel-Rodan syndrome, (LOMARS) a neurodevelopmental disorder with characteristics that overlap with Rett syndrome (Lopes et al., 2016; Rodan et al., 2016), and Ctnnb1 mutations have been linked to intellectual disabilities, craniofacial dysmorphology, and microcephaly (Dubruc et al., 2014; Winczewska-Wiktor et al., 2016), features that overlap with FASD. These data suggest that diverse congenital neurodevelopmental disabilities may share, at least in part, a common etiology. This hypothesis is further supported by our observation that PAE resulted in persistent downregulation of up to 26% of ASD-linked genes in females and up to 14% of ASD genes in males in both VZ cells and layer V/VI neurons.

In summary, our data show that a prevalent developmental perturbagen, alcohol, administered for a brief interval during the initiation of fetal cortical neurogenesis, has profound and persistent consequences on gene expression patterns, particularly for genes controlling RNA translation in progenitor cells and daughter neurons that were born during the exposure period. The effects documented in this study are the results of a single bout of ethanol exposure, intended to model a single episode of binge drinking early in pregnancy. The observed PAE effects on trajectory of cellular development and gene expression are likely to be dependent on the developmental state at exposure and assessment time points, which may be moderated by fetal sex (DiPietro and Voegtline, 2017) and other intrauterine sources of between-fetus variation (Ryan and Vandenbergh, 2002). Our analyses show that the effects of PAE were clearly sex dependent, with female cortical cells exhibiting an exaggerated transcriptomic response compared with male cells, which was linked to decreased X chromosome inactivation. The loss of Xist was linked to cis de-repression of X chromosome genes and was associated with repression of RNA splicing machinery. Finally, we identified an elevated neural hemoglobin gene response including elevated embryonic hemoglobin gene expression related to loss of the epigenetic repressive factor, Bcl11a. Moreover, Bcl11a was a member of a cluster of upstream regulators such as Htt and Ctnnb that mediate NDDs and were themselves dysregulated by PAE. These data suggest that the fetal response to PAE shares similarities with other congenital disorders of neural development.

Limitations of the study

The effects of prenatal ethanol exposure are well documented to be dependent on dose, timing, duration, and pattern of maternal ethanol exposure (Livy et al., 2003; Maier et al., 1999a; Maier and West, 2001). The current study modeled a single binge-like episode of maternal ethanol exposure, albeit at a critical fetal brain developmental stage of neurogenesis. This approach permitted an analysis of the heterogeneity in fetal cell-type response to a single perturbagen. Nevertheless, this study cannot capture the complexity associated with the effects of repeated and varied alcohol consumption patterns that is common in human populations. Additional scRNAseq studies at other temporal windows and using alternate exposure models will be needed to develop a complete picture of cell-type-specific vulnerabilities and resilience. Moreover, the effects of an exposure episode are likely to evolve over time, and the transcriptome response of individual cell types, layer V/VI neurons for example, are likely to evolve over developmental time scales, whereas other cell types, withing the cortical VZ for example, may disappear. These analyses also present evidence for a link between PAE and gene expression, for example, upregulation of neural hemoglobin genes, downregulation of a large complement of ASD-linked genes, or loss of X-inactivation in female fetal cells. However, at this time, these linkages between PAE and gene expression are associative. Further research will be needed to define mechanisms, possibly epigenetic pathways (Burrowes et al., 2017; Schaffner et al., 2020; Veazey et al., 2013, 2015), that link exposure to outcome. Further research will also be needed to determine both immediate- and long-term consequences of cell-type-specific perturbations in the transcriptome. Finally, even though this study analyzed a large number of cells, these ∼140,000 cells were derived from just six separate pregnancies, three control and three PAE. The analysis of scRNAseq data is rapidly evolving, and it is likely that future studies will consider the number of biological replicates insufficient to detect smaller but biologically significant effect sizes. We expect that the data from this study, deposited with NCBI (GEO: GSE158747), will be re-analyzed in combination with new data acquired by the research community, to facilitate more sophisticated insights into neuroteratology.

Resource availability

Lead contact

Rajesh C. Miranda, PhD. rmiranda@tamu.edu.

Materials availability

Not applicable.

Data and code availability

Raw and processed data files are deposited in NCBI GEO under accession number GSE158747.

Methods

All methods can be found in the accompanying transparent methods supplemental file.