Zbtb16 regulates social cognitive behaviors and neocortical development.

Zinc finger and BTB domain containing 16 (ZBTB16) play the roles in the neural progenitor cell proliferation and neuronal differentiation during development, however, how the function of ZBTB16 is involved in brain function and behaviors unknown. Here we show the deletion of Zbtb16 in mice leads to social impairment, repetitive behaviors, risk-taking behaviors, and cognitive impairment. To elucidate the mechanism underlying the behavioral phenotypes, we conducted histological analyses and observed impairments in thinning of neocortical layer 6 (L6) and a reduction of TBR1+ neurons in Zbtb16 KO mice. Furthermore, we found increased dendritic spines and microglia, as well as developmental defects in oligodendrocytes and neocortical myelination in the prefrontal cortex (PFC) of Zbtb16 KO mice. Using genomics approaches, we identified the Zbtb16 transcriptome that includes genes involved in neocortical maturation such as neurogenesis and myelination, and both autism spectrum disorder (ASD) and schizophrenia (SCZ) pathobiology. Co-expression networks further identified Zbtb16-correlated modules that are unique to ASD or SCZ, respectively. Our study provides insight into the novel roles of ZBTB16 in behaviors and neocortical development related to the disorders.

Introduction

ZBTB16 (PLZF) encodes a transcription factor, which contains a BTB/POZ protein-protein interaction domain in its N-terminal and a C2H2-type zinc finger DNA binding domain in its C-terminal, playing key roles in many biological processes such as stem cell maintenance and proliferation, cell differentiation, spermatogenesis, musculoskeletal development, hematopoiesis, apoptosis, chromatin remodeling, metabolism, and immunity[1,2]. Primary studies of spontaneous luxoid (lu) mutant of Zbtb16 in mice have shown that Zbtb16 is essential for skeletal development and germ cell self-renewal[3-5]. A single nucleotide variant (c.1849A>G; p.Met617Val) in the C2H2-type zinc finger domain of ZBTB16 has been identified as a causative mutation for skeletal defects, genital hypoplasia, and mental retardation (SGYMR)[6,7]. The SGYMR individual with this homozygous mutation in ZBTB16 showed intellectual disability, microcephaly, craniofacial dysmorphism, short stature, skeletal anomalies such as thumb deficits, and hypoplasia of the ulnae, retarded bone age, and hypoplastic genitalia[6,7].

A missense heterozygous mutation (c.1319G>A; p.Arg440Gln) of ZBTB16 was recently identified in brothers with ASD[8]. Moreover, other studies reported that a nonsense heterozygous mutation in ZBTB16 (c.1741A>T; p.Lys581*) in SCZ patients[9,10]. There are additional reports on the association between ZBTB16 and SCZ[11-13]. ASD is a heterogeneous neurodevelopmental disorder that causes pervasive abnormalities in social communication, as well as repetitive behaviors and restricted interests. The etiology of ASD is thought to involve complex, multigenic interactions, and possible environmental contributions[14]. SCZ is also a heterogeneous neuropsychiatric disorder characterized by positive symptoms (hallucinations and delusions), negative symptoms (flat affect, avolition, and social impairment), and cognitive impairment[15]. The biological mechanisms underlying ASD and SCZ are not fully understood. However, it is well-known that there are common overlapping mechanisms such as genetics, ethology, and brain dysfunction underlying the pathology of ASD and SCZ[16,17].

The involvement of ZBTB16 in brain development has been reported in several studies. Zbtb16 expression begins at embryonic day 7.5 in the neuroepithelium of the mouse embryonic brain and is eventually expressed in the entire neurectoderm at later stages[18]. ZBTB16 is expressed in human embryonic stem cell (ES)-derived neural stem cells (NSC) and primary neural plate tissue, playing a role in maintenance, proliferation[19], and neuronal differentiation[20]. A recent study has reported reduced cortical surface area and a number of deep-layer neurons in the neocortex of Zbtb16 mutant mice at neonatal stages[21]. In addition, Zbtb16 mutant mice showed an impairment of recognition memory in the novel object recognition test[21]. These studies indicate the involvement of Zbtb16 in neocortical development. However, these previous findings were primarily limited to the embryonic and neonatal periods, whereas the role of Zbtb16 in the adult brain, particularly the association with brain functions and behaviors is unknown.

In this study, we aimed to uncover the roles of ZBTB16 in brain functions and behaviors, as well as how ZBTB16 might be biologically relevant to disorders such as ASD and SCZ by focusing on behaviors and brain development. To address this question, we utilized a Zbtb16 homozygous mutant (Zbtb16 KO) mouse and characterized the behavioral and neocortical phenotypes. Zbtb16 KO mice displayed ASD-like and SCZ-like behaviors, and impairments in neocortical thickness, and a reduction of L6 neuronal numbers in the PFC. We also found increases in the numbers of dendritic spines and microglia, and oligodendrocyte developmental abnormalities resulting in impaired neocortical myelination. Finally, we characterized the Zbtb16 transcriptome in the PFC by conducting RNA-sequencing (RNA-seq) and identified that Zbtb16 regulates genes known to be involved in ASD, SCZ, and neocortical maturation including myelination. Co-expression gene networks identified the disorder-specific modules for ASD and SCZ, respectively. These results demonstrate that Zbtb16 plays an essential role in both ASD-like and SCZ-like behaviors via neocortical development, particularly deep layer formation, spinogenesis, and myelination. Taken together, our study demonstrates that Zbtb16 is involved in shared neurodevelopmental features that are at risk in both ASD and SCZ.

Materials and methods

Mice

