Gene Dosage- and Age-Dependent Differential Transcriptomic Changes in the Prefrontal Cortex of Shank2-Mutant Mice

Shank2 is an abundant postsynaptic scaffolding protein that is known to regulate excitatory synapse assembly and synaptic transmission and has been implicated in various neurodevelopmental disorders, including autism spectrum disorders (ASD). Previous studies on Shank2-mutant mice provided mechanistic insights into their autistic-like phenotypes, but it remains unclear how transcriptomic patterns are changed in brain regions of the mutant mice in age- and gene dosage-dependent manners. To this end, we performed RNA-Seq analyses of the transcripts from the prefrontal cortex (PFC) of heterozygous and homozygous Shank2-mutant mice lacking exons 6 and 7 at juvenile (week 3) and adult (week 12) stages. Juvenile heterozygous Shank2-mutant mice showed upregulation of glutamate synapse-related genes, downregulation of ribosomal and mitochondrial genes, and transcriptomic changes that are opposite to those observed in ASD (anti-ASD) such as upregulation of ASD_down (downregulated in ASD), GABA neuron-related, and oligodendrocyte-related genes. Juvenile homozygous Shank2 mice showed upregulation of chromatin-related genes and transcriptomic changes that are in line with those occurring in ASD (pro-ASD) such as downregulation of ASD_down, GABA neuron-related, and oligodendrocyte-related genes. Adult heterozygous and homozygous Shank2-mutant mice both exhibited downregulation of ribosomal and mitochondrial genes and pro-ASD transcriptomic changes. Therefore, the gene dosage- and age-dependent effects of Shank2 deletions in mice include differential transcriptomic changes across distinct functional contexts, including synapses, chromatin, ribosomes, mitochondria, GABA neurons, and oligodendrocytes.


Tissue Preparation for RNA-Seq Analysis
Brains isolated from Shank2 ± and Shank2 −/− mice (n = 5 brains from five mice for W3-HT, five for W3-HM, four for W12-HT, and four for W12-HM) were chilled in icecold phosphate-buffered saline solution. Mouse brains were not pooled, and samples from individual mouse brains were handled independently. Coronal brain sections (1 mm thickness) were used to dissect out PFC regions (AP + 2.8 to + 1.1 mm; ML −1 to +1 mm). The dissected tissues were kept in RNAlater solution (Ambion) at room temperature for 30 min and then kept at -70 • C. Total RNA concentration was calculated by Quant-IT RiboGreen (Invitrogen, R11490). To determine the integrity of total RNAs, samples were run on the TapeStation RNA screentape (Agilent). Only highquality RNA preparations, with RIN greater than 7.0, were used for RNA library construction. A library was prepared with 1 µg of total RNA for each sample by Illumina TruSeq mRNA Sample Prep kit (Illumina). The first step in the workflow involved purifying the poly-A-containing mRNA molecules using poly-T oligo-attached magnetic beads. Following purification, the mRNA is fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments are copied into first-strand cDNAs using SuperScript II reverse transcriptase (Invitrogen) and random primers, which was followed by second-strand cDNA synthesis using DNA Polymerase I and RNase H. These cDNA fragments then went through an end repair process, the addition of a single "A" base, and then ligation of the indexing adapters. The products were then purified and enriched with PCR to create the final cDNA library. The libraries were quantified using qPCR according to the qPCR Quantification Protocol Guide (KAPA Library Quantification kits for Illumina Sequencing platforms) and qualified using the TapeStation D1000 ScreenTape (Agilent Technologies). Indexed libraries were then submitted to an Illumina NovaSeq (Illumina), and the paired-end (2 × 100 bp) sequencing was performed by Macrogen.

