Sequential perturbations to mouse corticogenesis following in utero maternal immune activation

1. CesarPCanales1
2. MykaLEstes1
3. KarolCichewicz1
4. KartikAngara2
5. JohnPaulAboubechara1
6. ScottCameron1
7. KathrynPrendergast1
8. LindaSu-Feher1
9. IvaZdilar1
10. EllieJKreun1
11. EmmaCConnolly1
12. JinMyeongSeo1
13. JackBGoon1
14. KathleenFarrelly1
15. TylerWStradleigh1
16. DeborahvanderList1
17. LoriHaapanen3
18. JudyVandeWater3
19. DanielVogt2
20. AKimberleyMcAllister[email protected]1
21. AlexSNord[email protected]1

1. Center for Neuroscience, UC Davis
   DavisUnited States
2. Department of Pediatrics & Human Development, Michigan State University
   East LansingUnited States
3. Division of Rheumatology, Allergy and Clinical Immunology, UC Davis
   DavisUnited States

• Research Article

• Genetics and Genomics
• Neuroscience

• MIA
• RNA-Seq
• transcriptomic
• brain development
• WGCNA
• Mouse

• publisher-id60100
• doi10.7554/eLife.60100
• elocation-ide60100

Introduction

Epidemiological association between maternal infection and neurodevelopmental disorders (NDDs) has been found for autism spectrum disorder (ASD), schizophrenia (SZ), bipolar disorder (BPD), anxiety, and major depressive disorder (MDD) 10Brown201126Estes and McAllister201561Parboosing et al.201362Patterson2009. Indeed, maternal immune activation (MIA) itself is sufficient to produce NDD-relevant outcomes in offspring in animal models 12Canetta et al.201447Machado et al.2015. Offspring from female mice exposed to the toll-like receptor 3 (TLR3) agonist and viral mimic polyinosinic:polycytidylic acid [poly(I:C)], and offspring in maternal immune activation models using the bacterial mimic lipopolysaccharide (LPS) 17Depino201559O'Loughlin et al.201771Simões et al.2018, recapitulate neuropathologies and aberrant behaviors seen in offspring born to flu-infected dams 69Shi et al.200372Smith et al.2007. MIA in mice produces behavioral and cognitive abnormalities with relevance across multiple NDD diagnostic boundaries 27Estes and McAllister201649Malkova et al.201264Richetto et al.2014.

MIA-induced outcomes in adult offspring have been well characterized; however, mechanisms of initiation and progression of NDD-related brain pathology remain poorly defined. Poly(I:C) increases proinflammatory cytokines such as IL-1β, IL-6, CXCL1, and TNF-α in the maternal circulation 6Ballendine et al.201516Dahlgren et al.200632Hsiao and Patterson201153Meyer et al.200652Meyer et al.200672Smith et al.2007. Downstream pathological changes in the brains of offspring have been reported across cell types, from neurons to microglia, and these changes have been linked to perturbations in proliferation, migration, and maturation in the developing brain 29Haida et al.201970Shin Yim et al.201773Soumiya et al.201189Zhang et al.201991Zhao et al.2019. Recent studies have identified acute transcriptomic changes in fetal brain in mice and rats at single time-points (3–4 hr) following MIA, with significant overlap between transcriptome changes in the cortices of MIA-exposed offspring and those altered in the brains of children with ASD 46Lombardo et al.2018. Although MIA could cause a dynamic sequence of transcriptional changes in the fetal brain following acute exposure, no integrated picture of MIA-associated transcriptomic pathology exists that maps changes spanning from the initial acute response to MIA through birth. Here, we combined time-course transcriptomics with neuroanatomy to map changes in prenatal mouse cortical development following MIA induced by poly(I:C) injection at E12.5.

Results

We used paired transcriptomic and neuroanatomical analyses to map changes in the fetal brain following MIA in a poly(I:C) mouse model (Figure 1, Figure 1—figure supplement 1, Figure 1—figure supplement 2 and Figure 1—figure supplement 3). Induction of MIA via poly(I:C) injection at E12.5 produces relevant phenotypes in mice 14Choi et al.201670Shin Yim et al.201772Smith et al.2007 and the specific implementation of this model was recently validated in our hands to produce sufficient levels of MIA to mimic viral infection and cause aberrant behavioral outcomes in offspring 25Estes et al.2020. The dosage of 30 mg/ml poly(I:C) elicited substantial elevations in maternal serum IL-6 (on average 2200 pg/ml) 4 hr following injection (Figure 1—figure supplement 1). Pregnant mice were injected with saline or poly(I:C) at E12.5, and dorsal telencephalon (E12.5 + 6 hr, E14.5, E17.5, and birth (P0)) was microdissected (Figure 1a). RNA-seq datasets were generated from male and female embryos from 28 independent litters across control and MIA groups, typically with one to three embryos represented per independent litter (Supplementary file 20, Supplementary file 5, File source data 1). For control and MIA, 7 vs. 7 samples at E12.5, 13 vs. 13 samples at E14.5, 6 vs. 6 samples at E17.5, and 10 vs. 12 samples at P0, were compared in a differential expression (DE) analysis. For transcriptomic analysis, we used a strategy of first defining DE genes at each time-point using the general linear model (GLM) approach implemented in edgeR, followed by mapping DE signatures to systems-level expression patterns via module assignment using weighted gene co-expression network analysis (WGCNA).

Figure 1.

# Global R markdown code chunk options
knitr::opts_chunk$set(message=FALSE, 
                      warning = FALSE, 
                      error=FALSE, 
                      echo=TRUE, 
                      cache=TRUE, 
                      fig.width = 7, fig.height = 6, 
                      fig.align = 'left')

# Installs packages if missing

list.of.packages <- c("sva", "RUVnormalize", "pheatmap", "edgeR", "GenomicFeatures", "mclust", "parallel", "ggplot2", "reshape2", "dplyr", "data.table", "DT", "RColorBrewer", "knitr", "ggrepel", "gridExtra", "grid", "plyr", "plotly", "Hmisc")

new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)


lapply(list.of.packages, require, character.only = TRUE)



exp.data <- read.csv("./Extra_submission files/Gene_counts.csv")
rownames(exp.data) <- exp.data$X
exp.data$X <- NULL

sample.info <- read.csv("./Extra_submission files/metadata.csv")

load("./Extra_submission files/exonic.gene.sizes")

# Exonic gene sizes were calculated as follows:
# Generate mm9 list of exon sizes 
# Method description at: https://www.biostars.org/p/83901/

#txdb <- makeTxDbFromGFF("G:/Shared drives/Nord Lab - Personal #Folders/Karol/Genomes/Mus_musculus/UCSC/mm9/Annotation/genes.gtf", #format="gtf")
#exons.list.per.gene <- exonsBy(txdb,by="gene")

# I parallelized this increasing the speed of the process >2x on a 4-core (logical) machine. 
# Use mclapply instead of parLapply if you use a Mac.
#cl <- makeCluster(detectCores())
#exonic.gene.sizes <- parallel::parLapply(cl, #exons.list.per.gene,function(x){sum(width(reduce(x)))})
#stopCluster(cl)

sample.info_submission <- read.csv("Extra_submission files/Supplementary Table 5, RNA-seq sample metadata_03_02_2021.csv")

#datatable(sample.info_submission, 
#          rownames = FALSE, 
#          filter = list(position = 'top', clear = FALSE, plain = FALSE), 
#          options = list(paging = FALSE, scrollX=T, scrollY="300px", searching = FALSE))


