Clustering the autisms using glutamate synapse protein interaction networks from cortical and hippocampal tissue of seven mouse models

Abstract

Background

Autism spectrum disorders (ASDs) are a heterogeneous group of behaviorally defined disorders and are associated with hundreds of rare genetic mutations and several environmental risk factors. Mouse models of specific risk factors have been successful in identifying molecular mechanisms associated with a given factor. However, comparisons among different models to elucidate underlying common pathways or to define clusters of biologically relevant disease subtypes have been complicated by different methodological approaches or different brain regions examined by the labs that developed each model. Here, we use a novel proteomic technique, quantitative multiplex co-immunoprecipitation or QMI, to make a series of identical measurements of a synaptic protein interaction network in seven different animal models. We aim to identify molecular disruptions that are common to multiple models.

Methods

QMI was performed on 92 hippocampal and cortical samples taken from seven mouse models of ASD: Shank3B, Shank3Δex4-9, Ube3a^2xTG, TSC2, FMR1, and CNTNAP2 mutants, as well as E12.5 VPA (maternal valproic acid injection on day 12.5 post-conception). The QMI panel targeted a network of 16 interacting, ASD-linked, synaptic proteins, probing 240 potential co-associations. A custom non-parametric statistical test was used to call significant differences between ASD models and littermate controls, and Hierarchical Clustering by Principal Components was used to cluster the models using mean log₂ fold change values.

Results

Each model displayed a unique set of disrupted interactions, but some interactions were disrupted in multiple models. These tended to be interactions that are known to change with synaptic activity. Clustering revealed potential relationships among models and suggested deficits in AKT signaling in Ube3a^2xTG mice, which were confirmed by phospho-western blots.

Conclusions

These data highlight the great heterogeneity among models, but suggest that high-dimensional measures of a synaptic protein network may allow differentiation of subtypes of ASD with shared molecular pathology.

Background

As the incidence of autism spectrum disorder (ASD) has climbed over the past decades to 1 in 59 children [1], next-generation sequencing studies have described likely causative mutations in hundreds of genes, each accounting for < 0.1–1% of the total autistic population [2,3,4]. Additional factors such as maternal immune activation [5], maternal anti-brain antibodies [6], chemical exposures [7], and polygenetic inheritance of a susceptible genetic background [8] all likely contribute to the development of ASD on an individual-by-individual basis. Thus, much like cancer, ASD is an individually rare, collectively common disorder with a shared diagnostic phenotype: reduced interest in social interaction, reduced communication, and increased stereotyped or repetitive interests and behaviors [9].

The fact that ASD is a diagnostic entity, with a common set of behavioral impairments shared among patients, has led to the widespread hypothesis that other disease mechanisms must also be shared among patients at the level of anatomy [10], neural circuits [11], genetic networks [12, 13], or molecular pathways [14]. Along these lines, a few clear themes have emerged from combining diverse lines of evidence: the immune system is likely involved, with immune-mediated risk factors (reviewed in [15]), and abnormal peripheral [16] and central ([17, 18], but see [19]) inflammatory phenotypes present. Gene regulatory pathways are clearly implicated by genetic studies, as a large percentage of ASD-linked genes are transcription factors, chromatin remodelers, or translational regulators [4, 12]. Synaptic proteins have also been implicated by genetic studies, and by the fact that one unifying feature of animal models of ASD has been disrupted synaptic transmission [20] (although note that the specific nature of the disruption varies greatly between models, or even between brain regions in the same model, discussed below). Recently, unifying theories of ASD have proposed that disruptions to activity-dependent, homeostatic neuronal processes are an underlying characteristic of ASDs [21, 22]. Indeed, diverse ASD-linked genes can disrupt the complex molecular circuitry that translates synaptic ion currents into intracellular signal transduction cascades, traffics those messages to sites of translation and transcription, and converts protein-level modifications into long-term changes in gene expression.

Despite these hints at convergent mechanisms, heterogeneity is still the dominant theme when comparing different autism types [23], or even when comparing genetically similar autisms. The prototypical example is the gene Shank3, responsible for Phelan-McDermid syndrome-associated autism and implicated in ~ 1% of total ASD cases [24]. Shank3 encodes multiple alternatively spliced protein variants (at least six), which each contain different combinations of protein-interaction-mediating domains. No fewer than 13 different mutant mouse lines have been reported thus far, which disrupt different exons of Shank3. While the majority of lines show deficits in social (nine lines), repetitive (nine lines), or vocalization (four lines) behavior, each line shows a different combination of behavioral and molecular deficits, depending on which Shank3 isoforms are disrupted (reviewed in [24]). For example, in a complete knockout line, reducing mGluR5 activity normalized repetitive grooming [25], while in an exon 11 deletion line, enhancing mGluR5 activity rescued abnormal grooming [26]. Similarly, at the level of electrophysiology, reduced striatal mEPSP amplitude and frequency has been reported in adult Shank3B^−/− striatum [27], reduced amplitude but increased frequency in Shank3Δex4–9^+/− hippocampus [28], increased mEPSP frequency and amplitude in p14 Shank3B^−/− striatum [29], and increased activity in p14 Shank3B^−/− cortex [29]. Thus, even models targeting the same gene display different phenotypes dependent on mutation type, age, and brain region. For the majority of ASD genes, only a single model (typically a complete knockout) has been published, and the ages and brain regions targeted differ between labs, complicating attempts at directly comparing pathology between published studies.