RNA-Seq Analysis
Transcript abundance was estimated with Salmon (v1.1.0) (Patro et al., 2017) in Quasi-mapping-based mode onto the Mus musculus genome (GRCm38) with GC bias correction (-gcBias). Quantified gene-level abundance data was imported to R (v.3.5.3) with the tximport (Soneson et al., 2015) package, and differential gene expression analysis was carried out using R/Bioconductor DEseq2 (v1.22.2) (Love et al., 2014). Some genes were filtered out when their adjusted p values were null value, failing to pass the independent filtering step in DESeq2, and when an HGNC (HUGO Gene Nomenclature Committee) human gene did not exist for a given mouse gene. For independent quality control, the percentage of mapped reads (y-axis) was plotted against the number of mapped reads (x-axis). The scatter plots were drawn using MultiQC tools and the statistics data from Salmon. Principle component analysis (PCA) was performed for the regularized log transform (rlog) of the normalized counts using plotPCA tools implemented in DEseq2. Normalized read counts were computed by dividing the raw read counts by size factors and fitted to a negative binomial distribution. The P values were adjusted for multiple testing for each developmental stage and comparisons with the Benjamini-Hochberg correction. Genes with an adjusted P value of less than 0.05 were considered as differentially expressed. Volcano plots were generated using the R ggplot2 (v.3.1.1) package. The Gene Ontology (GO) enrichment analyses were performed using DAVID software (version 6.8) (Huang da et al., 2009) on the human background. Mouse gene names were converted to human homologs using the Mouse Genome Informatics (MGI) database 1 .