B6.C3-Zbtb16/J mice were purchased from The Jackson Laboratory (#000100)[3,4]. Genotyping was performed using the following primers; for Zbtb16: 23559, F-5’-CCACCTCTTTCGGTCTCTCA-3’; 23560, R-5’-CCCCTCTTTGCTCCTCTCTT-3’ to detect a point mutation (C > T) by Sanger sequencing. Mice were housed in the barrier facilities of Osaka University under a 12 h light–dark cycle and given ad libitum access to water and food. All procedures were performed according to ARRIVE guidelines and relevant official guidelines under the approval of the Animal Research Committee of Osaka University.

Behavioral overview

All mice used for behavioral testing and weight measurement were 7–8 weeks male littermate progeny of heterozygous Zbtb16 mutant crossings. The following numbers of mice were used for each experiment: wild-type (WT) = 18, KO = 15 for stereotyped behavior, locomotion, open field, elevated plus maze, 3-chamber social interaction, and marble-burying tests; WT = 14, KO = 8 for novel object recognition test. The minimum number of animals for biological replicates was based on previous experiments to enable the detection of a significant difference between groups at P < 0.05. An experimenter blind to genotypes performed all behavioral tests. All behavioral tests were performed between 10:00 to 16:00 h.

Stereotyped behavior test

Mice were placed in a novel home cage where they were habituated for 10 min followed by a 10 min recording period. Time spent and the number of grooming events was manually quantified from recorded movie.

Locomotion test

Mice were placed in a novel chamber (W700 × D700 × H400 mm, #OF-36(M)SQ, Muromachi Kikai Co., Ltd., Tokyo, Japan) and allowed to freely explore for 10 min. Horizontal locomotor activity was measured by ANY-maze behavior tracking software (Stoelting Co., Wood Dale, IL). Habituation to a novel camber was performed the day before the test.

Open field test

Mice were placed in one of the corners of a novel chamber (W700 × D700 × H400 mm, #OF-36(M)SQ, Muromachi Kikai Co., Ltd., Tokyo, Japan) and allowed to freely explore for 10 min. Time spent in the center of the arena (140 × 140 mm) and in all corners of the arena (140 × 140 mm × 4 corners), and locomotor activity were measured and tracked by ANY-maze behavior tracking software.

Elevated plus maze test

Mice were placed in the center of the maze (open arms W54 × D297 mm; closed arms W60 × D300 × H150 mm; Height from floor 400 mm, #EPM-04M, Muromachi Kikai Co., Ltd., Tokyo, Japan) and allowed to freely explore the maze for 5 min. Time and distance in each arm were measured and tracked by ANY-maze behavior tracking software.

Three-chamber social interaction test

The social interaction test consisted of three 5 min trials in the 3-chamber apparatus (W600 × D400 × H220 mm, SC-03M, Muromachi Kikai Co., Ltd., Tokyo, Japan). During the first trial, the mouse was allowed free exploration of the 3-chamber apparatus. Each end chamber contained an empty wire cage (φ90 × H185 mm) with the middle chamber being empty. In the second 5 min trial to examine social novelty, one of the end chambers contained a novel stranger mouse in a wire cage while the opposite end chamber contained an empty wire cage. In the third 5 min trial to examine social cognition, one of the end chambers kept a mouse used in the second trial in the same wire cage as a familiar mouse, while the opposite end chamber contained the other novel stranger mouse in a wire cage. The test mouse was also given a choice between an inanimate cage and a novel stranger mouse in the second trial, and a familiar mouse and a novel stranger mouse in the third trial. Interaction with the targets around a wire cage was tracked and measured by ANY-maze behavior tracking software.

Marble-burying test

Mice were placed in the corner of a novel home cage evenly placed with eighteen marbles and allowed to freely explore for 20 min. After 20 min, the number of marbles buried was recorded. A marble was defined as buried when less than one-third of the marble was visible.

Novel object recognition test

Mice were habituated the day before the test to a chamber (W700 × D700 × H400 mm, #OF-36(M)SQ, Muromachi Kikai Co., Ltd., Tokyo, Japan). On the second day, two same objects were placed on the two opposite corners of a chamber from approximately 50 mm from the closest wall. Then, mice were placed in the corner of a chamber and allowed to freely explore for 10 min. On the third day, one of the objects was replaced with a different shaped object like a novel object. Mice were placed in the corner of a chamber and allowed to freely explore for 10 min. Interaction with novel and familiar objects was tracked and measured by ANY-maze behavior tracking software. The difference score was calculated by subtracting the time exploring the familiar object from the time exploring the novel object. The discrimination ratio was calculated by dividing the time exploring the familiar object by the total time exploring both novel and familiar objects. A positive difference score or discrimination ratio >0.5 indicates that a mouse recognizes the novel object.

Immunohistochemistry

Mouse brains at 7–8 weeks were fixed with 4% PFA in PBS overnight at 4 °C, cryoprotected in 30% sucrose in PBS, then embedded in Tissue-Tek O.C.T. Compound (#4583, Sakura Finetek Japan Co., Ltd., Osaka, Japan) for cryosectioning. Cryosections (20 μm thick) were placed in PBS. Antigen retrieval pretreatment was performed by incubating sections in citrate buffer (10 mM citrate, 0.05% Tween-20, pH 6) at 95 °C for 10 min. Sections were stained with the following primary antibodies: mouse monoclonal anti-NeuN (1:200, #MAB377, Millipore, Billerica, MA), rat monoclonal anti-CTIP2 (1:500, #ab18465, Abcam, Cambridge, UK), rabbit polyclonal anti-TBR1 (1:250, #ab31940, Abcam, Cambridge, UK), mouse monoclonal anti-PDGFRα (CD140a) (1:200, #14-1401-82, Thermo Fisher Scientific, Waltham, MA), mouse monoclonal anti-APC (1:250, #OP80, Merck, Darmstadt, Deutschland), rabbit polyclonal anti-MBP (1:200, #ab40390, Abcam, Cambridge, UK), rabbit polyclonal anti-IBA1 (1:1000, #019-19741, FUJIFILM Wako pure chemical corporation, Osaka, Japan). For fluorescence immunostaining, species-specific antibodies conjugated to Alexa Fluor 488 and/or Alexa Fluor 597 (1:2,000; Invitrogen, Carlsbad, CA) were applied, and cover glasses were mounted with Fluoromount/Plus (#K048, Diagnostic BioSystems, Pleasanton, CA) or ProLong Diamond Antifade Mountant with DAPI (#P-36931 or #P36971, Thermo Fisher Scientific, Waltham, MA) for nuclear staining. DAPI (#11034-56, Nacalai Tesque, Kyoto, Japan) was also used to stain the nucleus. Images were collected using an Olympus microscope and digital camera system (BX53 and DP73, Olympus, Tokyo, Japan) and an All-in-One fluorescence microscope (BZ-X700, KEYENCE Corporation, Osaka, Japan). Cell numbers in the PFC (at bregma 2.22 to 0.86 mm) and cortical thickness in primary somatosensory (S1) cortex (at bregma −1.46 to −1.70 mm) were quantified manually or using KEYENCE analysis software with Hybrid cell count application (KEYENCE Corporation, Osaka, Japan). Myelinated PFC area (at bregma 2.22 to 0.86 mm) was quantified as described previously[22].

Golgi staining

Whole brains collected at 7–8 weeks were subjected to Golgi staining using superGolgi Kit (#003010, Bioenno Tech, LLC, Santa Ana, CA) according to the manufacturer’s instruction. Coronal sections (100 μm thick) were cut using a vibrating blade microtome (VT1000S, Leica Biosystems, Wetzlar, Germany), and mounted on the slides. Images were collected using KEYENCE analysis software with quick full focus (KEYENCE Corporation, Osaka, Japan).

RNA-seq

RNA-seq was performed as a service by Macrogen Japan Corp. (Kyoto, Japan). Briefly, total RNA was extracted from PFC of male mice at 7 weeks with the AllPrep DNA/RNA Mini Kit (#80204, Qiagen, Hilden, Germany) according to the manufacturer’s instruction. RNA integrity number (RIN) of total RNA was quantified by Agilent 2100 Bioanalyzer using Agilent RNA 6000 Pico Kit (#5067-1513, Agilent, Santa Clara, CA). Total RNA with RIN values of ≥8.1 was used for RNA-seq library preparation. mRNA was purified from 500 ng total RNA, and subjected to cDNA library making (fragmentation, first and second strand syntheses, adenylation, ligation, and amplification) by TruSeq Stranded mRNA Library Prep (#20020594, Illumina, San Diego, CA) according to the manufacturer’s instruction. cDNA library quality was quantified by 2100 Bioanalyzer using Agilent High Sensitivity DNA Kit (#5067-4626, Agilent, Santa Clara, CA). The library was sequenced as 101 bp paired-end on Illumina NovaSeq6000.

RNA-seq alignment and quality control

Reads were aligned to the mouse mm10 reference genome using STAR (v2.7.1a)[23]. For each sample, a BAM file including mapped and unmapped reads that spanned splice junctions was produced. Secondary alignment and multi-mapped reads were further removed using in-house scripts. Only uniquely mapped reads were retained for further analyses. Quality control metrics were assessed by the Picard tool (http://broadinstitute.github.io/picard/). Gencode annotation for mm10 (version M21) was used as reference alignment annotation and downstream quantification. Gene level expression was calculated using HTseq (v0.9.1)[24] using intersection-strict mode by exon. Counts were calculated based on protein-coding genes from the annotation file.

Differential expression

Counts were normalized using counts per million reads (CPM). Genes with no reads in either Zbtb16 KO or WT samples were removed. Surrogates variables were calculated using sva package in R[25]. Differential expression analysis was performed in R using linear modeling as following: lm(gene expression–Treatment + nSVs). We estimated log2 fold changes and P-values. P-values were adjusted for multiple comparisons using a Benjamini–Hochberg correction (FDR). Differentially expressed genes were considered for FDR < 0.05. Mouse Gene ID was translated into Human Gene ID using biomaRt package in R[26].

Gene ontology analyses

The functional annotation of differentially expressed and co-expressed genes was performed using GOstats[27]. A Benjamini-Hochberg FDR (FDR < 0.05) was applied as a multiple comparison adjustment.

Network analyses

We carried out weighted gene co-expression network analysis (WGCNA)[28]. Prior to the co-expression analysis, normalized RNA-seq data were residualized and balanced for the nSVs detected using a linear model. A soft-threshold power was automatically calculated to achieve approximate scale-free topology (R2 > 0.85). Networks were constructed with blockwiseConsensusModules function with biweight midcorrelation (bicor). We used corType = bicor, networkType = signed, TOMtype = signed, TOMDenom = mean, maxBlockSize = 16000, mergingThresh = 0.15, minCoreKME = 0.5, minKMEtoStay = 0.6, reassignThreshold = 1e−10, deepSplit = 4, detectCutHeight = 0.999, minModuleSize = 50. The modules were then determined using the dynamic tree-cutting algorithm. A deep split of 4 was used to more aggressively split the data and create more specific modules. Spearman’s rank correlation was used to compute module eigengene-treatment association.

GWAS data and enrichment

We used genome-wide gene-based association analysis implementing MAGMA v1.07[29]. We used the 19346 protein-coding genes from human gencode v19 as background for the gene-based association analysis. SNPs were selected within exonic, intronic, and UTR regions, as well as SNPs within 5 kb up/down-stream the protein-coding gene. SNP association revealed 18988 protein-coding genes with at least one SNP. Gene-based association tests were performed using linkage disequilibrium between SNPs. Benjamini-Hochberg correction was applied and significant enrichment is reported for FDR < 0.05. Summary statistics for GWAS studies on neuropsychiatric disorders and non-brain disorders were downloaded from Psychiatric Genomics Consortium and GIANT Consortium[30-40]. GWAS acronyms were used for the figures (ADHD = attention deficit hyperactivity disorder, AZ = Alzheimer’s disease, ASD = autism spectrum disorder, BD = bipolar disorder, Epilepsy = epilepsy, MDD = major depressive disorder, SCZ = schizophrenia, CognFunc = cognitive functions, EduAtt = educational attainment, Intelligence = Intelligence, BMI = body mass index, CAD = coronary artery disease, OSTEO = osteoporosis).

Gene set enrichment

Gene set enrichment was performed in R using Fisher’s exact test with the following parameters: alternative = “greater”, confidence level = 0.95. We reported Odds Ratio (OR) and Benjamini–Hochberg adjusted P-values (FDR).

Statistical analysis

All behavioral and histological data are represented as means of biological independent experiments with ±standard error of the mean (SEM). Statistical analysis (unpaired t-test) was performed using Prism 7. Asterisks indicate P-values (****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05). P < 0.05 was considered to indicate statistical significance.

Accession number

The NCBI Gene Expression Omnibus (GEO) accession number for the RNA-seq data reported in this manuscript is GSE155424 (token: qzqfauwmdnmprsv).

Results

Zbtb16 KO mice display ASD-like behaviors

Zbtb16 KO mice showed skeletal dysplasia and smaller body size compared with wild-type (WT) mice (Fig. 1a), due to a single nucleotide (lu) mutation in Zbtb16, resulting in a nonsense mutation (p.Arg234*) (Fig. 1b, Supplementary Fig. 1). We first investigated locomotion activity in the mice. Zbtb16 KO mice showed decreased body weight (P < 0.0001) (Fig. 1c), but no difference in normal locomotion activity (P = 0.87) (Fig. 1d).

Fig. 1

Zbtb16 KO mice display social impairment and repetitive behaviors.

a Zbtb16 KO mice display small bodies and hind limb dysplasia. WT: wild-type, KO: Zbtb16 knockout. b A single nucleotide (luxoid) mutation in Zbtb16 results in a nonsense mutation (p.Arg234*). c Weight loss in Zbtb16 KO mice. d No significant difference in locomotion activity in Zbtb16 KO mice. e Representative heatmaps of social novelty behavior in the 3-chamber social interaction test. E: empty, S1: stranger mouse one. f, g Quantification of interaction time (f) and distance during interaction (distance traveled around the targeting cage in white circle) (g) in social novelty section in WT mice. h, i Quantification of interaction time (h) and distance during interaction (i) in social novelty section in Zbtb16 KO mice. j Representative heatmaps of social cognitive behavior in the 3-chamber social interaction test. S2: stranger mouse two, F: familiar mouse. k, l Quantification of interaction time (k) and distance traveled during interaction (l) in social cognition section in WT mice. m, n Quantification of interaction time (m) and distance traveled during interaction (n) in social cognition section in Zbtb16 KO mice. WT mice spent more time and distance interacting with a stranger mouse than a familiar mouse, compared to Zbtb16 KO mice. o Quantification of grooming events. p Quantification of the number of grooming times. q Representative images after a marble-burying test. r Quantification of the number of buried marbles. Zbtb16 KO mice showed an increase in repetitive behaviors like burying marbles. Data are represented as means (±SEM). Asterisks indicate ****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05, unpaired t-test. n = 15–18/condition for weight, locomotion, 3-chamber social interaction, stereotyped behavior, and marble-burying tests.

Since there was no significant difference in the locomotion activity of Zbtb16 KO mice that could confound conducting behavioral analyses (Fig. 1d), we investigated whether Zbtb16 KO mice exhibit behaviors relevant to ASD. We first examined social behaviors using a 3-chamber social interaction test. In the social novelty trial, WT mice preferred a novel mouse to an inanimate empty cage and spent more time (P < 0.0001) and distance (P < 0.0001) interacting with a social target (Fig. 1e–g). On the other hand, Zbtb16 KO mice also preferred a novel mouse over an inanimate empty cage; however, Zbtb16 KO mice spent increased time (P = 0.0090), but not distance (P = 0.37) interacting with a social target (Fig. 1e, i, h). In the social cognition trial, WT mice spent more time (P = 0.0178) and distance (P = 0.0322) interacting with a novel mouse than a familiar mouse (Fig. 1j–l). In contrast, Zbtb16 KO mice spent approximately similar time (P = 0.31) and distance (P = 0.24) for interacting with both a novel mouse and a familiar mouse (Fig. 1j, m, n). These results indicate that Zbtb16 KO mice show decreased social novelty and impaired social cognition.

We next examined repetitive behaviors and found that the number of grooming events (P = 0.0006) and grooming time (P = 0.0014) were significantly increased in Zbtb16 KO mice (Fig. 1o, p). In addition, we found that Zbtb16 KO buried more marbles compared with WT mice (P < 0.0001) (Fig. 1q, r). These results indicate that Zbtb16 KO mice show increased repetitive behaviors.

Together, these results demonstrate that Zbtb16 KO mice display ASD-relevant behaviors.

Zbtb16 KO mice show SCZ-like behaviors

We next investigated whether Zbtb16 KO mice show behaviors relevant to SCZ. Zbtb16 KO mice displayed social impairment, which is similar to one of the negative symptoms in SCZ. It has been reported that impulsive risk-taking behaviors are common in patients with SCZ[41,42]. Thus, we analyzed risk-taking (antianxiety) behaviors in Zbtb16 KO mice using an open field test[43]. We found that Zbtb16 KO mice spent more time (P = 0.0217) and distance (P = 0.0468), but not the number of entries (P = 0.08) in the center of the arena compared to WT mice (Fig. 2a–d). On the other hand, there were no differences in the number of entries (P = 0.69), time (P = 0.27), distance (P = 0.63) in the corners of the arena (Fig. 2a, e–g). These results indicate that Zbtb16 KO mice explore the field more than WT mice, suggesting risk-taking behavior is increased in Zbtb16 KO mice. To clarify the risk-taking behavior of Zbtb16 KO mice, we performed elevated plus maze[43]. Similarly, we found significantly increased time (P = 0.0004) and distance (P = 0.0048), but not the number of entries (P = 0.33) in open arms of the maze in Zbtb16 KO mice compared with WT mice (Fig. 2h–k). In addition, we found significantly decreased time (P = 0.0035), but not the number of entries (P = 0.41) and distance (P = 0.50) in closed arms of the maze in Zbtb16 KO mice (Fig. 2h, l–n). These results demonstrate that Zbtb16 KO mice show risk-taking behaviors.

Fig. 2

Zbtb16 KO mice show risk-taking behaviors due to cognitive impairment.

a Representative heatmaps of anxiety-like behavior in the open field test. b–g Quantification of the number of entries in the center (b), time in the center (c), distance in the center (d), the number of entries in the corner (e), time in the corner (f), and distance in the corner (g) during the open field test. Zbtb16 KO mice spent more time in the center than WT mice. h Representative heatmaps of anxiety-like behavior in the elevated plus-maze test. i–n Quantification of the number of entries in open arm (i), time in open arm (j), distance in open arm (k), the number of entries in the closed arm (l), time in the closed arm (m), and distance in the closed arm (n) during the elevated plus-maze test. Zbtb16 KO mice also spent more time in the open arm than the closed arm compared with WT mice, indicating that impulsive risk-taking behaviors were increased in Zbtb16 KO mice. o Representative heatmaps of learning and memory behaviors in the novel object recognition test. N: Novel mouse, F: familiar mouse. p, q Quantification of difference score (p) and discrimination ratio (q) in the novel object recognition test. Cognitive function was significantly impaired in Zbtb16 KO mice. A positive difference score or discrimination ratio >0.5 means a mouse recognizes the novel object. Data are represented as means (±SEM). Asterisks indicate ***P < 0.001, **P < 0.01, *P < 0.05, unpaired t-test. n = 15–18/condition for open field and elevated plus-maze tests, n = 8–14/condition for a novel object recognition test.

Since cognitive impairment in the novel object recognition test in Zbtb16 KO mice (without age information) has been reported[21], we also investigated the cognitive function of Zbtb16 KO mice in the same test. As we expected, we observed significant reductions in difference score (P = 0.0122) and discriminant ratio (P = 0.0032) of the novel object recognition test in Zbtb16 KO mice (Fig. 2o–q). These results indicate that cognitive function, in particular learning and memory is impaired in Zbtb16 KO mice.

Together, these results demonstrate that Zbtb16 KO mice also display behaviors relevant to SCZ.

Impairments of neocortical deep layer formation in Zbtb16 KO mice

To place the behavioral abnormalities in a biological context, we analyzed the development of the neocortex of Zbtb16 KO mice. Zbtb16 is mainly expressed in the neocortex, striatum, amygdala, hippocampus, midbrain, and cerebellum (Supplementary Fig. 2). Since reduced neocortical area and abnormal deep layer formation of primary motor (M1) cortex in neonatal Zbtb16 KO mice have been reported[21], we examined neocortical thickness and lamination in the S1 cortex of Zbtb16 KO mice at 7–8 weeks. The reduced brain size of Zbtb16 KO mice was observed as previously reported (Fig. 3a). We quantified neocortical thickness with DAPI staining and found a significant reduction of neocortical thickness in the S1 cortex of Zbtb16 KO mice (P = 0.0087) (Fig. 3b, d). Moreover, we identified that L6 (P = 0.0486) was specifically thinner in Zbtb16 KO mice compared to WT mice, but not L5 (P = 0.39) as measured by immunostaining using CTIP2 and TBR1, deep layer markers for L5 and L6, respectively (Fig. 3c, e, f). We also found a significant reduction of TBR1+ cells (P < 0.0001) in the S1 cortex of Zbtb16 KO mice, but not CTIP2+ cells (P = 0.72) (Fig. 3c, g, h). These results demonstrate that deep layer formation, in particular, L6 formation was impaired in Zbtb16 KO mice.

Fig. 3

Impairments of neocortical development in Zbtb16 KO mice.

a Zbtb16 KO mice show smaller brains than WT mice. b Representative images showing thinner neocortical thickness in primary somatosensory (S1) cortex of Zbtb16 KO mice. c Representative fluorescent images of CTIP2 + and TBR1 + neurons in the S1 cortex. d Quantifications of neocortical thickness by DAPI staining. A significant reduction of neocortical thickness was observed in the S1 cortex of Zbtb16 KO mice. e, f Quantification of deep layer thickness by co-immunostaining of CTIP2, a marker for layer 5 (L5) (e) and TBR1, a marker for L6 (f). Neocortical L6 was specifically thinner in Zbtb16 KO mice. g, h Quantification of CTIP2+ (g) and TBR1+ (h) cells in the mouse S1 cortex. The number of TBR1 + L6 neurons were specifically decreased in the S1 cortex of Zbtb16 KO mice. i Representative images of dendritic spines by Golgi staining. j Quantification of dendritic spines of L5 neurons in the mouse PFC (M1 cortex). k Representative fluorescent images of IBA1+ microglia in the PFC. l Quantification of IBA1 + microglia the upper layer (L1-4) of mouse PFC. Significant increases of dendritic spines (j) and IBA1+ microglia (l) were observed in the PFC of Zbtb16 KO mice. m Representative fluorescent images of myelination in the PFC by MBP immunostaining. n, o Representative fluorescent images of PDGFRα+ or APC+ oligodendrocytes in the PFC. p Quantification of the proportion of myelinated neocortex. Immature myelin formation was observed in the PFC of Zbtb16 KO mice. q, r Quantifications of oligodendrocyte in the PFC. Significant reductions of PDGFRα + oligodendrocyte progenitors (q) and APC+ mature oligodendrocytes (r) were observed in the PFC of Zbtb16 KO mice. Insets show a higher magnification of the boxed area depicted in each fluorescent image. Data are represented as means (±SEM). Asterisks indicate ****P < 0.0001, **P < 0.01, *P < 0.05, unpaired t-test. n = 5–6/condition for length and cell count, n = 10 cells/condition from 5–6 mice for the dendritic spine. Scale bars: 1000 μm in a, 300 μm in b, c, k, m–o, 75 μm in insets of c, k, n, o, 10 μm in i.

Increased numbers of dendritic spines and microglia in Zbtb16 KO mice

We further investigated the dendritic spine and microglia, because it has been reported that the numbers of spines and microglia are linked to ASD and SCZ[44,45]. Thus, we quantified the number of dendritic spines of L5 neurons in the PFC of Zbtb16 KO mice by Golgi staining. The number of dendritic spines was significantly increased in Zbtb16 KO mice (P = 0.0178) (Fig. 3i, j). Since microglia are responsible for spine pruning, we analyzed microglia in the upper layer (L1-4) of PFC of Zbtb16 KO mice by immunostaining. The number of IBA1+ microglia was significantly increased in Zbtb16 KO mice (P = 0.0175) (Fig. 3k, l). These results indicate that increased numbers of dendritic spines of L5 neurons and upper layer microglia in Zbtb16 KO mice also underlie functional abnormalities in the PFC.

Immature myelination occurs due to oligodendrocyte loss in Zbtb16 KO mice

Since myelination is associated with ASD and SCZ, as well as cognitive functions[46], we, therefore, examined myelination in the PFC of Zbtb16 KO mice by immunostaining with MBP, a structural component of myelin. The ratio of myelinated area to the total area of the PFC was measured as previously reported[22]. The proportion of myelinated cortex was significantly decreased in Zbtb16 KO mice (P = 0.0304) (Fig. 3m, p). To clarify the cause of the myelination defect, we investigated oligodendrocyte development by immunostaining with oligodendrocyte differentiation markers and found significant reductions of PDGFRα+ immature (P = 0.0269) and APC+ mature (P = 0.0409) oligodendrocytes in the PFC of Zbtb16 KO mice (Fig. 3n, o, q, r). These results suggest that decreased myelination is due to abnormal oligodendrocyte development in Zbtb16 KO mice.

Together, our results suggest that the histological abnormalities in neocortical cytoarchitectures such as L6 formation, dendritic spines, microglia, and myelination may underlie the behavioral deficits of Zbtb16 KO mice.

Zbtb16 regulates neurodevelopmental genes and myelination-associated genes

To understand the molecular mechanisms underlying behavioral and histological phenotypes, we characterized the Zbtb16 transcriptome by RNA-seq. Transcriptome profiles were clearly separated between WT and Zbtb16 KO mice (Fig. 4a, Supplementary Fig. 3). Differential expression analysis of the RNA-seq data uncovered 533 differentially expressed genes (DEGs) (FDR < 0.05) in PFC of Zbtb16 KO mice compared to WT mice (Fig. 4b, Supplementary Table 1).

Fig. 4

The Zbtb16 transcriptome is enriched for ASD, SCZ, and cortical development-related genes.

a Heatmap of differential expressed genes (DEGs) in mouse PFC. b Volcano plot showing Zbtb16 DEGs with the top 20 gene names indicated. FDR < 0.05, Y-axis = −log10(FDR), X-axis = log2(Fold Change). c GO analyses of Zbtb16 DEGs in biological process (BP). Scatterplots represent the top 3 functions in each module. Y-axis = Odds Ratio, X-axis = −log10(p-value). AllGene: all DEGs, DownReg: downregulated DEGs, UpReg: upregulated DEGs. d Zbtb16 DEGs are enriched in neurons and oligodendrocyte. Dot plots of Zbtb16 DEGs showing the cell-type-specific enrichments. Zbtb16 DEGs, in particular, the downregulated DEGs are highly enriched in Slc17a7+, Slc17a6+, and Gad1 + Gad2 + neurons, interneurons, and oligodendrocytes. In contrast, the upregulated DEGs are enriched in Slc17a7+ neurons and endothelial cells. Dot plots represent the Odds Ratio. Colors represent the −log10(FDR). e Dot plots of Zbtb16 DEGs showing the enrichment for disorder-specific human transcriptome. Downregulated Zbtb16 DEGs are enriched in ASD and SCZ downregulated genes. In contrast, upregulated Zbtb16 DEGs are enriched in bipolar disorder and schizophrenia upregulated genes. SCZ: schizophrenia, BD: bipolar disorder, ASD: autism spectrum disorder. f Venn diagram of Zbtb16 DEGs overlapping with SFARI ASD genes. The overlapping top 10 genes are highlighted. n = 5/condition for RNA-seq.

We also performed gene ontology (GO) analysis to identify the functions of DEGs. DEGs are involved in cellular localization, nervous system development, and neurogenesis (Fig. 4c, Supplementary Fig. 4). Interestingly, we found cellular localization, myelination, and axon ensheathment enriched among the downregulated DEGs (Fig. 4c, Supplementary Fig. 4). The negative regulation of the biological process, long-term memory, and regulation of cellular component organization were found in the upregulated DEGs (Fig. 4c). Together, these results suggest that Zbtb16 regulated genes are involved in neurogenesis, myelination, and memory.

To further characterize DEGs, we conducted cell-type-specific enrichment analysis using scRNA-seq data sets (see Methods) and Zbtb16 transcriptome data. We found that the downregulated Zbtb16 DEGs were highly enriched in Slc17a7+, Slc17a6+, and Gad1 + Gad2 + neurons, and interneurons, as well as oligodendrocytes (Fig. 4d, Supplementary Table 2). In contrast, the upregulated Zbtb16 DEGs were enriched in Slc17a7+ neurons and endothelial cells (Fig. 4d, Supplementary Table 2). These results suggest that Zbtb16 plays a role in the neurogenesis of the deep layers and oligodendrogenesis.

Zbtb16-regulated genes are associated with ASD and SCZ

We next examined whether the DEGs are associated with human diseases by enrichment analyses using disorder-specific human transcriptomes[16]. We found that DEGs were enriched for ASD-specific and SCZ-specific downregulated DEGs (OR = 1.9, FDR = 3.1 × 10−03; OR = 1.8 FDR = 2.3 × 10−05, respectively), in particular, downregulated DEGs were highly enriched in ASD-specific and SCZ-specific downregulated DEGs (OR = 3.2, FDR = 3.1 × 10−06; OR = 2.6, FDR = 1.2 × 10−06, respectively) (Fig. 4e). In contrast, upregulated DEGs were enriched in SCZ-specific and bipolar disorder (BD)-specific upregulated DEGs (OR = 2.1, FDR = 0.02; OR = 2.0, FDR = 6.9 × 10−05, respectively) (Fig. 4e). We also investigated how many DEGs overlapped with ASD genes from the SFARI database, and found that 49 DEGs (approximately 10% of the DEGs) (OR = 2.4, FDR = 1.48 × 10−06) overlapped with SFARI ASD genes (Fig. 4f). These analyses indicate that the Zbtb16-regulated transcriptome is related to both ASD and SCZ. These transcriptomic findings point to molecular mechanisms that could underlie the behavioral and histological phenotypes of Zbtb16 KO mice.

Co-expressed gene networks regulated by Zbtb16 are related to ASD and SCZ genes

To identify the individual molecular networks for the unique roles of Zbtb16, we conducted co-expression network analyses using the Zbtb16 transcriptome. Weighted gene co-expression network analysis (WGCNA) identified 35 modules (Supplemental Table S1). Among these modules, the brown, light yellow, and royal blue modules were associated with Zbtb16 genotype (Fig. 5a–d). The brown module had Zbtb16 as one of its hub genes (Fig. 5b). A hub gene is a key modulator in the co-expression network. The majority of hub genes in these 3 modules were also enriched in Zbtb16 DEGs (Supplementary Fig. 5).

Fig. 5

Co-expressed gene networks uncovered gene networks unique to ASD and SCZ.

a Dot plots with standard errors (SEs) showing module eigengene demonstrate the association of the modules detected by parsimony for Zbtb16. Y-axis = rho, X-axis shows the modules highlighted. b–d Zbtb16-specific modules showing the top 10 hub genes ranked by weighted topological overlap values in the brown (b), light yellow (c), and royal blue (d) modules. Zbtb16 DEGs are shown in boldface. e Zbtb16-specific modules showing the enrichment for disorder-specific human transcriptome. The brown module is enriched for ASD and schizophrenia downregulated genes. The light yellow and royal blue modules are enriched in schizophrenia upregulated genes. SCZ: schizophrenia, BD: bipolar disorder, ASD: autism spectrum disorder. f Zbtb16-specific modules are associated with BMI, cognitive traits, ADHD, BD, and epilepsy. Dot plots represent Betas from MAGMA. Dot colors represent the −log10(FDR) from MAGMA. Y-axis shows the acronyms for the GWAS data utilized for this analysis (see Methods). The X-axis shows the modules of this study. ADHD: attention deficit hyperactivity disorder, AZ: Alzheimer’s disease, ASD: autism spectrum disorder, BD: bipolar disorder, Epilepsy: epilepsy, MDD: major depressive disorder, SCZ: schizophrenia, CognFunc: cognitive functions, EduAtt: educational attainment, Intelligence: Intelligence, BMI: body mass index, CAD: coronary artery disease, OSTEO: osteoporosis. n = 5/condition for WGCNA.

We next examined whether the Zbtb16-specific modules show specific enrichment for ASD-related or SCZ-related genes, and found that the brown module was enriched for ASD-specific downregulated DEGs (OR = 2.9, FDR = 6.5 × 10−06) (Fig. 5e). In contrast, the light yellow and royal blue modules were enriched in SCZ-specific upregulated DEGs (OR = 2.4, FDR = 1.2 × 10−05; OR = 2.2, FDR = 1.9 × 10−04, respectively) (Fig. 5e). These results indicate that the brown module is highly enriched for ASD genes, and the other modules are enriched for SCZ genes. Finally, we examined whether the Zbtb16-specific modules were enriched for human genetic variants. The brown module was enriched in educational attainment (EduAtt) and body mass index (BMI) (FDR = 0.03; FDR = 0.02, respectively) (Fig. 5f, Supplementary Table 3). The light yellow module was enriched in ADHD, BD, cognitive functions (CognFun), EduAtt, and intelligence (FDR = 0.04; FDR = 0.02; FDR = 0.02; FDR = 0.003; FDR = 0.01, respectively) (Fig. 5f, Supplementary Table 3). The royal blue module was enriched in BD, epilepsy, and BMI (FDR = 0.01; FDR = 0.02; FDR = 0.05, respectively) (Fig. 5f, Supplementary Table 3). Our results demonstrate that these Zbtb16-specific modules are key gene networks for uncovering the molecular mechanisms of brain development at risk in both ASD and SCZ.

Discussion

In this study, we demonstrate the role of Zbtb16 in ASD-like and SCZ-like behaviors such as social impairment, repetitive behaviors, risk-taking behaviors, and cognitive impairment. We also found histological impairments in the neocortical development of Zbtb16 KO mice. Transcriptomic analyses identified 533 DEGs involved in neocortical maturation such as neurogenesis and myelination. Co-expression network analyses suggest that Zbtb16-specific modules may be distinct with respect to biological pathways underlying ASD and SCZ. Our study proposes the novel role of Zbtb16 in social cognitive behaviors and neocortical development.

In this study, we uncovered behavioral signatures of Zbtb16 KO mice. Zbtb16 KO exhibited abnormal behaviors relevant to the core symptoms of ASD, social impairment, and repetitive behaviors, that are also consistent with the reported symptoms in the brothers with ASD and ZBTB16 mutations[8]. On the other hand, Zbtb16 KO mice showed social and cognitive impairment, as well as risk-taking behaviors, which are all behaviors relevant to SCZ. A previous study has reported that Zbtb16 mutant mice showed an impairment of recognition memory in the novel object recognition test[21]. The PFC plays distinct roles in social behaviors and social cognition by coordinating with the amygdala[47,48], as well as short-term and long-term memory by coordinating with the hippocampus, and decision making[49]. Moreover, the L6 neurons of the visual cortex are also important for processing the object-recognition memory[50].

Risk-taking behavior is a key component in neuropsychiatric disorders such as SCZ[42,51], bipolar disorder[52], and attention deficit hyperactivity disorder (ADHD)[53,54], but also in drug[55] and alcohol[56] abuses. The associations of Zbtb16 with neuropsychiatric disorders are further supported by the gene expression profiles in Zbtb16 KO mice. Interestingly, Zbtb16 is known to be one of opioid-response and alcohol-response genes[57,58]. Thus, it will be interesting to explore the roles of Zbtb16 in risk-taking behavior and cognitive function as future studies. Together, our behavioral findings demonstrate that Zbtb16 plays essential roles in social and cognitive behaviors.

Previous studies have focused on the role of Zbtb16 in the maintenance and proliferation of NSC in embryonic stages[19,20,59] and cortical surface area and M1 cortical thickness in neonatal stages[21]. Compared to those studies, we characterized the role of Zbtb16 in the young adult brain, and found significantly decreased neocortical thickness of L6 and the number of TBR1+ cells in the S1 cortex of Zbtb16 KO mice[21]. Previous work has reported that a smaller cortex and loss of TBR1+ cells in Zbtb16 mutant mice are due to a decrease in proliferating mitotic cell numbers at early embryonic stages[21]. These findings indicate a decrease in L6 thickness of regions at least M1 and S1 cortex, suggesting the decreased neocortical thickness is a common phenotype in Zbtb16 KO mice. However, other brain regions should also be evaluated comprehensively. In addition, the mechanism that affects only deep neurons has not been clarified. Interestingly, it has been reported that the genes associated with ASD are enriched in the deep layers[60,61].

We also found increased numbers of dendritic spines and IBA1+ microglia in Zbtb16 KO mice. The number of spines is well studied in neuropsychiatric disorders. It is frequently reported that spines increase in ASD and decrease in SCZ, but opposite phenotypes have also been reported[45,62]. Changes in the number of spines are not only driven by genetic factors, but pruning by microglia is also an important factor. An increased number of microglia in postmortem brain[63,64] and activated intracerebral microglia[65,66] have been reported in patients with ASD. Microglial autophagy plays an essential role in dendritic spine pruning and social behaviors[67]. In summary, these findings suggest that an increase in the number of dendrite spines in Zbtb16 KO mice could be a responsible factor for the observed impairment of social behavior. However, unlike the general theory that microglia act on synaptic pruning, we found increases in both dendritic spine and microglia in the particular brain region of Zbtb16 KO mice, indicating that an increased the number of IBA1+ microglia in the upper layer, but not in the deep layer (L5-6) (data not shown). Thus, we acknowledge that it could not conclude whether these are the regional or the general PFC phenotypes of Zbtb16 KO mice. As another possibility, an increased dendritic spine of L5 neuron may be the function of Zbtb16 as a cell-autonomous manner in neuronal differentiation. In fact, it is suggested that GO of the royal blue module was involved in the dendritic spine (Supplementary Fig. 6), and the genes such as Slc29a4, Rnf169, Calb1 in the royal blue module were upregulated DEGs in the PFC of Zbtb16 KO mice (Figs. 4b, 5d, Supplementary Table 1), however further investigation is needed to address those questions.

Moreover, we found developmental defects of oligodendrocytes, resulting in impaired neocortical myelination in Zbtb16 KO mice. This myelination-relevant phenotype was predicted by the RNA-seq results. To explain the impairment of neocortical myelination, we identified abnormal oligodendrocyte development via the reduced numbers of PDGFRα+ immature and APC+ mature oligodendrocytes in Zbtb16 KO mice. Zbtb16 is expressed in the OLIG2+ neural progenitor cells in the early developmental stages and regulates oligodendrocyte differentiation by suppressing neurogenesis[59]. Thus, these data suggest that the deletion of Zbtb16 results in abnormal oligodendrogenesis, and eventually causes behavioral abnormalities such as ASD-like and SCZ-like phenotypes. A recent study examined the role of the ASD-relevant gene TCF4 and found that Tcf4 mutant mice show impairments of oligodendrocyte development and myelination, supporting the significance of oligodendrocytes implicating ASD etiology[68]. In addition, reductions in white matter or corpus callosum volumes in the brains of ASD[47,69] and SCZ[70,71] subjects using diffusion tensor imaging (DTI) have been reported. Collectively, our study demonstrates a novel role for Zbtb16 in neocortical development such as abnormalities in deep layer formation, spinogenesis, and oligodendrogenesis, which are similar to both ASD and SCZ pathology.

Here, we describe the Zbtb16 transcriptome in the adult mouse brain, which is also a resource for understanding the role of Zbtb16 in brain development, behavior, and disease. In the RNA-seq results, we identified 533 DEGs in Zbtb16 KO mice involved in essential biological functions such as neurogenesis, myelination, and memory. Overlapping the Zbtb16 DEGs with the list of the SFARI ASD genes also provides ASD-relevant Zbtb16-mediated signaling pathways. Among 533 DEGs, downregulated targets of Zbtb16 were significantly enriched for ASD-specific and SCZ-specific downregulated DEGs. For example, CHRM2 plays a role in the communication between the cortex and hippocampus[72] to modulate cognitive functions such as behavioral flexibility, working memory, and hippocampal plasticity[73]. A recent study has reported that variants in DDHD1 are associated with hereditary spastic paraplegia and ASD[74], and patients show clinical features such as cerebellar impairment, axonal neuropathy, distal sensory loss, and/or mitochondrial impairment[75,76]. HSPA4L is significantly downregulated in lymphoblastoid cells from patients with SCZ[77], however, the function of HSPA4L in the brain is unknown. A recent in silico study has suggested that HSPA4L is upregulated in the corpus callosum of patients with multiple sclerosis, and is associated with myelination and the immune system[78].

Lastly, we identified the ASD-specific (brown) and SCZ-specific (light yellow and royal blue) Zbtb16 modules using WGCNA. In the brown module, the top hub genes were enriched for genes severely downregulated in ASD patients. Moreover, the brown module is enriched for genetic variants associated with EduAtt and BMI, and involved in mitochondrial functions (Supplementary Fig. 6). In contrast, the light yellow and royal blue modules showed specific enrichment for SCZ downregulated genes. At the genetic level, the light yellow module is enriched for genetic variants associated with ADHD, BD, CognFun, EduAtt, and intelligence. The hub genes of the light yellow module are involved in behavior, memory, and synaptic membrane, further confirming the role of Zbtb16 in these synaptic etiologies (Supplementary Fig. 6). Furthermore, the royal blue module is enriched for genetic variants associated with epilepsy and BD. The hub genes of the royal blue module are involved in the dendritic spine, GTP binding, and purine ribonucleoside binding (Supplementary Fig. 6). Consequently, detailed analyses of these Zbtb16 target genes will lead to a deeper understanding that how Zbtb16 regulates the pathological mechanism underlying ASD, SCZ, and other neuropsychiatric disorders, as one of the risk factors.

In summary, our study demonstrates a role for Zbtb16 in social, repetitive, and risk-taking behaviors and cognitive function via neocortical development, particularly deep layer formation and myelination. To the best of our knowledge, this is the first study to directly show the association of Zbtb16 with ASD-like and SCZ-like behaviors, as well as neocortical phenotypes and genomic associations. Our genomics results suggest that ZBTB16 potentially regulates different gene networks that are at risk in ASD or SCZ. However, since only Zbtb16 KO mice were analyzed in this study, the behavioral and brain phenotypes of Zbtb16 heterozygous mice are unknown. Previous studies have shown that there was no statistical difference in the dorsal cortical surface area and cortical length between WT and heterozygous mice[21], but internal structures such as deep layer formation have not been examined. Therefore, we acknowledge that analysis of Zbtb16 heterozygous mice may provide findings that reflect the pathology of patients with ASD or SCZ. Further analysis of the mutations of ZBTB16 (p.Arg440Gln in ASD and p.Lys581* in SCZ) would be the key to unlocking such molecular mechanisms. Presumably, alterations in protein conformation could play a role in the interaction of ZBTB16 with different binding factors, leading to differential downstream transcriptional targets. Further understanding of the differential roles of Zbtb16 in the brain should give rise to novel insights and targets for understanding the molecular mechanisms underlying ASD and SCZ.

Supplementary information

Supplementary Table1

Supplementary Table2

Supplementary Table3