This study was designed to make a series of identical, directly comparable molecular measurements in several mouse models of autism, in order to address the question of molecular convergence among models. We compared measurements of synaptic proteins in two brain regions (frontal cortex and hippocampus), in age-and-sex-matched adult (postnatal day 60) animals from six genetic and one environmental model of ASD. We used a novel proteomic technique, Quantitative multiplex co-immunoprecipitation (QMI), that compares the abundance of, and interactions among, a panel of native proteins in mutant animals vs. a matched wildytype littermate control. In QMI, protein complexes are immunoprecipitated onto 5 um polystyrene latex beads and probed with fluorophore-coupled antibodies to quantitatively measure the amount of proteins in shared complexes. The resulting fluorescent signals are read on a flow cytometer, and raw abundance measures are normalized to wildtype controls run on the same plate to cancel out batch effects; only fold-change values compared to control are reported [30].

QMI is a candidate-based approach that targets carefully selected networks of interacting proteins. The high-dimensional data produced is linear over a large dynamic range and is several-fold more sensitive than traditional Western blotting techniques [31]. We used a previously published QMI panel that targets 16 synaptic proteins and measures 240 binary proteins in shared complexes by exposed surface epitopes (PiSCES). This panel consists of ASD-linked proteins that are known to physically interact at the synapse [32]. In each mouse model, we identified a unique combination of disrupted PiSCES, with occasional overlap of disruptions that were common to multiple models. We then clustered the data by model and brain region to reveal possible higher-level relationships among the seven animal models, and we confirmed a previously unreported molecular deficit in one model that was predicted by our clustering. Our approach has the potential to identify unexpected commonalities among genetic autisms and to suggest novel treatments based on shared molecular pathology.

Methods

Animal models