Gene Set Enrichment Analysis
Gene Set Enrichment Analysis (GSEA) 2 (Subramanian et al., 2005) was used to determine whether a priori-defined gene sets would show statistically significant differences in expression between WT and Shank2-mutant mice. Enrichment analysis was performed using GSEAPreranked (gsea-3.0.jar) module on gene set collections downloaded from Molecular Signature Database (MSigDB) v7.0 3 . GSEAPreranked was applied using the list of all genes expressed, ranked by the fold change and multiplied by the inverse of the P value with recommended default settings (1,000 permutations and a classic scoring scheme). The False Discovery Rate (FDR), which adjusts different gene set sizes and multiple testings, was estimated to control the false-positive finding of a given Normalized Enrichment Score (NES) by comparing the tails of the observed and null distributions derived from 1,000 gene set permutations. The gene sets with an FDR of less than 0.05 were considered significantly enriched.
At week 3, eight and 50 DEGs were identified in the PFC of W3-HT and W3-HM mice, respectively, compared to WT mice at 3 weeks ( A numerical comparison of the four DEG sets ( Figure 1G) highlighted that there were age-dependent increases in the total numbers of DEGs in both heterozygous and homozygous Shank2 mice, with a greater increase seen in heterozygous mice. A Venn diagram plotting the DEG numbers showed that there were small overlaps between W3 and W12 DEGs within the HT and HM groups [0.2% (HT) and 5.2% (HM), respectively], compared to the greater overlaps of HT and HM DEGs at W3 (3.6%) and W12 (20.3%) ( Figure 1H).
DAVID analysis of W12-HT transcripts (up-and downregulated together) revealed ribosome-and translationrelated GO terms in the cellular component (CC), biological process (BP), molecular function (MF), and KEGG (Kyoto Encyclopedia of Genes and Genomes) domains  were pooled for these analyses to increase the strength of our analytic results.

Shank2-Mutant Transcriptomes at 3 and 12 Weeks for Biological Functions
We additionally performed GSEA, which uses the entire list of genes with transcriptional changes (ranked by p values) rather than the relatively small number of DEGs above a particular p value or fold change, in order to identify altered biological pathways contributed by modest but coordinate changes in a more reproducible manner (Mootha et al., 2003;Subramanian et al., 2005) (see text footnote 4).
In W3-HT, the transcript pattern of ASD-related genes was strongly opposite that seen in ASD patients and model mice: W3-HT transcripts were negatively and moderately enriched for gene sets that are upregulated in ASD (DEG_Up_Voineagu and Co-Exp_Up_M16_Voineagu) and positively enriched for gene sets that are downregulated in ASD (DEG_Down_Voineagu and Co-Exp_Down_M12_Voineagu) ( Figure 5A). Moreover, W3-HT transcripts were positively enriched for ASD-risk gene sets that are usually downregulated in ASD likely through missense, non-sense, splice-site, frame-shift, and deletion mutations; these included SFARI genes (all), SFARI (high confidence), FMRP targets, DeNovoMis, DeNovoVariants, and AutismKB ( Figure 5A). Therefore, W3-HT transcripts show patterns that are strongly opposite to those observed in ASD (here, called "anti-ASD").
W3-HM transcripts were negatively enriched for gene sets that are upregulated in ASD (DEG_Up_ Voineagu and Co-Exp_Up_M16_ Voineagu) (Figure 5A), similar to the pattern in W3-HT transcripts (anti-ASD). However, W3-HM transcripts were negatively and strongly enriched for gene sets downregulated in ASD (DEG_Down_Voineagu and Co-Exp_Down_M12_ Voineagu) (Figure 5A), and thus also exhibited some "pro-ASD" transcription. W3-HM transcripts were positively enriched for SFARI genes (all), SFARI genes (high confidence), DeNovoMis, DeNovoVariants, and AutismKB but negatively enriched for FMRP targets ( Figure 5A). These results from W3-HM mice partly deviate from the W3-HT pattern and suggest that an increase in the dosage of Shank2 deletion dampens the anti-ASD transcriptomic changes at W3.
The four transcript groups (W3-HT, W3-HM, W12-HT, and W12-HM) were minimally enriched for other brain disorderrelated gene sets, compared with ASD-related/risk gene sets, in DisGeNet 5 (Pinero et al., 2020; Figure 5B). One exception was showing enrichment patterns for ASD-and brain disorder-related gene sets. The entire lists of genes with transcriptional changes ranked by p values were used as inputs. Note that the ASD-related gene sets include those that are upregulated in ASD (DEG_Up_ Voineagu and Co-Exp_Up_M16_ Voineagu), those that are downregulated in ASD (DEG_Down_Voineagu and Co-Exp_Down_M12_ Voineagu), and that ASD-risk genes are likely to be downregulated in ASD [SFARI (all), SFARI (high confidence), FMRP targets, DeNovoMis, DeNovoVariants, and AutismKB] (see Supplementary Table 9 for gene-set details) (n = 5 mice for W3-WT/HT/HM and four mice for W12-WT/HT/HM, FDR > 0.05). (C) GSEA results for all four transcript groups (W3-HT, W3-HM, W12-HT, and W12-HM) showing enrichment patterns for specific cell type-related gene sets (cell-type-specific; see Supplementary  Table 9 for gene-set details). Note that ASD has been associated with downregulation of neuron (glutamate and GABA)-and oligodendrocyte-related genes and upregulation of glia (astrocyte and microglia)-related genes (n = 5 mice for W3-WT/HT/HM and four mice for W12-WT/HT/HM). (D) GSEA results for all four transcript groups (W3-HT, W3-HM, W12-HT, and W12-HM) showing enrichment patterns for cell-type-specific genes validated by single-cell RNA-Seq results of human cortical tissues (single-cell-type-specific; see Supplementary Table 9 for gene-set details) (n = 5 mice for W3-WT/HT/HM and four mice for W12-WT/HT/HM). the relatively strong enrichment of W3-HT, but not the other three transcripts, for genes related to epilepsy, in line with that (1) epilepsy is one of the comorbidities of ASD, (2) epilepsy symptoms were observed in individuals with SHANK2-related genetic variations (although rare) (Leblond et al., 2012), and (3) Shank2 expression is altered in epileptic humans and rat models (Fu et al., 2020). These results collectively suggest that W3-HT transcripts show unique and strong anti-ASD transcriptomic changes, whereas the other three transcript sets (W3-HM, W12-HT, and W12-HM) show both anti-ASD and pro-ASD changes.

GSEA of Shank2-Mutant Transcriptomes for Cell Type-Specific Gene Sets
Autism spectrum disorders has been associated with celltype-specific changes in gene expression, including decreased neuronal gene expression and increased/decreased glial cell expression (increased astrocyte/microglia-related genes and decreased oligodendrocyte-related genes) (Voineagu et al., 2011;Werling et al., 2016). We thus tested if Shank2 transcripts from W3/12-HT/HM mice show any cell-type-specific enrichments, using cell-type-specific gene sets (Kang et al., 2011) and singlenucleus cell-type-specific gene sets (Velmeshev et al., 2019) (see Supplementary Table 9 for further details).
W3-HT transcripts were positively enriched for cortical glutamate neuron-related gene sets identified across different cortical layers (Figure 5C; Kang et al., 2011). Similarly, W3-HT transcripts were largely positively enriched for GABA neuronrelated gene sets, including those for parvalbumin-positive neurons ( Figure 5C; Kang et al., 2011). For glia-related gene sets, W3-HT transcripts were positively and strongly enriched for oligodendrocyte (mature)-, but not for astrocyte-or microgliarelated gene sets.
W3-HM transcripts were also positively enriched for cortical glutamate neuron-related gene sets, similar to W3-HT transcripts; however, they were largely negatively enriched for GABA neuron-related gene sets, including those associated with parvalbumin-positive neurons ( Figure 5C). W3-HM transcripts were negatively enriched for all three glia-related gene sets (astrocyte, microglia, and oligodendrocytes).
W12-HT and W12-HM transcripts were positively enriched for cortical glutamate neuron-related gene sets and largely negatively enriched for GABA neuron-related gene sets ( Figure 5C). In addition, W12-HT and W12-HM transcripts were negatively enriched for glia-related gene sets (astrocytes, microglia, and oligodendrocytes).
Lastly, the four transcript groups (W3-HT, W3-HM, W12-HT, and W12-HM) were tested cell-type-specific genes validated by single-cell RNA-Seq results of human cortical tissues (Velmeshev et al., 2019). The overall enrichment patterns were similar to those described above (Figure 5C) in that W3-HT transcripts were uniquely and strongly anti-ASD, with upregulated neuronal and oligodendrocyte genes, and that the other three transcripts were less strongly anti-ASD, as shown by moderately upregulated neuronal and oligodendrocyte genes and largely downregulated astrocyte/microglia genes ( Figure 5D).

DISCUSSION
In the present study, we analyzed and compared the transcriptomic changes occurring in the PFC of four different Shank2-mutant mouse groups (W3-HT, W3-HM, W12-HT, and W12-HM) using DEG and GSEA analyses. Our goal was to provide insight into the molecular changes induced in the PFC by Shank2 heterozygous and homozygous deletions across juvenile and adult stages in mice.
Our analyses identified differential transcriptomic changes in heterozygous and homozygous Shank2-mutant mice across juvenile and adult stages and interesting biological functions, including those associated with synapse, ribosome, mitochondria, chromatin, GABA neuron, and oligodendrocyte. Before interpreting these results, the following limitations of the study should be pointed out. First, the fact that some biological functions are missing in a subset of the four transcript groups does not necessarily mean that they do not exist. They could be present, but the statistical powers may be low; we used only 4-5 mice per group. Second, the biological functions identified in the present study involving up-or downregulations of the genes are those that work only at the transcript level and are not supported by any functional results. In addition, the transcriptional changes may not reflect functional up-or downregulations but merely reflect compensatory transcriptional changes, which often move toward opposite directions, although additional hints could come from whether the changes are anti-or pro-ASD (see below).
One interesting transcriptomic change was the increase in synapse (particularly excitatory postsynapse)-related gene expression in W3-HT mice, which was not observed in W3-HM, W12-HT, or W12-HM mice (Figures 3, 4). Given that Shank2 is a core component of excitatory postsynaptic protein complexes (Kim and Sheng, 2004;Grabrucker et al., 2011;Sala et al., 2015;Mossa et al., 2017), and the loss of Shank2 exons 6 and 7 suppresses excitatory synaptic transmission mediated by NMDARs in the mPFC, hippocampus, and amygdala (Won et al., 2012;Lee et al., 2015;Chung et al., 2019), and that the W3-HT transcriptome has anti-ASD nature, the upregulation of excitatory synapse-related genes in W3-HT mice is likely to represent indirect changes that compensate for reduced excitatory transmission. However, although W3-HM mice show decreased NMDAR function in the mPFC, the previous studies did not directly test if W3-HT mice show reduced excitatory synaptic transmission in the PFC or hippocampus (Won et al., 2012;Chung et al., 2019).
Although W3-HM mice did not display altered synapserelated gene expression, they did show increased chromatinrelated gene expression. This suggests that a stronger Shank2 deletion (HM relative to HT) converts the "color" of the transcriptomic changes from synapse-to chromatin-related gene expression changes. It is tempting to speculate that broader changes in the expression levels of chromatin remodeling-related genes are required to compensate for the deficits induced by a stronger Shank2 deletion. In line with this, human cortical neurons carrying SHANK2 mutations with excitatory synaptic hyperconnectivity show chromatin-related changes in gene expression (Zaslavsky et al., 2019). In addition, chromatin-related changes have frequently been associated with ASD (De Rubeis et al., 2014) and observed in mice lacking Shank3 (a relative of Shank2) (Zhu et al., 2014;Ma et al., 2018;Qin et al., 2018;Wang et al., 2019;Zhang et al., 2021).
The increased synaptic gene expression observed in W3-HT mice is markedly weakened in W12-HT mice, where only moderate increases in synapse-and chromatin-related gene expression were observed. It is possible that the continued absence of Shank2 expression across juvenile and adult stages dampens the compensatory gene expression changes. In addition, it could be that the period of active synapse development and flexible synapse-related transcriptomic changes comes to an end roughly around the third postnatal week and does not persist into adulthood (Sheng and Kim, 2011). Moreover, the prolonged heterozygous Shank2 deletion might have changed the nature of the compensatory responses from synapse-only changes to those involving both synapse-and chromatin-related gene expressions. Notably, W12-HM mice show transcript patterns that are similar to those of W12-HT mice (moderate increases in synapse-and chromatin-related gene expressions), suggesting that prolonged Shank2 deletion weakens the dosage-dependent effects of Shank2 loss at W12, which were strong at W3 between HT and HM mice. This could be attributable to age-dependent changes or some compensatory changes such as altered expression of other members of the Shank family (i.e., Shank3), as reported previously (Schmeisser et al., 2012).
Another notable transcriptomic change observed in W3-HT mice was the downregulations of ribosome-and mitochondriarelated genes. This change was not observed in W3-HM mice, where only moderate decreases in lipid-and cation transportrelated gene expression levels were observed. This suggests that Shank2 deletions of different dosages induce distinct changes in gene downregulations at W3, similar to the distinct upregulatory patterns seen at W3 (synapse-and chromatin-related genes in W3-HT and W3-HM mice, respectively). Notably, equally strong decreases in ribosome/mitochondria-related gene expression levels were observed in W12-HT mice, suggesting that, unlike the substantial dampening of synapse gene upregulations across W3 and W12, there is no age-dependent dampening of transcriptomic downregulations. In addition, W12-HM mice also showed strong decreases in ribosome-related gene expressions. It thus seems that decreases in ribosome/mitochondria-related gene expression levels represent the strongest and most widespread transcriptomic changes observed in three out of the four transcript sets (W3-HT, W12-HT, and W12-HM, but not W3-HM). Although we named some of the genes associated with ribosome/mitochondrial functions in Results based on the relative contributions to the enrichments, caution should be taken in narrowing down to specific genes because there are multiple gene sets associated with a single biological function (i.e., ribosome), and there are multiple genes in a single gene set that contributed to the enrichments. It is possible that moderate but coordinated transcriptional changes in many genes can lead to strong gene-set enrichments, which might be occurring in ribosome/mitochondria genes in Shank2-mutant mice.
It is unclear how ribosome/translation-related genes are downregulated in Shank2-mutant mice. Protein synthesis is thought to be reciprocally connected with synaptic proteins and functions under the context of ASD-risk gene functions (Santini and Klann, 2014). In addition, Shank2 is known to coordinate excitatory synaptic signaling pathways (Grabrucker et al., 2011;Sheng and Kim, 2011;Sala et al., 2015;Mossa et al., 2017). It is, therefore, possible that suppressed NMDAR or excitatory synaptic functions in Shank2-mutant mice are associated with the decreased ribosome/translation-related transcriptomic activities, perhaps through mTOR signaling, a central regulator of translational initiation and synaptic and neuronal functions implicated in various brain disorders, including ASD (Hoeffer and Klann, 2010;Costa-Mattioli and Monteggia, 2013;Borrie et al., 2017;Saxton and Sabatini, 2017;Switon et al., 2017;Winden et al., 2018;Bagni and Zukin, 2019;Sossin and Costa-Mattioli, 2019) and known to be regulated by synaptic receptors, including NMDARs and metabotropic glutamate receptors (mGluRs) (Hoeffer and Klann, 2010). However, our GSEA on Shank2mutant transcriptomes did not detect any significant changes in mTOR signaling-related gene expressions, although altered mTOR signaling should be directly tested at the protein level using protein phosphorylation/activity as readouts. Intriguingly, bi-allelic but not mono-allelic deletion of Shank2 has been shown to increase Akt phosphorylation without changes in total protein levels during early neuronal differentiation in SH-SY5Y cells, which involved decreased apoptosis and increased cell proliferation (Unsicker et al., 2021), suggesting that Shank2 regulates mTOR signaling downstream of receptor tyrosine kinase activation.
GSEA for ASD-related gene sets indicate that W3-HT transcripts show enrichment patterns that are strongly opposite to those observed in ASD (Figure 5). The gene sets that are downregulated in ASD are negatively enriched (pro-ASD) in W3-HM, W12-HT, and W12-HM, but positively enriched in W3-HT (anti-ASD). In addition, the positive enrichment of W3-HT transcripts for synapse-related gene sets is opposite to the finding that synapse-related gene functions are downregulated in ASD post-mortem cortex (Voineagu et al., 2011). In these contexts, therefore, W3-HT transcripts again show the strongest and most unique anti-ASD patterns.
The four transcript groups are also similarly enriched for the ASD-related cell-type gene sets. W3-HT transcripts show positive enrichment for glutamate neuron-related gene sets (anti-ASD), and similar results are seen in the other three groups (W3-HM, W12-HT, and W12-HM). However, W3-HT transcripts are positively enriched for GABA neuron-related gene sets (anti-ASD), whereas the other three groups show negative enrichments, or less strong positive enrichments, for GABA neuronal genes (pro-ASD). Given the decreased neuronal gene expression in ASD (Voineagu et al., 2011) and mouse models of ASD (Jung et al., 2018;Yook et al., 2019), these results, together with the abovementioned results, indicate that W3-HT transcripts are uniquely and strongly anti-ASD. These early and strong transcriptomic changes might underlie the lack of strong autistic-like behavioral phenotypes in heterozygous Shank2-mutant mice, although hyperactivity was observed in these mice (Won et al., 2012). In addition, these results further suggest that decreased GABA neuron-related gene expression and consequent GABA neuronal dysfunction might promote ASD-like phenotypes in Shank2-mutant mice. Consistent with this possibility, parvalbumin-positive neuron-specific deletion of Shank2 (exon 6 + 7) leads to moderate hyperactivity, enhanced self-grooming, and suppressed seizure susceptibility in mice . In addition, promoting GABAergic synaptic transmission rescues cognitive deficits observed in adult homozygous Shank2-mutant mice (Lim et al., 2017).
W3-HT transcripts show strong positive enrichment for oligodendrocyte-related gene sets, whereas the other three groups show negative, or only moderate, enrichments. Given that downregulation of oligodendrocyte-related gene sets, but not those related to astrocytes or microglia, has been implicated in ASD (Voineagu et al., 2011), our present results suggest that W3-HT transcripts but not the other three are anti-ASD patterns. Whether this difference differentially affects the density and function of oligodendrocytes in Shank2-mutant mice remains to be determined, although abnormalities in oligodendrocytes and axonal myelination have been associated with ASD (Graciarena et al., 2018;Galvez-Contreras et al., 2020;Kawamura et al., 2020a,b).
In summary, our study demonstrates that heterozygous and homozygous Shank2 (exons 6 + 7) deletions in mice at W3 and W12 induce differential transcriptomic changes, and that these are distinct at W3 but become more similar at W12. Our study also identifies distinct functional groups of genes whose expression changes are associated with Shank2 deletions in mice, such as those related to synapses, chromatin, ribosomes, mitochondria, GABA neurons, and oligodendrocytes.

DATA AVAILABILITY STATEMENT
The RNA-Seq datasets presented in this study can be found in the GEO (Gene Expression Omnibus) repository under Accession Number (GSE168211).

ETHICS STATEMENT
The animal study was reviewed and approved by the committee of animal research at KAIST (KA2020-99).

AUTHOR CONTRIBUTIONS
SL, EL, HJ, and HK performed RNA-Seq analysis. EK and EL designed research and wrote the manuscript. All authors contributed to the article and approved the submitted version.  subdomains (CC, BP, and MF) for W3-HT, W3-HM, W12-HT, and W12-HM transcript groups, showing information on individual genes contributed to the enrichments such as original rank and running enrichment score.