# Qualifies F or M based on the Xist expression
sex.by.rna <- c(ifelse(exp.data["Xist",]>1000,"F","M")) #

Xist_exp <- as.data.frame(reshape2::melt(exp.data["Xist",]))
Xist_exp <- cbind(Xist_exp, sex.by.rna)
Xist_exp <- arrange(Xist_exp, value)
colnames(Xist_exp) <- c("variable", "counts", "sex.by.rna")

#ggplot(Xist_exp, aes(x=sex.by.rna, y=counts, colours=sex.by.rna))+
#  geom_jitter()+
#  geom_boxplot(alpha=0.2)+
#  theme_bw()+
#  labs(title="Xist read counts", x = "", y = "Read counts")+
#  theme(plot.title = element_text(size = rel(2), hjust=0.5))+
#  scale_x_discrete(breaks= c("F", "M"), labels=c("Females", "Males"))


# Calculates RPKM values
# Gene lengths calculated with lapply
gene.lengths <- as.numeric(lapply(1:nrow(exp.data), function(x) FUN= as.numeric(exonic.gene.sizes[rownames(exp.data)[x]])))

rpkm.data <- rpkm(exp.data, gene.length=gene.lengths, log=T, prior.count=.25)
rpkm.data_linear <- rpkm(exp.data, gene.length=gene.lengths, log=F)
rpkm.data_linear_all <- rpkm.data_linear

#write.csv(rpkm.data, file = "GEO_submission #files/rpkm_log2_24015.csv")
#write.csv(rpkm.data_linear, file = "GEO_submission #files/rpkm_linear_24015.csv")


### Removes genes with low expression ###

# Lets plot a histogram of all rpkm values in a dataset.
#dim(rpkm.data) #24015 rows and 74 columns

df <- reshape2::melt(rpkm.data)
df$rpkm_linear <- 2^df$value # This converts the log2 RPKM values back to linear scale
colnames(df) <- c("gene_name", "sample", "rpkm_log2","rpkm_linear")

# Histogram of all rpkm values in a dataset.
#ggplot(df, aes(x=rpkm_log2))+ 
#  geom_histogram(bins = 1000, color="black")+
#  theme_bw()+
#  labs(title="Unfiltered dataset of 24015 genes", 
#       x="log2 RPKM", y = "Read counts")+
#  theme(plot.title = element_text(hjust = 0.5))

# Setting the threshold for minimum rpkm value in at least 2 samples. 
threshold = -2

# Number of genes after filtering
#sum(reshape2::melt(rowSums(rpkm.data > threshold) >= 2)$value) # 17051

# Let's filter the dataset and plot the histogram distribution
keep <- as.data.frame(reshape2::melt(rowSums(rpkm.data > threshold) >= 2))
keep$gene_name <- rownames(keep)
keep <- filter(keep, value == "TRUE")$gene_name

# NOTE: The set of included genes reported in the manuscript was determined using expression criteria from an expanded set of samples including a condition that was not used for the manuscript. The difference between the gene counts is 17195 in the manuscript from the larger set of samples, compared to 17051 genes that would pass in the set of samples used in the MIA vs. control comparison.  Note that there is no difference among DE genes passing FDR < 0.05 using either gene set.  

original_exp.data <- read.csv("Extra_submission files/Count_gene_set_from_initial_submission.csv")

# length(keep) # 17051
# length(original_exp.data$gene_name) # 17195

keep <- original_exp.data$gene_name

# 
datExpr_test <- as.data.frame(rpkm.data)
datExpr_test$gene_name <- rownames(datExpr_test)

# Filters the "keep" genes 
datExpr_test <- filter(datExpr_test, gene_name %in% keep)

# Generates a matrix-type data frame
rownames(datExpr_test) <- datExpr_test$gene_name
datExpr_test$gene_name <- NULL

# Overwrites the original datExpr object.
datExpr <- datExpr_test


# Plots the histogram after filtering
df <- reshape2::melt(datExpr)
df$rpkm_linear <- 2^df$value # This converts the log2 RPKM values back to linear scale
colnames(df) <- c("sample", "rpkm_log2","rpkm_linear")

#ggplot(df, aes(x=rpkm_log2))+
#  geom_histogram(bins = 1000, color="black")+
#  theme_bw()+
#  labs(title = "Dataset of 17195 genes \n after excluding genes expressed at a very low level", x="log2 #RPKM", y = "Read counts")+
#  theme(plot.title = element_text(hjust = 0.5))


# Removing low expression genes. Overwrites exp.data object used for DE analysis ### 
#dim(exp.data)  #exp.data contains counts, datExpr contains rpkms 

exp.data$gene_name <- rownames(exp.data)
exp.data <- filter(exp.data, gene_name %in% rownames(datExpr))

rownames(exp.data) <- exp.data$gene_name
exp.data$gene_name <- NULL 


# Recalculates rpkm values
# Gene lengths calculated with lapply
gene.lengths <- as.numeric(lapply(1:nrow(exp.data), function(x) FUN= as.numeric(exonic.gene.sizes[rownames(exp.data)[x]])))

# RPKM calculation
rpkm.data <- rpkm(exp.data, gene.length=gene.lengths, log=T, prior.count=.25)
rpkm.data_linear <- rpkm(exp.data, gene.length=gene.lengths, log=F)


#write.csv(rpkm.data, file = "GEO_submission #files/rpkm_log2_17051.csv")
#write.csv(rpkm.data_linear, file = "GEO_submission #files/rpkm_linear_17051.csv")



group <- ifelse(sample.info[,"Condition"]=="Saline",1,2)
rRNA <- as.numeric(sample.info$rRNA)
sex.by.rna <- factor(sample.info$sex.by.rna)
dpc <- factor(sample.info$DPC)
response <- factor(sample.info$Response)
lane <- factor(sample.info$Lane)

############################################

min.cpm.criteria <- 0.1 # really relaxed
test.data <- exp.data
#colnames(test.data) <- paste(sample.info$ExperimentalDesign, 1:102, sep = "_")
test.samples <- 1:nrow(sample.info)
min.cpm <- min.cpm.criteria
y <- DGEList(counts=test.data, group=ifelse(sample.info[,"Condition"]=="Saline",1,2))
keep <- rowSums(cpm(y)>min.cpm) >=2 #keeps only genes expressed in above min.cpm in at least 2 libraries in each group 
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y) # Normalizes for RNA composition
y <- estimateCommonDisp(y) # Estimates common dispersions. Calculates pseudo-counts, a type of normalized counts. Don't misteke them with 0+x type counts from other packages. Also "users are advised not to interpret the psuedo-counts as general-purpose normalized counts".
y <- estimateTagwiseDisp(y) #Estimates dispersions. Applicable only to experiments with single factor.
#Alternatively the two commands from above can be replaced with y <- estimateDisp(y)


#PCA plot

global.pseudo <- y$pseudo.counts
rownames(global.pseudo) <- rownames(y$counts) # this doesn't feel necessary because all(rownames(global.pseudo) == rownames(y$counts))
pca.results <- prcomp(scale(log(global.pseudo+1), center=T, scale=T)) # It is a good idea to scale your variables. Otherwise the magnitude to certain variables dominates the associations between the variables in the sample.

b <- data.frame(Samples=rownames(pca.results$rotation), 
                DPC=as.character(sample.info$DPC), 
                ExperimentalDesign=sample.info$ExperimentalDesign, 
                Condition=sample.info$Condition,
                sex.by.rna=sample.info$sex.by.rna,
                rRNA=sample.info$rRNA,
                lane=as.factor(sample.info$Lane),
                pca.results$rotation)