The specific identity of all mouse strains used is shown in Table 1. Littermate mice were co-housed in groups of 2–5 under standard laboratory conditions. At 60 days of age, mice were deeply anesthetized with isofluorane, decapitated, and brains were removed. We chose day 60 because we wanted to focus on adult animals, since the majority of animal models have been behaviorally tested as adults (see Table 1). We used two males and two females for each genotype to focus our study on robust, non-sex-dependent effects, except FMR1^−/y mice and controls, which were all male since FMR1 is an X-linked gene. We used frontal cortex for all models and hippocampus for some models because these two brain regions have been frequently analyzed in electrophysiology and biochemical studies of ASD models (see Table 1). For frontal cortex, the rostral 3 mm of cortex was cut with a razor blade in a metal brain mold, making sure not to include any striatal tissue in the section, and the olfactory bulb was removed; for hippocampus, bilateral hippocampi were removed with curved forceps. Tissue was frozen in liquid nitrogen and stored at − 80 until homogenization. All work was performed under an approved animal protocol at Seattle Children’s Research Institute (#15580).

Table 1 Animal models used in this study

VPA mice were prepared at McMaster University in compliance with standards of the Canadian Council on Animal Care and with approval from the McMaster University Animal Research Ethics Board. CD-1 female mice were mated until a sperm plug was detected (E0). On day 12.5 after conception (E12.5), pregnant mice received a single intraperitoneal (i.p.) injection of 500 mg/kg sodium valproate (VPA; Sigma, Oakville, ON, Canada) dissolved in 0.9% NaCl solution, while controls were injected with only saline. E12.5 was chosen to match previous reports from our group and others (see [32]). Pups were weaned on postnatal day (PD) 21 and subjected to behavioral assays (three-chamber sociability, elevated plus maze, and marble-burying assays for social behavior, anxiety, and repetitive behavior, respectively) on PDs 29–34. Animals were killed by decapitation on PD 35, and brains were rapidly dissected and stored at − 80 °C.

QMI analysis

Tissue was homogenized in 0.32 M Sucrose in HEPES buffer, pH 7.4 with Sigma protease (Cat # P8340) and phosphatase (Cat # P5726) inhibitors (Sigma Aldrich), using 12 strokes of a glass-Teflon homogenizer. Samples were spun at 1000×g for 5 min to pellet membranes, then spun at 10,000×g for 15 min to pellet P2 synaptosomes. Synaptosomes were solubilized on ice in 200 ul lysis buffer (150 mM NaCl, 50 mM Tris (pH 7.4), 1% NP-40, 10 mM NaF, 2 mM sodium orthovanadate + protease/phosphatase inhibitor cocktails [Sigma]) for 15 min, spun at 4 °C at 10,000×g for 15 min to remove insoluble material, and protein concentration was measured by BCA assay (Thermo-Fisher).

QMI beads (Luminex) were prepared as previously described [32], with each bead color-class coupled to a distinct immunoprecipitating antibody, as shown in Table 2. Equal amounts of protein from each matched pair of animals (transgenic vs. wild type littermate or VPA- vs. saline-treated control) were incubated with QMI beads overnight at 4 °C, with constant rotation. Beads from each sample were then distributed into 32 wells of a 96-well plate, approximately 250 beads of each class per well, and each of 16 probe antibodies was added, in duplicate, to individual wells. Beads were then washed with ice-cold Fly-P buffer [50 mM Tris (pH 7.4), 100 mM NaCl, 1% bovine serum albumin, and 0.01% sodium azide], incubated for 30 min with streptavidin-PE (1:200, BioLegend), washed again, and read on a custom refrigerated Bioplex 200 flow cytometer (BioRad), which recorded the bead classification (corresponding to IP’d protein, X) and PE fluorescence (corresponding to the amount of probe antibody target protein, Y) of each bead. An above-background reading for IP:X Probe:Y indicates the occurrence of a protein complex containing both X and Y [33].

Table 2 QMI targets, autism linkage, and antibody information

Data analysis

Data were exported in .xml files containing all data on a bead-by-bead, well-by-well basis. A custom Javascript was written to generate histograms showing bead distributions for a given bead class in a given well and to extract the median fluorescent intensity of each bead class in each well for export to Excel and R (faculty.washington.edu/seps/program). A custom MatLab script, “Adaptive Non-parametric statistical test with an adjustable alpha Cutoff” (abbreviated ANC), previously described in detail [30], was used to identify interactions that changed significantly in > 70% of experiments; these interactions are referred to as “hits.” ANC first uses a K-S test to compare histogram distributions of technical replicates to both discard duplicate wells that are significantly different from each other (presumed manual error) and to adjust the alpha value based on technical error. K-S test results from comparisons between an experimental sample and a matched control are then corrected for multiple comparisons and technical errors to obtain a final p value. “Hits” were interactions with p < 0.05. Please see [30] for details. Prior QMI analysis in both T cells [30] and neural tissue [32] found that N of four biological replicates are sufficient to produce a consistent number of significant hits, so an N of at least four matched pairs was used. To eliminate batch effects due to both technical and biological variation, we limit comparisons to ASD model animals and co-housed, littermate controls euthanized on the same day and run on the same assay plate; ANC statistics are therefore based on consistent differences in paired comparisons for N = 4 experiments (each run with technical replicates). Workflow and examples of smoothed histograms are shown in Fig. 1.

Fig. 1

Workflow. a 3-mm sections of frontal cortex or bilateral hippocampi were collected from matched pairs of wildtype and mutant littermates. b P2 fractions were prepared to enrich for synaptic proteins. Shown here is typical enrichment of Homer1, PSD95, and NMDAR1 in P2 fractions, compared to equal amounts of total protein (by BCA assay) from brain homogenate (HO), P1 membrane pellet, and S2 soluble protein. c A panel of IP beads, each conjugated to a different antibody, is incubated with lysate, probed with fluorophore-conjugated antibodies, and read on a flow cytometer. d Known protein-protein interactions among the targeted protein network, in mouse, from the BioGRID database. Red lines indicate IP-western interactions, black lines IP-mass spectrometry. e, f Example histograms and corresponding node-edge visualizations. e Reduced IP: Shank3 Probe: Shank3 in a Shank3B^−/− animal. Blue loop on Shank3 indicates a negative log₂FC of an IP_Probe for the same target. f Increased Homer_PSD95 in VPA cortex lysate (red) and matched wildtype littermate control (black). Red line between nodes indicates positive log₂FC of an interaction

Data matrices for each matched pair were exported from Java to Excel. For each matrix position, we divided the median fluorescence value (of the two technical replicates) of each ASD model animal by its wildtype littermate control and log₂-transformed the result. Then, log₂ fold change (log₂FC) matrices from N = 4 experiments were averaged to generate a single mean log₂FC matrix per genotype/tissue type, shown in Additional file 1: Table S1. ANC significant hits were identified and imported into Cytoscape for visualization. Significant interactions are represented by an edge connecting two protein nodes; the color and width of the edge corresponds to the direction (red = up, blue = down) and magnitude of the change (Fig. 1e, f). Changes in protein abundance (IP probe of same target protein) are represented by loops.

To cluster samples by shared ANC-significant hits, we used the hclust(dist()) and heatmap.2 functions in R. To cluster samples by average fold-change matrices shown in Additional file 1: Table S1, we first performed principal component analysis to reduce noise due to nonspecific background fluctuations using the “PCA” function in the “FactomineR” package for R; then we used the “HCPC” function in the same package to cluster genotypes/tissue types by principal components. To test the robustness of clustering, we used the “pvclust” function in R. All options were used in the default settings.

Western blots were run on cortical tissue using standard protocols. Briefly, cortical P2 fractions were lysed in lysis buffer, protein concentrations were normalized using BCA assays, equal amounts of protein were loaded into each well and run at 110 V. Protein was transferred onto PVDF membranes, blocked with 5% milk in TBS-T, primary antibodies were incubated overnight at 4C, followed by washes, species-specific secondary antibody incubation (anti-mouse or rabbit, 1:10,000, Jackson Immunoresearch), and luminol detection (Pierce Femto reagent). Antibodies used (all 1:1000 dilutions) the following: Ube3a clone E-4 (Santa Cruz), pAKTs473 clone D9E, pAKTt308 clone 244F9, panAKT clone 40D4, pMTOR polyclonal Cat #2971, and pS6 clone D57.2.2 (all from Cell Signaling).

Results

Additional file 1: Table S1 shows the median log₂ fold change values for the complete dataset (N = 92 samples from 56 animals; 7 ASD models, with 4 ASD model animals and 4 WT controls per group, except N = 6 for Shank3B hippocampus; some animals contributed both cortical and hippocampal tissue). Numbers in bold case indicate fold changes > 1.19 or < 0.84 (which corresponds to ± 0.25 in log₂ scale), while red highlighting indicates that a value was statistically significant by ANC statistical analysis. Note that while some ANC-significant values are smaller than ± 0.25, indicating a small but high-confidence change, several bolded cells are not ANC-significant due to biological or technical variation and the stringent requirements of our statistical test. Below, we first focus our analysis on only significant ANC hits, then we perform inter-model comparisons using the entire data matrix to attempt to cluster models into biologically relevant groups.

Cortex

Overall, we found 32 statistically significant differences across the 7 mouse models (Additional file 1: Table S1 Sheet 2, and Fig. 2). Of 240 total IP_Probe combinations, 9 proteins (IP probe for the same target) and 18 proteins in shared complexes (PiSCES—IP probe for different proteins) showed differences in abundance across the 7 models. Four PiSCES (Homer1_PSD95, Homer_NMDAR1, SynGAP_PSD95 and NL3_FYN) and three abundance measures (FYN, SynGAP and PSD95) were significantly different in multiple models, while the remaining differences were unique to a single model.

Fig. 2

Cortical QMI Diagrams for seven autism models. Edges indicate ANC-significant (p < 0.05) changes in the connected nodes; red = increased, blue = decreased in mutant/wildtype comparisons. Node thickness and color indicates the magnitude of the change. a Shank3B^−/−, N = 4 pairs. b Shank3Δex4–9^+/−, N = 4 pairs. c Ube3a^2xTG, N = 4 pairs. d E12.5 VPA exposed animals, N = 4 pairs. e TSC2^+/−, N = 4 pairs. f FMR1^−/y, N = 4 pairs. g Cntnap2^−/−, N = 4 pairs

Shank3B^−/− animals (Fig. 2a) showed a reduction in Shank3 protein levels and an increased co-association of Homer with PSD95 and NMDAR1. This is counterintuitive, since Homer-PSD interactions are likely mediated by Shank proteins, and may be expected to be reduced in Shank3 animals. These results may reflect changes in Shank1/2 vs. 3 scaffolding, or an increase in these activity-labile interactions [32] may be downstream of reduced synaptic activity in mutant animals [27]. Interestingly, these same two interactions were changed, but in the opposite direction, in FMR^−/y mice (see below), which have increased basal activity, potentially supporting an activity-dependent mechanism (see discussion). We also observed a small increase in SynGAP_PSD95 (also activity-labile) and a small decrease in FYN_NL3, interactions that were both observed in the opposite direction in the TSC2^+/− model.

Shank3Δex4–9^+/− animals (Fig. 2b) were tested as heterozygotes because the heterozygotes showed abnormal behavior in the original publication [28] and more accurately represent the human condition, a heterozygous deletion/mutation. Consistent with only moderately reduced Shank3 levels, a modest reduction in Shank3 (log₂FC = − 0.227) was not significant by ANC. Remarkably, there was no overlap in significant hits with Shank3B knockouts; in fact, SynGAP_PSD95 was significant in the opposite direction in the two models. Similar to the TSC2^+/− animals, Shank3Δex4–9^+/− animals showed significant decreases in SynGAP_SynGAP and SynGAP_PSD95, in addition to reductions in SynGAP_Homer1A and SynGAP_NMDAR1 that were unique to this model. Finally, a moderate increase in FYN_Shank3 was observed. While published electrophysiology revealed reduced excitatory transmission in both the Shank3B^−/− and Shank3Δex4–9^+/− animals [27, 28], the experiments were performed in different brain areas (striatum and hippocampus, respectively) and showed small but important differences, such as reduced vs. increased miniature EPSP frequency, respectively. In summary, while QMI data from the Shank3Δex4–9^+/− animals highlight reduced SynGAP associations with NMDARs and scaffolds, the Shank3B animals show differences in Homer-PSD-NMDAR complexes but no changes in SynGAP. These data suggest that the molecular deficits in the two animal models may be quite different, consistent with the different isoforms that are affected in the two models [24].

Ube3a^2xTg mice (Fig. 2c) showed an expected increase in the amount of Ube3a and reduced co-association between GluR1_PSD95 and NL3_GluR2. Prior work in the cortex of Ube3a animals showed reduced glutamatergic transmission [34] and reduced scaffolding of GluRs is therefore consistent with prior observations.

E12.5 VPA mice (Fig. 2d) are the only non-genetic model analyzed here. Mice were generated by injection of VPA on E12.5, and the efficacy of the treatment was confirmed by behavioral testing (as in [32]) of the adult offspring before dissection and QMI analysis. We observed a large increase in the amount of co-associated Homer_PSD95. In all other models, the amount of Homer_NMDAR1 correlated with Homer_PSD95, and VPA mice were trending towards an increase in this interaction as well (log₂FC = 0.48, NS). In addition, levels of Fyn were increased, PSD95 were decreased, and interactions between SAP97_Homer1A and Shank3_NL3 were increased. Prior reports in VPA-treated rat cortex showed enhanced NMDAR-mediated synaptic currents and enhanced LTP [35], consistent with the observed increase in Homer-NMDAR scaffolding. Decreased levels of PSD95 have also been reported by Western blotting [52].

TSC2^+/− mice (Fig. 2e) showed large reductions in the abundance of SynGAP and SynGAP_PSD95. The mTOR activator Rheb (ras homolog enriched in brain), which is directly suppressed by the TSC1/2 complex, is activated by SynGAP following NMDAR stimulation [36], so the reduction of SynGAP may be a homeostatic response to chronically activated Rheb. Reduced levels of NMDAR2B were also observed, along with increased abundance of complexes containing SAP97_mGluR5, Homer1_SAP97, and Fyn_NL3. Taken together, these data indicate reduced NMDAR2B and SynGAP expression, abnormal scaffolding of mGluR5 to Sap97, and abnormal FYN signaling, which could contribute to the altered LTP phenotype reported in the hippocampus of TSC mice [37].

Fragile X mice (Fig. 2f) showed increased abundance of complexes containing NMDAR2A_NL3 and NMDAR2B_NL3. Both NMDA receptors [38] and Neuroligins [39] bind PSD95, which could mediate this observed interaction. FragileX mice also showed reduced complexes with Homer1_PSD95 and Homer1_NMDAR1, demonstrating disrupted Homer-Shank-PSD95-NMDAR complexes, consistent with previous reports [40,41,42]. These activity-dependent interactions were also significant hits in the Shank3B and E12.5 VPA models, but in the opposite direction, possibly reflecting hyper- vs. hypo-activity of cortical neurons in these models. Finally, reduced levels of PI3K and increased Fyn were detected, consistent with disrupted kinase cascades downstream of mGluR5 in FMR1 mice [43,44,45].

CNTNAP 2 KO mice (Fig. 2g) showed the greatest number of ANC hits (7) as well as many large but non-significant changes. The abundance of GluR2 was reduced, accompanied by reduced GluR2_Shank3, NMDAR1_Shank3, and NMDAR2B_Homer1A, consistent with reduced scaffolding and expression of glutamate receptors. In addition, the detected levels of Sap97 and SynGAP were reduced, while Fyn was increased, a change also observed in the Fragile X model. While the CNTNAP2 gene product CASPR is known to cluster at the nodes of Ranvier following myelination [46], acute CASPR knockdown acts cell-autonomously to reduce both AMPA and NMDA-mediated EPSPs [47], congruent with reduced NMDAR and AMPAR levels and scaffolding observed here.

Hippocampus

In four models, we also isolated P2 fractions from the hippocampi of the same animals that supplied cortical tissue. We identified 45 statistically significant differences across the 4 mouse models, with the majority of differences, 21, found in the VPA hippocampus (Fig. 3). Nine proteins showed differences in protein abundance (IP probe for the same target), and 32 protein interactions showed differences (IP probe for different targets). However, only three complexes (Homer_mGluR5, Sap97_NMDAR1, and SynGAP_NMDAR2A) and 1 abundance measure (PSD95) were detected in multiple models, while the remainder was unique to a single model. Below, we describe the findings from each model, compared with prior data from the cortex of the same model.

Fig. 3

Hippocampal QMI diagrams for four autism models. Edges indicate ANC-significant (p < 0.05) changes in the connected nodes; red = increased, blue = decreased in mutant/wildtype comparisons. Node thickness and color indicates the magnitude of the change. a Shank3B^−/−, N = 6 pairs. b Shank3Δex4–9^+/−, N = 4 pairs. c Ube3a^2xTG, N = 4 pairs. d E12.5 VPA exposed animals, N = 4 pairs

Shank3B^−/− hippocampal tissue (Fig. 3a) showed an increase in PSD95 levels and a decrease in FYN levels, neither of which were observed in cortical tissue. A decrease in NMDAR1_mGluR5 likely reflects disrupted scaffolding linking the two receptor types via PSD95/Shank3/Homer linkages [24]. We observed an increase in Homer1_SYNGAP, PSD95_GluR1, and, counter-intuitively, PSD95_Shank3. The latter interaction may reflect elevated expression of an alternative isoform of Shank3 that lacks the PDZ domains in complex with PSD95, possibly mediated via another protein such as Homer. Besides this interaction, no major changes in Shank3 were detected, likely due to the fact that very little Shank3 was detected from Shank3 IPs or probes, consistent with low hippocampal Shank3 expression. We are unable to relate these changes to known electrophysiological abnormalities in these animals, since to our knowledge, hippocampal electrophysiology has not been reported in this model.

Shank3Δex4–9^+/− animals (Fig. 3b) showed reduced levels of SynGAP, consistent with cortical tissue from these animals. Interactions involving SynGAP_NMDAR2A were reduced, while SynGAP_GluR2 were increased. Complexes containing NL3_NMDAR2A and _FYN were reduced. Complexes containing PI3K_Sap97, _GluR2, and _SynGAP were all increased. Hippocampal electrophysiology in this model indicated reduced basal AMPA-mediated transmission, and a failure of hippocampal LTP that was correlated with failure to maintain spine expansion following a tetanizing stimulation [28]. Our results indicate that SynGAP, a critical mediator of signal transduction downstream of NMDARs [36, 48], is dysregulated in hippocampal tissue prior to any type of stimulation. Further, changes in FYN and PI3K suggest downstream disruption of signaling cascades.

Ube3a^2xTG hippocampus (Fig. 3c) showed the expected increase in Ube3a expression, the only change that was consistent between hippocampus and cortex. A reduction in mGluR5 levels, Homer_mGluR5, and Homer_GlurR2 suggest reduced Homer-mediated scaffolding. Ube3A_Homer interactions were strongly increased, although the significance of this increase is unclear since Ube3a has not been documented to bind directly to or ubiquinate/degrade Homer proteins. The amount of PSD95_Shank3 was increased, as was SAP97_NMDAR1 and SAP97_Shank3. Finally, the amount of Homer1A was increased. These data demonstrate complex changes in scaffolding of AMPA, NMDA, and metabotropic glutamate receptors mediated by both Homer and DLG scaffolds in the Ube3a^2xTG animal. Hippocampal electrophysiology has not been reported in these animals, although LTP disruptions due to lack of small conductance potassium channel 2 (SK2) channel regulation have been reported in the Ube3a knockout animal [49].

VPA hippocampus yielded 21 significant QMI hits, the most of any sample tested, all in the positive direction. PI3K was involved in six significant interactions, with _mGluR5, _NMDAR1, _NMDAR2A, _PSD95, _HOMER1A, and _PI3K. These disruptions in PI3K, which controls AKT/mTOR signaling, is consistent with several reports implicating dysregulated mTOR signaling in the VPA model [32]. The amount of Homer_PSD95, Homer_NMDAR2B, and Homer_NMDAR1 were each increased by almost twofold, reflecting increased NMDAR scaffolding and/or expression. Levels of detected NMDAR1 and NMDAR2B were also increased. These data support prior studies showing increased NMDAR expression in rats following VPA exposure in the cortex [35], although note that a separate study did not find differences in mRNA expression in the cortex or hippocampus [50]. Other notable hits included SynGAP_NMDAR2A and B, SAP97_NMDAR1 and 2B, and Homer_mGLUR5. Comparing these results with VPA cortex, only 2/5 QMI hits in the cortex were shared with the hippocampus, Homer_PSD95 and SAP97_HOMER1A. However, several other interactions that were significant in hippocampus were trending towards significance in cortex; for example, Homer_NMDAR1 and _NMDAR2B were increased by 1.29 and 1.37-fold in the cortex, respectively, but were not significant by ANC criteria (see Additional file 1: Table S1).

Comparisons between models

For the most stringent possible clustering analysis between models, we set all non-ANC-significant measurements to 0 and performed unsupervised clustering using the “complete” method, based on the Euclidian distance matrix of all samples. Because interactions that were significant in a single sample are irrelevant for clustering using this method, we only included the 16 interactions that were significant in two or more samples (Fig. 4). The plot highlights the correlation between certain interactions, such as Homer_PSD95 and Homer_NMDAR1, or SynGAP_PSD95 and SynGAP_SynGAP. However, it is clear from this plot that because there were relatively few interactions that reached ANC significance in multiple models, the clustering is not robust; for example, FragileX and CNTNAP2 mice are shown associated with each other on the basis of a single shared ANC-significant hit, Fyn_Fyn.

Fig. 4

Summary of ANC-significant interactions present ≥ 2 models. Columns are clustered by genotype/tissue type, while rows are clustered by each protein interaction/abundance measure. While this format is useful to give an overview of shared ANC hits, so few hits are shared by multiple models that clustering occurs based on only 1–2 common hits, making the clustering unreliable. Model identifiers in blue represent hippocampal tissue, red cortical tissue. Gray bars indicate the potential confounding factors of age and background strain (see Table 2)

To overcome this limitation, we repeated our cluster analysis with all log₂FC data, reasoning that smaller changes that did not reach the high bar for ANC significance could still be informative for clustering analysis. However, we were concerned about noise contributed by interactions that did not change, but fluctuated randomly around 0, so we first performed principal component analysis (PCA) to focus on factors that contributed the most variation to the dataset. PCA was performed on the mean log₂FC values of each interaction for all genotypes/tissue types using default settings. Plotting the data by principal components 1 and 2, which accounted 30.1% and 12.1% of total variation, respectively (Fig. 5a), revealed clear clustering of tissue types within models; in all cases, the hippocampal and cortical tissue shared similar coordinates in PCA space. Both Shank3 models and Fragile X animals were in close proximity in PCA space, and Ube3a^2xTG were near VPA animals. To mathematically determine the relationships between models in PCA space, we used a hierarchical clustering on principal component (HCPC) analysis using default settings in the FactoMineR package and cutting the HCPC tree at the recommended level to maximize inertia gain (Fig. 5b). HCPC yielded four clusters: CNTNAP animals were an outgroup (group1). Group 2 contained all Shank3 models, and FMR1 animals. TSC2 animals, alone in group 3, were clustered on a branch adjacent to group 4, which contained cortical and hippocampal tissue from both VPA and Ube3a models. We calculated approximately unbiased p values for the clustering based on multiscale bootstrap resampling. The co-clustering of Shank3B hippocampus with Shank3Δex4–9^+/− tissues, and the clustering of VPA tissue with Ube3a cortical tissue reached statistical significance (AU < 0.95); AU values for all other branches are shown in Fig. 5b.

Fig. 5

Clustering models by log₂FC matrices. a Principal component analysis of all genotypes/tissue types. Each group is plotted by its PC1 and PC2 values. Points are colored by HCPC clustering show in B. b HCPC clustering of ASD models. Based on the inertia gained by cutting at each level (inset graph), the HCPC program suggested clustering into four groups as shown. Numbers at the branch points show the approximately unbiased (AU) p value calculated by multiscale bootstrap resampling; clusters with AU greater than 95 are strongly supported by the data. Model identifiers in blue represent hippocampal tissue, red cortical tissue. Gray bars indicate the potential confounding factors of age and background strain (see Table 2)

Shared molecular pathology in cluster 4

We noticed from the structure of the clustering that groups 3 and 4 contained two models with known abnormalities in the AKT/mTOR signaling pathway; in fact, mTOR inhibitors have been reported to rescue behavioral deficits in both models [37, 51]. TSC2^+/− mice are heterozygous for a critical inhibitor of the mTOR complex and show sustained mTOR activation and abnormalities throughout the pathway [37]. VPA animals also show abnormal AKT signaling, with a recent report showing reduced levels of AKT and mTOR, as well as reduced ratios of phospho-to-total AKT and mTOR in VPA exposed rats [52]. Ube3a^2xTG mice clustered closely with VPA mice, but AKT/mTOR has never been implicated in this model. Indeed, mining the factors that differentiated HCPC clusters indicated that PI3K was a significant factor that differentiated group 4, and Ube3a hippocampal issue showed a large, but non-ANC-significant increase in PI3K_PI3K (log₂FC = 0.42, NS). We therefore performed phospho-Western blots on cortical samples from an independent cohort of Ube3a^2xTG animals (Fig. 6). AKT phosphorylation was reduced by 41% at p-Ser473, while no difference was observed at p-Thr308. Total AKT levels were similar. Downstream of AKT, mTOR phosphorylation was also similar, as were levels of p~S6. We confirmed previous reports of altered p-AKT levels in cortical tissue from VPA animals and found that p-AKT and total AKT levels were normal in cortical tissue from all other models examined (Additional file 2: Figure S1). These data confirm the predictions of our clustering that Ube3a^2xTG mice share a core deficit in the PI3K/AKT/mTOR pathway with the VPA mice that share the same branch of the HCPC cluster tree.

Fig. 6

AKT phosphorylation is reduced in Ube3a^2xTG mice. a Representative western blots of synaptosomal fractions from adult mice probed with the indicated antibodies and b quantification. N (WT, Ube3a^2xTG) = 5, 6 for all blots except 11, 12 for p-AKTs473 and panAKT. *p < 0.0001 by two-tailed t test

Discussion

The goal of these experiments was to perform a series of identical protein measurements of brain tissue from multiple mouse models of autism, with the aim of cutting through the immense heterogeneity of the diagnostic entity and identifying some underlying points of convergence. We did not expect every animal model to show an identical set of protein network disturbances. Rather, we hypothesized that a set of interactions might be disrupted in more than one model; perhaps we would be able to identify subtypes of genetic autisms that share distinct sets of disrupted interactions. Indeed, our stringent ANC criteria identified several interactions that were common to multiple genetic models, but clustering by only ANC-significant interactions was not robust. Bioinformatics analysis using PCA and HCPC clustered the models by both genotype and tissue type indicated generally similar changes in cortical and hippocampal tissue from the same models. For such a clustering approach to be broadly useful, it would need to make testable predictions about pathologic mechanisms. Indeed, analysis of the interactions that contributed to clustering suggested that Ube3a^2xTG mice might share a molecular deficit with the other models sharing its branch of the tree, namely disrupted AKT/mTOR signaling. Western blots revealed that AKT signaling was disrupted in Ube3a^2xTG and VPA mice, but not other models, confirming clustering predictions. We have therefore successfully identified a set of protein-protein interactions that are disrupted in multiple animal models of autism, clustered models based on high-dimensional QMI data, and used our clusters to make testable predictions about the molecular pathology of closely clustered models.

Proteins or protein interactions that were ANC-significant in multiple models identified here share striking similarity to a set of interactions that we recently reported to be activity-dependent. In response to 5 min of acute stimulation with glutamate, QMI identified significant changes in 26 protein-protein interactions [32]. Homer, Shank and SynGAP were the most connected nodes, each changing its interactions with several other members of the network. Many of these activity-dependent interactions were also identified here as significantly different in ASD models vs. wildtype controls. For example, Homer1_PSD95 and the abundance of the Ras GTPase SYNGAP were each altered in four sample types, and interactions between Homer1_NMDAR1 and SYNGAP_PSD95 were each altered in three sample types. Of the 26 Glutamate-significant interactions, only 15 were included in the QMI panel presented in this paper; but of those 15 interactions, 9 were “hits” in the ASD models (Additional file 1: Table S1, sheet 2). Notably, the directionality of these changes was variable. For example, Homer_PSD95 levels were increased in VPA cortex, but decreased in Fragile X cortex.

Differences in synaptic activity are a defining and unifying characteristic of animal models of autism—virtually, every report of an autism model includes electrophysiological characterization showing altered synaptic transmission. The directionality of change in synaptic activity is also variable between models; for example Ube3a^2xTG mice show reduced cortical excitability [32] while Fragile X [53, 54] mice show increased. Viewed through this activity-dependent lens, the bidirectionality of our data makes more sense. Glutamate stimulation results in dissociation of Homer-PSD95 complexes [32]; thus, the reduced amount of this interaction seen in the fragile X model could reflect the hyperactive tonic signaling that has been previously reported [53, 54]. Conversely, in the VPA model, a reduction in intrinsic cortical activity has been reported [55, 56], which would be predicted to cause increased levels of this activity-labile interaction. Future studies could manipulate activity in ASD models and measure the resulting QMI profiles to directly test this hypothesis and disentangle activity-dependent from activity-independent processes.

However, activity-dependent interactions were not uniformly altered within models; for example, while Homer1_PSD95 and Homer1_NMDAR1 were reduced by activity and increased in the VPA model, SynGAP_PSD95 was also reduced by activity [32], but unchanged in the VPA model. This could imply an underlying dysregulation in the network response to activity, or a de-coupling of normally correlated molecular processes due to differences in the cell’s ability to compensate for some long-term changes better than others. An analogous network-level dysregulation has been observed in transcriptomic analysis of postmortem autism brain tissue, where individual mRNAs show normal levels of abundance, but the coordinated expression of mRNAs is dysregulated, reflecting disrupted regulatory mechanisms [18]. In the future, it will be informative to design experiments that can de-couple acute, activity-dependent changes from long-term, genotype-dependent changes in PPI networks. The stimulus-dependent dynamics of protein interaction networks encode cellular information, such that different cellular inputs lead to different rearrangements of the interactome, encoding different cellular responses [57]. Understanding how information processing through this synaptic network differs in ASD models could lead to further insights into disease pathogenesis.

To our knowledge, only one other study has attempted to subtype a large number of mouse models of autism [58]. This MRI-based study found great heterogeneity in the relative size of many brain areas in ASD models vs. matched controls, but was able to identify clusters of animal models that shared similar patterns of changes. Three models were analyzed by both the current study and the Ellegood et al. study, FMR1, CNTNAP2, and Shank3B. Both Ellesgood et al. and our study clustered the FMR1 and Shank3B mice as neighbors on the same branch of the dendrogram, suggesting both structural and molecular convergence between the two models. Indeed, prior work has shown that Shank3 mRNA is posttranscriptionally regulated by FMR1 [59], and that FMR1 mice show deficits in mGluR signaling [25, 26] that is mediated by Homer and Shank-containing scaffolds [42]. More generally, genetic studies [60, 61] have implicated several genes or gene regulatory networks related to mGluR signaling in autism. It is plausible that our “cluster 2” may represent a subtype encompassing Shank3 and Fragile X models, previously and independently identified in the Ellesgood study. However, the models may have co-clustered in both studies by chance.

The PI3K/AKT/mTOR pathway has also been implicated in many diverse models of ASD [14, 52, 62], so it was noteworthy that two models with known disruptions to the mTOR pathway appeared together in clusters 3/4. TSC2 is directly involved in regulating mTORC1 downstream of AKT, and autism-linked mutations in TSC2 cause increased mTOR activation and a de-coupling of mTOR from AKT [63]. Prenatal VPA exposure causes reduced mTOR pathway protein expression and phosphorylation [52], which we confirmed here (Additional file 2: Figure S1). After our clustering results suggested a potential mTOR deficit in Ube3a^2xTG animals, we found by phospho-Western blots that Ube3a^2xTG animals showed reduced AKT S473 phosphorylation, but normal levels of T308 phosphorylation and normal phosphorylation of other components of the mTOR pathway. Rapamycin treatment has been shown to rescue behavior in both TSC2 and VPA models [51, 63], and future work could explore if correction of AKT phosphorylation in the Ube3a model might similarly correct behavioral deficits.

Limitations

Several limitations of our study should be noted. The background strain was different among several of the models, the age of the VPA mice was different from the other six models, and mice of both sexes were used. Our analysis approach, in which each mutant animal was normalized to a matched littermate control, was designed to cancel out these effects, as well as assay-dependent batch effects, to identify differences caused by each mutation. However, this experimental design prevents us from making wildtype-to-wildtype comparisons (since batch effects cannot be normalized for mice run on different assay plates), so we are unable to unambiguously demonstrate that our clustering was not driven by uncorrected effects stemming from these differences in background, age, or sex. QMI is a candidate-based approach and shares limitations with all antibody-based assays, including potential antibody cross-reactivity and issues of binding epitope access in native protein complexes. The absence of a detected interaction cannot be interpreted as unambiguously indicating that the interaction does not exist in vivo, since occlusion of binding sites could lead to false negative results. We carefully selected and screened all antibodies used in the QMI panel [32], but antibody caveats are unavoidable. We used NP40 detergent after pilot data showed that it produced higher mean matrix MFIs than TritonX100, Digitonin, or Deoxycholate [32]. However, NP40 does not fully solubilize the core postsynaptic density, where several of our protein targets are enriched (discussed in [32]). Since detergents that solubilize the PSD also disrupt protein interactions, detergent selection is necessarily a trade-off, and further studies could more thoroughly quantify differences in synaptic QMI networks due to different detergent conditions. Finally, many of our interactions vary with neuronal activity or by brain region. Small variations in microdissection (e.g., inclusion of small amounts of striatal tissue in cortical samples) or euthanasia protocols (i.e., animal sleep/wake state prior to euthanasia) could have large effects on protein detection (e.g., Shank3, which is highly expressed in striatum, or PSD95_SynGAP, which is activity-dependent). While we were careful to perform our dissections as consistently as possible, at a similar time of day and using metal brain molds to ensure consistent slicing, thinner vibratome slicing followed by a period of controlled slice recovery in ACSF, as for electrophysiology, may be a more optimal experimental strategy to ensure normalization of both activity and location.

Conclusions

In conclusion, we performed a series of identical QMI experiments to measure differences in the abundance of, and binary interactions among, 16 synaptic proteins in 7 mouse models of autism. Employing a mutant-littermate control design, we found a unique combination of disrupted protein interactions in each model and tissue type measured. Many of the disrupted interactions were identified as activity-dependent interactions in a separate study, highlighting the complex relationships between ASD risk genes and activity-dependent homeostatic processes [21]. PCA and cluster analysis of models revealed two identifiable sub-groups, with VPA and TSC2 mice comprising a hypothetical “mTOR” cluster, and Shank3 and FragileX mice comprising a second cluster; the latter co-clustering was consistent with a prior MRI study [58]. The inclusion of Ube3a^2xTG mice in the mTOR cluster led to our identification of AKT phosphorylation deficits in this model. Our data highlight the heterogeneity of ASD models, while offering hope that high-dimensional measures of biologically relevant molecular processes may allow differentiation of subtypes of ASD amenable to common treatment strategies. Future work to expand the number of ASD models analyzed and to perform similar QMI experiments in human iPS-derived neurons could offer further insights into the complex relationships among autism risk factors.