rownames(b) <- NULL

b$DPC <- ifelse(b$DPC == 19.5, "P0", b$DPC)

#PCA plot
PCA_plot <- 
  ggplot(b, aes(x=PC1, y=PC2, fill = DPC, colour=DPC, shape=Condition))+
  geom_point(size=3, alpha=0.6, aes(shape=Condition))+
  theme_bw()+
  labs(title = "PCA plot")+
  theme(plot.title = element_text(size = rel(1.5), hjust=0.5))

#PCA_plot


#PCA plot lacking lane 12 and the inhibitor samples - slightly different version
PCA_plot_sex <- ggplot(b, aes(x=PC1, y=PC2, colour=sex.by.rna, shape=Condition))+
  geom_point(size=3, alpha=0.6)+
  theme_bw()+
  labs(title = "PCA plot")+
  theme(plot.title = element_text(size = rel(1.5), hjust=0.5))
  

#PCA_plot_sex

#source("Dimensionality reduction plots.r")

grid.arrange(PCA_plot, PCA_plot_sex, ncol = 2)

Differential expression across embryonic brain development following maternal immune activation via poly(I:C) injection at E12.5.

(a) Schematic representation of MIA model and study design. E: embryonic day, P: postnatal day, IHC: immunohistochemistry, WB: western blot. Poly(I:C) was injected at E12.5. Samples for RNA-seq were collected at E12.5 + 6 hr, E14.5, E17.5, and at birth (P0). IHC and WB analysis were conducted on a separate animal cohort at E17.5 and P0. (b, c) Principal component analysis (PCA) of RNA-seq data indicates that developmental age accounts for the majority of variance across samples. Age (DPC) or sex are represented as colored symbols in (b) and (c) respectively, and poly(I:C) or saline treatment indicated by circles or triangles, respectively (d) Heatmap representing relative gene expression changes between control and MIA samples across time-points. Hierarchical clustering by relative fold changes shows stage-specific differential expression signatures (DE genes shown have FDR < 0.05 and log₂ fold change (log₂FC)>1 or < −1). (e) Volcano plots of DE effect size and significance show stage-specific differences in number of DE genes and DE effect size, with the strongest dysregulation at E17.5. Colors represent directionality and statistical significance. (f) Numbers of upregulated and downregulated DE genes at FDR < 0.05 and p<0.05 thresholds again show varying DE genes numbers with a peak of dysregulated genes at E17.5. (g) Intersection of stage-specific DE genes with the 82 SFARI autism-associated mouse gene orthologs that were expressed at measurable levels in the RNA-seq data. E17.5 DE genes (FDR < 0.05) were enriched for SFARI genes (p=3.9e-04, hypergeometric test). Concentric circles represent developmental time-points. Light red, upregulated DE genes (p<0.05); dark red, upregulated DE genes (FDR < 0.05); light blue, downregulated DE genes (p<0.05); dark blue, downregulated DE genes (FDR < 0.05). ::: {#fig1}

Figure 1—figure supplement 1.

Figure 1—figure supplement 2.

Figure 1—figure supplement 3.

Principal component analysis (PCA) of MIA and control RNA-seq samples showed developmental age accounts for the majority of variance across samples, with some additional separation between MIA and control groups (Figure 1b–c). We performed stage-specific differential expression analysis for male and female samples and for both sexes with sex as a covariate. For each time-point, DE genes were defined using a stringent false discovery rate (FDR) < 0.05 threshold and a more inclusive p < 0.05 threshold (Figure 1d,e, Supplementary file 1–!number(5)). Sex-stratified analysis indicated generally concordant DE differences between male and female samples following MIA, although DE effect sizes appeared stronger in females at E14.5 and E17.5, leading to a larger number of DE genes identified independently in females compared to males (Figure 1—figure supplement 2, Supplementary files 6–!number(15)). Overall, sex-stratified differences in DE genes were generally subtle, as demonstrated by high DE correlation and shared DE gene sets between sexes, with consistent findings across sexes among key DE genes (Figure 1—figure supplement 3). Samples from both male and female offspring were used for overall DE analysis models, with sex included as a covariate.

For the DE gene set passing the stringent FDR < 0.05 threshold, there were varying numbers of DE genes across time-points (Figure 1d). The strongest transcriptional signature was observed at E17.5 (2621 up and 3058 down), suggesting dramatic impact of MIA on cortical development at this time-point. P0 represented the subtlest transcriptomic signature, with no genes passing the stringent FDR < 0.05 threshold, and with the p < 0.05 DE genes showing a strong upregulation bias. We tested for overlap between DE RNA-seq genes and high confidence autism-associated genes in the Simons Foundation Autism Research Initiative (SFARI) gene database 8Basu et al.2009. The 82 mouse orthologs of SFARI ASD genes that were expressed at measurable levels were significantly enriched among E17.5 DE genes, with 25 upregulated and 19 downregulated DE genes at FDR < 0.05 (p=3.9e-04, hypergeometric test), and 28 upregulated and 28 downregulated DE genes at p < 0.05 (p = 2.8e-06, hypergeometric test) (Figure 1e, Supplementary files 16–!number(17)). This enrichment during peak transcriptomic dysregulation suggests involvement of ASD-relevant pathways in MIA etiology here and demonstrates general NDD relevance of our model. Detailed results of DE analysis are reported in Supplementary files 1–!number(15) and can be visualized using our interactive online browser.

We next sought to capture MIA-induced gene regulatory changes at the systems and network level across embryonic cortex development using WGCNA (39Langfelder and Horvath2008; 90Zhang and Horvath2005; Figure 2, Figure 2—figure supplement 1, Figure 2—figure supplement 2). In this approach, genes with correlated expression patterns are assigned into modules, enabling identification of gene sets with shared expression and function in neurodevelopment. WGCNA co-expression modules were arbitrarily named after colors, with Grey reserved for genes with no strongly correlated expression patterns (Figure 2a). Our initial analysis identified 10 modules (Figure 2—figure supplement 1a), of which correlated modules were combined to produce the final set of five modules and the unassigned Grey set. Modules were tested for association with sample age, MIA versus saline condition, and sex (Figure 2b). The modules captured dynamic trajectories of neurodevelopmental gene expression (Figure 2c). The Blue and Turquoise modules represented genes that increased or decreased in expression respectively, during neurodevelopment (Blue: p=9e-34; Turquoise: p=1e-29, Figure 2—figure supplement 1b). The other modules showed more complex patterns of expression. Each WGCNA module included DE genes, with differences in what time-point had the most extensive DE, indicating module-specific timing of perturbation associated with MIA (Figure 2—figure supplement 1c). DE genes were enriched for module-specific GO terms, with GO enrichment findings similar using all module genes and using an alternative rank-based gene set enrichment that is not dependent on FDR or p-value (56Mootha et al.2003; 75Subramanian et al.2005; Figure 2d, Supplementary file 18; Figure 2—figure supplement 2, Supplementary file 19).

Figure 2.

Figure 2—figure supplement 1.

Figure 2—figure supplement 2.

The WGCNA-resolved modules enabled mapping of DE signatures to neurodevelopmental processes that were acutely induced by MIA (Figure 3, Figure 3—figure supplement 1). There was significant association between genes in Grey (p = 0.007) and Green (p = 0.03) modules and MIA treatment, and marginal significance with MIA for the YeMaBl (p = 0.05) module (Figure 2—figure supplement 1b). The acute, initiating signaling pathways, found 6 hr after poly(I:C) injection captured by the Green module and by a subset of genes within the Grey module, show a strong signature of upregulated DE genes from E12.5 to E14.5 (Figure 3a). Analysis of the GO BP enrichment of upregulated DE genes in these modules identified activation of immune-related pathways, including defense response to virus, as well as angiogenesis and Vascular Endothelial Growth Factor A (VEGFA) signaling (Figures 2c and 3b, Figure 2—figure supplement 2a,c). We tested for protein-protein association networks among the acute E12.5 DE (FDR < 0.05) genes using STRING 77Szklarczyk et al.2019, and found significantly more interactions than expected by chance (observed edges = 161, expected edges = 46, enrichment = 3.5, STRING p < 1.0e-16) (Figure 3c). This intersecting protein network was associated with metabolism, hypoxia, and stress (Figure 3b). The network includes a highly interacting core gene set with hub nodes, such as Vegfa, that have many spokes and may be signaling factors that direct the MIA response (Figure 3b,c). Developmental expression profiles for DE genes that were associated with this interaction network, are shown in Figure 3d, and stratified by sex in Figure 3—figure supplement 1a. These signatures capture a complex but coherent transcriptional response involving hypoxia, immune, metabolic, and angiogenesis pathways that are strongly and transiently induced in fetal brain within 6 hr following mid-gestational MIA.

Figure 3.

Figure 3—figure supplement 1.

The most pervasive WGCNA-resolved DE signatures following MIA impacted a large proportion of genes in the Blue and Turquoise modules by E14.5 and peaked at E17.5 (Figure 1c). GO terms related to proliferation and cell cycle were enriched in the Blue module among downregulated DE genes, and terms related to neuronal differentiation and synapses were enriched in the Turquoise module among upregulated DE genes (Figure 2c). Considering the changes in proliferative and lamination markers associated with Blue and Turquoise genes, we sought to validate such changes at the protein level in independent samples (Figure 4, Figure 4—figure supplement 1). Four out of five DE genes tested validated with concordant changes at the protein level as determined by western blot (WB) analysis at E17.5, showing either statistical significance or trends in the same direction as observed in our transcriptomic analysis (Figure 4a–b, Figure 4—figure supplement 1). The DE signatures associated with Blue and Turquoise modules were evidenced by separation between MIA and saline samples in PCA plots at E14.5 and E17.5, with E17.5 MIA samples clustered closer to the P0 samples than saline age-matched counterparts (Figure 1b). To test for MIA-associated perturbation to temporal full transcriptome expression patterns, we generated a linear regression model using RNA-seq data principal components 1–5 and RNA-seq covariates to model age in saline samples. We then used this model to predict age in MIA samples (Figure 4c). E17.5 MIA samples exhibited older predicted age values relative to E17.5 saline controls (p = 4.5e-06, T-test), and E14.5 MIA samples showed a similar trend (p = 0.07, T-test). These findings suggest that mid-gestational MIA at E12.5 perturbs the timing of major downstream neurodevelopmental processes at E14.5 and E17.5 in embryonic cerebral cortex.

Figure 4.

Figure 4—figure supplement 1.

To characterize the neuroanatomical specificity of Blue and Turquoise dysregulation transcriptional signatures (Figure 4d), we performed a static immunofluorescent labeling comparison of late neurogenesis, cortical lamination and cell specification at E17.5 in an independent MIA cohort. E17.5 represents the end of the peak of neurogenesis during cortical development, before the neuronal-glial transition that occurs around E18.5. At this stage, it is expected that progenitor markers, such as PAX6, SOX9 and more specifically SOX2, will identify a small population of progenitors that are maintained and retain the ability to produce neurons, as well as subpopulations of proliferative cells commencing gliogenesis 50Martynoga et al.2012. PAX6+, SOX9+ and SOX2+ cell distribution patterns showed that progenitors localized correctly around the cortical ventricular zone in coronal brain sections; all analyzed cell populations were significantly reduced in MIA brains compared to controls (Figure 4e–f). The PAX6 finding in particular is consistent with our WB analysis at E17.5 (PAX6, p=0.039, two tailed Student’s t-test) (Figure 4b, Figure 4—figure supplement 1a), with protein levels that either resolved or were too subtle to be distinguishable by birth (p=0.766, two tailed Student’s t-test) (Figure 4a, Figure 4—figure supplement 1b–c,i). MIA brains showed a decrease in SOX9+ cells in the ventricular zone (VZ), which is expressed in progenitors 68Scott et al.2010. We also observed increased SOX9+ cells above the VZ in the neocortex, which can remain expressed in a subset of cells, including oligodendrocyte progenitors and maturing astrocytes 37Klum et al.201876Sun et al.2017 (SOX9 VZ, p=0.0198; nCtx, p=0.0058, two tailed Student’s t-test). Reduced expression of SOX2, a progenitor marker that is restricted to the VZ, was also observed in MIA brains (SOX2 VZ, p=0.0006 two tailed Student’s t-test) (Figure 4f). Overall, these data are consistent with the transcriptomic downregulation signature of proliferating cells in MIA brains.

To determine if the decrease in progenitor markers might coincide with decreased proliferation, we next labeled cells undergoing proliferation, using the marker KI67 (any stage of cell proliferation) and the Phospho Histone 3 (PH3) marker, which labels a subset of cycling cells in mitosis (M) phase 30Hendzel et al.199784Walton et al.2019. We also examined the post-mitotic intermediate progenitor marker TBR2. Analysis of KI67+ and PH3+ cells showed reduced numbers of labeled cells along with a qualitative reduction in numbers and ectopic positioning of postmitotic intermediate progenitor cells in the MIA group (TBR2+). Reduction in Ki67+ cells was observed in several regions of the developing cortex, including zones defined as S1 and S1DZ 14Choi et al.201670Shin Yim et al.2017, as well as the ganglionic eminences (Figure 4e). Quantification of Ki67+ cells and PH3+ cells in the neocortex confirmed the MIA-induced reduction in progenitor and mitotic populations at E17.5 (Ki67, p=0.0198; PH3, p=0.0003, two tailed Student’s t-test). The ratio of PH3+ (mitosis) cells to Ki67+ (all stages of cell division) cells was then used to determine whether cell cycle kinetics or the number of actively dividing cells in M phase was different in the MIA-exposed group. Ki67/PH3 ratios indicated a reduced proportion of progenitor cells in M phase in MIA brains at this specific time-point (p=0.0005 Chi-squared) (Figure 4f). These results provide evidence that the reduced proliferation signatures following MIA that were observed in the transcriptomic profiles are associated with changes at the protein level. Moreover, within the reduced progenitor cell populations at E17.5, the number of cells in the mitosis stage of the cell cycle are decreased to a greater extent.

MIA has been associated with altered cortical thickness and reported to cause cortical dysplasia in young adult MIA offspring 72Smith et al.200773Soumiya et al.201182Tsukada et al.2015, with stereotypic dysplasia impacting lamination localized to the somatosensory cortex reported in one MIA model 14Choi et al.201670Shin Yim et al.2017. A number of cortical layer markers were DE at E17.5 in RNA, with validated changes at the protein level (Figure 4a–b, Figure 3—figure supplement 1b). To examine cortical lamination, distribution of cortical layer-specific markers was analyzed via IHC in the same samples as above (Figure 5, Figure 4—figure supplement figure 1 h-i).

Figure 5.

Examination of the distribution of the lamination markers TBR1 (layer 6), CTIP2 (layer five and a subset of cells in layer 6), SATB2 (layers 2–4, and a subset of cells in layer 5), and CUX1 (layers 2 and 3), revealed MIA-induced cortical layer-specific alterations in lamination, thickness, and cell density at E17.5 (Figure 5a–b,e–f). Overall, MIA brains appeared to be at a more mature stage of lamination at E17.5, showing a full separation of cortical deep layers more typical of later ages and in contrast to the partial layer overlap found in control offspring at E17.5 (Figure 5a–b). Similarly, SATB2 and TBR1 co-staining suggested differences in lamination and migration state between the two groups. Oval cell body shapes were observed, which resembled radially migrating cells 15Cooper2013 were observed in both SATB2 and TBR1 cells in control tissue, whereas these cells showed a more rounded/post-migratory cell body shape in the MIA brains (Figure 5c–d). The number of TBR1+ and CTIP2+ cells, but not SATB2+ cells, were also significantly reduced in the neocortex from MIA offspring (TBR1+ cells, p=0.0138; CTIP2+ cells, p=0.0119; SATB2+ cells, p=0.306, two tailed Student t-test) (Figure 5e). The thickness of these layers and also a more external layer represented by CUX1, were measured by comparing the radial spread of stained cells within the neocortex. We found a significant reduction in CTIP2 layer thickness and a similar trend in TBR1 layer thickness in the MIA brains. In contrast, SATB2 showed an MIA-induced significant increase in layer thickness with no change in the more external layer represented by CUX1 (TBR1 thickness, p=0.081; CTIP2 thickness, 0.024; SATB2 thickness, p=0.039; CUX1, p=0.430, two tailed Student’s t-test). (Figure 5f). These findings are consistent with individual trajectories of DE genes and WB analyses (Figure 4a–b, Figure 4—figure supplement 1h). Finally, to assess cortical gross anatomy, we looked for MIA-associated dysplasic alterations across the brain 14Choi et al.2016 and measured the thickness of the developing somatosensory cortex at two anteroposterior regions at E17.5. In contrast to previous reports 14Choi et al.2016, we did not find obvious cortical dysplasia or gross alterations in brain anatomy in the analyzed rostral and caudal regions, including somatosensory cortex and S1DZ (Figure 5e–f).

To investigate other cell populations that may be impacted by MIA during development, we analyzed astrocytes, oligodendrocytes and GABAergic interneuron cell populations (Figure 6, Figure 4—figure supplement 1k). Concordant with RNA-seq results, the number of cells expressing GFAP, a protein found in radial glia and a subset of astrocytes 23Eng et al.2000, was increased in MIA offspring, mapping specifically to cells that seemed to emerge ventral from the neocortical VZ (Figure 6a–b,q, Figure 4—figure supplement 1k). Since astrocytes also express SOX9 37Klum et al.201876Sun et al.2017 and radial glia also express PAX6, which is decreased by MIA, this finding is consistent with the possibility that the increased SOX9+ cells in the dorsal regions of the neocortex, described earlier, may be indicative of early differentiated astrocytes (Figure 4f). OLIG2, which is expressed in oligodendrocyte progenitors and mature cells as well as transiently in astrocyte progenitors 18Dimou et al.2008, was also upregulated at E17.5 in the MIA group. OLIG2+ cells were increased in MIA brains compared to controls (n = 3 per group, p=0.0134, Student’s t-test) (Figure 6c–d,k, Figure 4—figure supplement 1k). PDGFRα, a marker of early developing oligodendrocytes, was also increased in MIA brains compared to controls (n = 7 per group, zone 1 p<0.0001, zone 2 p=0.0063, zone 3 p=0.0125, Student’s t-test) (Figure 6e–f,l). Thus, these experiments provide evidence that MIA leads to increased numbers of additional cell types (putative astrocytes and oligodendrocytes) that are typically born at later ages during corticogenesis, reinforcing findings of accelerated cell development that is supported by an increased representation of later-born cell types at E17.5 in brains exposed to MIA.

Figure 6.

Finally, we examined DLX2 and GAD67, markers of committed and maturing cortical interneurons. Although Dlx2 expression was downregulated in the MIA group from bulk tissue dissections at E17.5, the number of cells positive for DLX2 protein in the neocortex were comparable when assessed via IHC in the MIA and control groups (n = 3 per group, total DLX2 p=0.5494; by zones DLX2, zone 1 P=0.9848, zone 2 P=0.6903, zone 3 P=0.2714, Student’s t-test) (Figure 6g–h and m,o). In contrast, GAD protein in this DLX2+ subpopulation, assessed by the pan-GABAergic marker GAD67+, was increased across the neocortex of MIA brains, and was not restricted to any specific zones of the neocortex (n = 3 per group, total DLX2+/GAD67 p=0.0004; by zones DLX2+/GAD67, zone 1 P=0.2277, zone 2 P=0.1352, zone 3 P=0.2862, Student’s t-test) (Figure 6i–j and n,p). Since GAD67 expression increases in developing cortical interneurons, the decrease in Dlx2 transcript expression but no difference in cells expressing DLX2 protein could indicate increased interneuron maturation but no change in interneuron numbers at this age. While neuroanatomical characterization here was limited to E17.5 and interrogated a subset of genes and processes of interest, our IHC findings correlate with transcriptional profiles, providing validation of MIA-associated DE in developing cortex.

Discussion

Studies leveraging MIA models in mice have demonstrated the role of maternal immune signaling as a risk factor for cellular, anatomical, and behavioral outcomes in offspring. Nonetheless, a systems-level perspective of the early neurodevelopmental trajectory of pathology associated with MIA has been lacking. Our results offer a time course resolved map of molecular targets aberrantly transcribed in developing cerebral cortex following MIA. MIA-associated transcriptional signatures were robust, changed greatly across the surveyed time- points, and revealed novel pathways and processes involved in pathology. Among MIA-associated DE up-regulated and down-regulated genes, there was strong enrichment for ASD-relevant loci, indicating the general translational relevance of our MIA model to ASD. Both male and female samples were examined, and we did not find evidence for strong sex-specific DE signatures or impacts on neuroanatomy. Elucidating the trajectory of transcriptional changes across embryonic cortical development following mid-gestational MIA represents a valuable step toward understanding the progression of MIA-induced pathology in the developing brain of offspring.

Our results suggest two waves of transcriptional pathology following MIA. First, we found initial induction of immune, metabolic, and angiogenesis gene expression modules specific to MIA exposure and largely restricted to the days following MIA. There is overlap in the acute induction signatures identified here and acute changes reported in other mid-gestational MIA models, including those using a range of injection gestational ages and rat models of LPS 13Carpentier et al.201354Meyer201452Meyer et al.200659O'Loughlin et al.201764Richetto et al.201472Smith et al.200787Wu et al.2017. Previous studies linked maternal immune challenge, prenatal brain hypoxia, and impacts on proliferation 74Stolp et al.2011, neurogenesis, and cortical lamination 13Carpentier et al.2013. Thus, it is possible that the acute signatures identified here and mapped to the Green and YeMaBl modules as well as a subset of genes in the unassigned Grey module, capture a generalized response in fetal brain across mid-gestational MIA models. Notably, induction of hypoxia, cellular metabolism and neuroimmune transcriptional signatures, particularly captured here by the Green module, persist at reduced levels at least through P0, and may indicate lasting perturbation to these pathways. Significant association of the Grey unassigned module with the MIA phenotype was driven by a subset of genes that showed strong acute induction at E12.5 and E14.5. By altering WGCNA parameters, we could assign the Grey module genes to co-expression modules but those alterations increased splitting, resulting in a high number of modules and reduced module coherence without clarifying additional MIA-specific modules made up by the relevant Grey module MIA-associated genes. Thus, we opted for reduced splitting to capture major co-expression trends.

Identifying the early biological processes triggered by MIA are of critical interest for understanding how MIA initiates pathology. Our analysis identified potential drivers of these induced signaling pathways, including targets of substantial interest such as Vegfa, that link hypoxia, neurogenesis and cortical angiogenesis 4Argaw et al.201278Takahashi and Shibuya2005. Vegfa was upregulated in MIA-exposed animal at all time-points and was among the strongest DE genes at E12.5 and P0. Vegfa is critical for cortical development promoting neurogenesis and neurite outgrowth, as well as angiogenesis 9Bayless et al.201648Mackenzie and Ruhrberg201266Rosenstein et al.2010. Upregulation of Vegfa is found in several NDDs and psychiatric disorders, including SZ and depression, and was also reported in a different MIA mouse model 5Arrode-Brusés and Brusés201231Hohman et al.201533Iga et al.200734Isung et al.201236Kahl et al.200941Lee and Kim201245Lizano et al.201888Xie et al.2017. Induction of angiogenesis in the developing cortex also independently causes a transition from proliferative state to neuronal maturation 9Bayless et al.201635Javaherian and Kriegstein200979Tan et al.2016. Other genes of interest include Pdk1 (pyruvate dehydrogenase kinase 1), which regulates glucose metabolism 58Newington et al.2012, Ldha (lactate dehydrogenase A), which facilitates glycolysis and participates in brain angiogenesis 43Lin et al.2018, Bnip3 (BCL2 interacting protein 3), which is involved in cell death and mitochondrial autophagy 44Liu and Frazier2015, and Il20rb (interleukin 20 receptor subunit beta), a receptor subunit for pro-inflammatory interleukins IL-19, IL-20 and IL-24 1Akdis et al.2016. The early perturbed pathways observable within hours of MIA induction may be critical initiating steps leading to the altered general timing of fetal brain development described here. Our findings expand understanding of the molecular basis of these MIA-induced initial perturbations and reveal potentially key drivers of these changes.

The second component of the neurodevelopmental MIA signature defined here was perturbation to normal patterns of both proliferative and maturational neurodevelopmental modules that manifest several days after MIA induction. Initially, at six hours after MIA, there is a shift in expression of differentiation and neuronal maturation processes (e.g. synaptogenesis) to increased proliferative processes (e.g. positive regulation of cell cycle). This pattern parallels results described for the acute MIA response from other mouse and rat MIA models, where similar pathways were down- or up-regulated in fetal brain in the hours following MIA 55Money et al.201860Oskvig et al.2012. At this developmental stage, our data suggest that MIA induces an acute expenditure of the neocortical, and potentially other, progenitor pools, which in turn leads to an expedited maturation of cells that would normally be generated later during development. Specifically, by E17.5, multiple progenitor and cell cycle markers were significantly decreased in the MIA brains while concordantly, we observed an increase in maturational markers and laminar position that resembled a more developed cortex. Finally, astrocytes, oligodendrocytes, GABAergic interneurons and radially migrating excitatory neuron precursors were effected, suggesting an unbiased impact on progenitors of multiple brain cell types. Although these findings are indicative of overall decreased neurogenesis by E17.5, specific birth-dating and time-course experiments are needed to test whether progenitor cells exit the cell cycle differently and if overall neurogenesis is perturbed beyond the E17.5 time-point reported here. Nonetheless, our time course RNA-seq design revealed that the transcriptional changes following MIA can be separated into discrete biological and temporal expression modules that support altered proliferation and maturation of cell types following MIA.

Our findings of disrupted cortical lamination following the acute MIA signature are consistent with, but not identical to previous studies showing cortical dysplasia associated with MIA-induced aberrant behavior in offspring 14Choi et al.201670Shin Yim et al.2017. In those studies, dysplasia resolved at later postnatal stages but defects in lamination persisted into adulthood. Although we did not find evidence of aberrant cortical patches in the somatosensory cortex, secondary motor cortex or S1 and S1DZ, as previously reported, our results do show altered cortical lamination at E17.5. Our findings suggest that perturbation to neuronal migration, potentially due to differences in an early exit of pyramidal neuron precursors from the cell cycle that led to cells further along in their radial migration/lamination, developmental progression, at least in part contribute to MIA-associated lamination differences identified at E17.5 here. Thus, results from multiple studies using MIA models to date are consistent in showing a detrimental impact of MIA on cortical patterning; however, we also acknowledge that our and other studies may have differential variables, including onset of MIA, dose, mouse background, and endpoint measurements, which will be important factors to assess in future studies. Finally, bulk tissue RNA-seq is limited in cell type-specific resolution and signatures could additionally be in part driven by cortex-adjacent structures, such as primordial hippocampus, that may have also been captured at the youngest age. Future work using single-cell approaches is needed to resolve the impacted cell types and compare transcriptional pathology following MIA.

While significantly subtler in overall effects compared to earlier time-points, we also see evidence for altered neuronal transcription and lasting immune activation at P0. By P0, these signatures were reduced, which suggests a transient or attenuating effect on cortical patterning. It is also possible that biological variability in neonatal mice or technical variability impacted detection of gene expression changes in these P0 samples. Nonetheless, the subtle changes in gene expression at P0 are consistent with deficits in synapse formation that require changes in expression of immune molecules in neurons from newborn MIA offspring 22Elmer et al.2013.

Consistent with links between NDDs and developmental timing, atypical neurodevelopmental gene expression networks have been reported in a cortical neurodevelopmental model using autism-patient-derived induced pluripotent stem cells 9Bayless et al.201667Schafer et al.2019. The majority of the ASD-relevant genes that were DE following MIA (Figure 1e) were associated with perturbation to the modules associated with neurodevelopmental processes. These ASD-relevant genes map to functional and biological pathways generally associated with ASD and NDDs, including gene regulation, neurogenesis, and synapse development, structure, and function. The broad pathway overlap and presence of up- and down-regulation among DE ASD-relevant genes supports their generalized relevance. These overlapping genes represent targets of interest for future studies investigating parallel mechanisms between genetic and MIA-associated NDD risk. Our findings suggest some aspects of the transition from neurogenesis to astrogenesis appear to occur at an earlier stage following MIA based on similar signatures of SOX9+ and GFAP+ cells reported in literature 7Bansod et al.2017, and increased OLIG2+ cells and more mature GAD67+ DLX2+ interneurons in the cortex following MIA. These findings are additionally consistent with recent reports using poly(I:C) MIA models that identified perturbed GABAergic interneuron development 81Thion et al.2019, and maturation of GAD+ neuroblast after birth 83Vasistha et al.2020. Future work is needed to determine if and how MIA-associated dynamic transcriptional changes manifest postnatally and whether perturbation to embryonic development impacts cortical cytoarchitecture and function at later postnatal stages.

Altogether, our analyses indicate that MIA drives acute and lasting transcriptional changes that impact specific cell populations during forebrain development, disrupting normal developmental timing and causing increased numbers of astrocytes, oligodendrocytes, and altered molecular properties of GABAergic interneurons. The timing of MIA during gestation is known to cause a range of phenotypes in offspring, and thus the findings here regarding mid-gestational MIA at E12.5 may not generalize to other models with different timing of MIA induction, especially for later gestational MIA after neurogenesis is largely complete. Further, the strong embryonic signatures we observe were generally consistent across sexes, raising the question of whether secondary insults or sex-specific protective effects contribute to some reports of sex differences in pathology following MIA 29Haida et al.201957Naviaux et al.201372Smith et al.2007. Nonetheless, the overlapping results between our study and published single time-point studies of the acute MIA fetal brain transcriptional response suggest some consistency across bacterial and viral mimics, and rat and mouse models. Our study substantially expands previous work, revealing new molecular pathways and providing a temporal map of MIA-induced changes across fetal neurodevelopment that will be relevant to understanding the pathology of, and informing future treatments for, NDDs.

Materials and methods

Animal care and use

All studies were conducted in compliance with NIH guidelines and approved protocols from the University of California Davis Animal Care and Use Committee. C57BL/6N females (Charles River, Kingston, NY), were bred in house and maintained on a 12:12 hr light dark cycle at 20 ± 1°C, food and water available ad libitum. Mice utilized in this study do harbor Segmented Filamentous Bacteria (SFB) 24Estes et al.2019. Males and females embryos were analyzed across experiments, sex was determined as previously described 51McFarlane et al.2013.

Maternal immune activation

Assessment of females baseline immunoreactivity before pregnancy and maternal immune activation were performed as previously described 25Estes et al.2020. Since female mice display a wide range of baseline immunoreactivity to poly(I:C) that dictates susceptibility and resilience of offspring from later pregnancies to MIA, virgin females in the lowest 25th percentile of immune responsiveness to poly(I:C) (IL-6 serum levels) were excluded from the study to reduce variability and ensure sufficient immune responsiveness of all included dams. E12.5 pregnancies were determined by visualizing vaginal seminal plug (noted as E0.5) and by body weight increase. IL-6 serum levels were measured at 2.5 hr post-injection. Temperature, weight, and sickness behavior were tracked as previously described 25Estes et al.2020.

Tissue dissections for RNA-seq and immunoblotting

Dissections of dorsal telencephalon/pallium or developing cortex were performed at E12.5 + 6 hr, E14.5, E17.5, and P0 from control (saline) and poly(I:C)-exposed MIA litters. At E12.5, due to the smaller embryo size, dissections were made via a diagonal cut that maximized collection of dorsal telencephalon, but could have also included some subpallium/ventral telencephalon and superficial tissues. For later time-points, dissection was of pallium/developing cortex. Potential variance in the individual dissections should be mitigated by the inclusion of multiple replicates across both groups. Dissections included both hemispheres from males and females, with exact numbers and sample details reported in Supplementary file 20. Independent litters were used for RNA-seq and immunoblotting.

RNA-seq

Following tissue collection as described above, total RNA was isolated using Ambion RNAqueous total RNA Isolation Kit and assayed using Agilent RNA 6000 Nano Bioanalyzer kit/instrument. Stranded mRNA libraries were prepared using TruSeq Stranded mRNA kits. Eight to 12 samples per lane were pooled and sequenced on an Illumina HiSeq 4000 instrument using a single-end 50 bp protocol. Reads were aligned to mouse genome (mm9) using STAR (version 2.4.2a) 19Dobin et al.2013 and gene counts produced using featureCounts 42Liao et al.2014 and mm9 knownGenes. Quality assayed using FastQC 3Andrews2010 and RSeQC 85Wang et al.2012, and samples that exhibited 3′ bias or poor exon distribution were discarded. Raw RNA-seq fastq files and a gene count matrix is available on GEO GSE166376.

Immunoblotting

Tissue dissections from male and female littermates from at least three litters per condition were performed as described above and collected in HBSS, immediately frozen on dry ice, and stored at −80°C until processed. Samples were lysed in 50 mM Tris HCl, pH 8, 140 nM NaCl, 1 mM EDTA, 10% glycerol, 0.5% NP40 and 0.25% Triton with protease inhibitor cocktail (Roche). After sonication, samples were spun down and the supernatant was used for a BCA Bradford assay using the Spectramax 190 plate reader to assess protein concentration using a standard curve. Twelve or 18 μg of protein were run on a 10 or 12% Tris acetate gel using the Mini-PROTEAN system (BioRad). Membranes were blocked in Odyssey blocking buffer (TBS; LiCor) and probed with the indicated primary antibodies overnight at 4°C. Resolved proteins were visualized in two channels using fluorescent secondary antibodies at 680 and 800 nm on an Odyssey Clx infrared imaging system (LiCor). Specific band intensities for all detected antibodies were quantified using the manufacturer’s software ImageSoft (LiCor) and normalized to GAPDH loading control. Antibodies used were anti-PAX6 (1:250; cat #PRB-278P-100; Covance, Princeton NJ.), anti-TBR1 (1:500; cat#ab31940; Abcam, Cambridge, MA), anti-CTIP2 (1:250; cat# ab18465; Abcam, Cambridge, MA), anti-CUX1 (1:200; cat#A2213; Abclonal, Woburn, MA) and anti-SATB2 (1:500; cat#ab51502; Abcam, Cambridge, MA).

Gross anatomy, immunohistochemistry, and image quantification analyses

Histology was performed in triplicate sections in brains using male and female pups from at least two litters. Sections were selected to focus on the developing somatosensory cortex. Details on sample individual sample usage for each IHC experiment is detailed in Supplementary files 21 and !number(22). Embryos used in IHC experiments came from independent litters versus RNA-seq and immunoblotting. Morphological parameters were measured using FIJI ImageJ (NIH), by an experimenter blinded to group and sex. Brain regions were identified based on anatomical landmarks as previously described 28Gompers et al.2017. No data points were excluded. For IHC staining, fetal brains were dissected and fixed with 4% paraformaldehyde/PBS solution overnight at 4°C. Tissues were then equilibrated in approximately 15 mL of 30% sucrose/PBS solution until they sank to the bottom of a conical tube. Equilibrated brains were embedded in Optimum Cutting Temperature (OCT) compound (Tissue-Tek, Torrance, CA) and frozen on dry-ice. OCT-embedded brain blocks were cryo-sectioned on a coronal plane (30 µm). Immunostaining was performed in free-floating sections with agitation. Sections were washed five times in PBST (PBS with 0.05% Triton X-100, 5 min each) and antigen retrieval was performed using 1x Citrate buffer pH6.0 antigen retriever solution (cat# C9999; Millipore-Sigma, Burlington, MA), at 60°C for 1 hr. Sections were washed five times in PBST (5 min each), permeabilized in PBS containing 0.5% Triton for 20 min and blocked for 1 hr at room temperature in 5% milk/PBST. Primary antibodies were incubated overnight at 4°C with orbital agitation (40–50 rpm). All antibodies used for this study were validated and their specificity validated. The following primary antibodies were used: anti-PAX6 (1:250; cat #PRB-278P-100; Covance, Princeton NJ.), anti-KI67 (1:500; cat#12202; Cell Signaling, Danvers, MA), anti-SOX9 (1:500 cat# AF3075, R and D systems, Minneapolis, MN), anti-PH3 (1:500 cat# 9701, Cell Signaling, Danvers, MA), anti-TBR2 (1:500; cat#14-4875-82; Thermo Fisher Scientific, Waltham, MA.), anti-TBR1 (1:500; cat#ab31940; Abcam, Cambridge, MA), anti-CTIP2 (1:250; cat# ab18465; Abcam, Cambridge, MA), anti-CUX1 (1:200; cat#A2213; Abclonal, Woburn, MA.), anti-SATB2 (1:500; cat#ab51502; Abcam, Cambridge, MA.), anti-DLX2 (1:200; generous gift from John Rubenstein, UCSF), anti-GFAP (1:250, cat#Z0334; Agilent Dako, Santa Clara, CA), anti-OLIG2 (1:500, cat#AB9610, MilliporeSigma, Burlington, MA). After primary antibody incubation, free-floating sections were washed five times in PBST (5 min each). Species-specific fluorophores-conjugated IgG (1:1000; Invitrogen-Thermo Fisher Scientific, Waltham, MA) were used as secondary antibodies (45 min, RT). 40,6-Diamidino-2-phenylindole (DAPI) (1:10000; Millipore-Sigma, Burlington, MA) was used for nuclear staining (20 min, RT). For cell counting, boxes were drawn on specific areas of the neocortex for progenitor and mitosis markers; SOX2, SOX9, and PH3. This box was 300 × 300 pixels on an image of the neocortex obtained using a ×20 objective. To be consistent, the ventral aspect of the box was drawn along the VZ and one corner touched the medial aspect of the VZ. The same strategy was used to draw boxes from Ki67 images except the box was increased in size to 500 × 500 pixels to account for the additional Ki67+ cells seen in neocortical areas above the SVZ. SOX9 boxes were drawn along the dorsal border edge (600 pixels in width) of the neocortex and extended ventrally two thirds the length of the neocortex in order to avoid the VZ/SVZ region. For laminar zone counts, a box spanning 300 pixels was drawn along the VZ and the dorsal/ventral proportions were equally divided into three zones. These box sizes allowed for a sufficient number of cells to be counted from each tissue. Laminar thickness analyses were performed by measuring the spread of all positive cells within the cortical plate per individual marker. For all quantification analyses, one to three sections per brain utilizing the developing hippocampus as a landmark for positioning were used.

Bioinformatics analysis

Bioinformatic analysis was performed using R programming language version 3.5.1 63R Development Core Team2015 run in RStudio integrated development environment version 1.2.1269 80Team R2018. Plots were generated using ggplot2 R package version 3.1.0 86Wickham2009. Heatmaps were generated using pheatmap R package 1.0.10 38Kolde2018. The analysis scripts are available at: https://github.com/NordNeurogenomicsLab/Publications/tree/master/Canales_eLife_2021; 11Canales2021; copy archived at swh:1:rev:28836c8758908e130f1f6d8bfb3b6a112c6cbd1b.

Differential expression analysis

Raw count data for all samples were used as input along with sample information for differential expression analysis using edgeR 65Robinson et al.2010. Genes with minimum log₂reads per kilobase per million (RPKM) expression of −2 in at least two samples were included for analysis, resulting in a final set of 17,195 genes for differential testing. For time-point differential expression analysis genes with count per million (CPM) >0.1 in at least two samples were considered, resulting in E12.5, 16,396; E14.5, 16,658; E17.5, 16,229; P0, 16,613 genes in differential expression analyses. Principal component analysis indicated that the strongest driver of variance across samples was developmental age. Tagwise dispersion estimates were generated, and differential expression analysis was performed with edgeR using a generalized linear model, separately for each developmental time-point, including sex, sequencing run factor, and treatment as the variable for testing. Normalized expression levels were generated using the edgeR rpkm function. Normalized log₂ RPKM values were used for plotting of summary heatmaps and of expression data for individual genes.

SFARI gene set enrichment analysis

The autism risk gene-set was downloaded from https://gene.sfari.org/ on 02-26-2019 and is included in this manuscript as Supplementary file 17. High confidence risk genes, annotated as ‘gene-score’ 1 and 2 were selected, and their orthologs were found using getLDS function from the biomaRt R package 21Durinck et al.200920Durinck et al.2005. Overlap of up- and downregulated DE genes with SFARI ortholog genes was calculated and plotted using a custom R script. Statistical significance of overlap was tested with hypergeometric test using the following R script: sum(dhyper(t:b, a, n - a, b)):

• t = number of overlapping DE genes with SFARI gene orthologs with gene_score 1 or 2
• b = number of SFARI gene orthologs with gene_score 1 or 2
• a = number of DE genes at a developmental stage
• n = number of genes at a developmental stage sum returns the sum of all the values present in its arguments. dhyper calculates density for the hypergeometric distribution.

WGCNA

We used the WGCNA R package, version 1.66 40Langfelder and Horvath201239Langfelder and Horvath2008 to construct signed co-expression networks using the entire dataset containing 24,015 genes. After the network construction, the gene set was filtered for minimal gene expression at an RPKM value of 0.25 or higher in at least two sample, resulting in a dataset consisting 17,195 genes. A correlation matrix using the biweight midcorrelation between all genes was computed for all relevant samples. The soft thresholding power was estimated and used to derive an adjacency matrix exhibiting approximate scale-free topology (R² >0.85). The adjacency matrix was transformed to a topological overlap matrix (TOM). The matrix 1-TOM was used as the input to calculate co-expression modules using hierarchical clustering. Modules were branches of the hierarchical cluster tree base, with minimum module size set to 10 genes. Pearson’s correlation coefficients were used to calculate correlation between sample traits (e.g., sex, treatment) and modules. The expression profile of a given module was summarized by the module eigengene (ME). Modules with highly correlated MEs (correlation >0.80) were merged together. The module connectivity (kME) of each gene was calculated by correlating the gene expression profile with module eigengenes. The module connectivity (kME) of each gene was calculated by correlating the gene expression profile with module eigengenes. Genes with no network correlation were placed into the module Grey. Following manual data inspection further highly correlated modules were merged: BrRePi: Brown, Red, and Pink, YeMaBl: Yellow, Magenta, and Black.

Gene Ontology enrichment analysis

Mouse Gene Ontology (GO) data was downloaded from Bioconductor (org.Mm.eg.db). We used the TopGO R package version 2.34.0 2Alexa and Rahnenfuhrer2019 to test for enrichment of GO terms. For the analysis presented here, we restricted our testing to GO Biological Process annotations and required a minimal node size (number of genes annotated to GO terms) of 20. We used the internal ‘weight01’ testing framework and the Fisher test, a strategy recommended for gene set analysis that generally accounts for multiple testing comparisons. For GO BP analysis, we reported terms with p-value<0.05. For all enrichment analysis, the test set of DE genes was compared against the background set of genes expressed in our study based on minimum read-count cutoffs described above. Heatmaps showing positive log₂ (expected/observed) values were plotted for GO terms of interest.

Protein–protein interaction

Protein–protein interaction enrichment and network generation for E12.5 DE (FDR < 0.05) gene sets was performed using STRING 77Szklarczyk et al.2019, version 11, considering only experimentally and text mining interactions, with at least medium interaction confidence score. Disconnected nodes were removed from the network.

Data analysis

For all experiments, data was collected from at least two (usually three, and sometimes more, as indicated in each analysis) experiments and is presented as the mean ± SEM. Protein validation data were analyzed by Student’s t-test or one-way ANOVA, followed where appropriate by Tukey’s honestly significant difference post hoc test (Graphpad Prism v.7). Significance was defined as: *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001. Samples for RNA-seq were randomly collected across litters and processed blind to experimental condition. For differential gene expression analysis using edgeR, differences were considered statistically significant at high stringency for FDR < 0.05, or with reduced stringency for p values < 0.05.