Multi-influential genetic interactions alter behaviour and cognition through six main biological cascades in Down syndrome mouse models.

Down syndrome (DS) is the most common genetic form of intellectual disability caused by the presence of an additional copy of human chromosome 21 (Hsa21). To provide novel insights into genotype-phenotype correlations, we used standardized behavioural tests, magnetic resonance imaging and hippocampal gene expression to screen several DS mouse models for the mouse chromosome 16 region homologous to Hsa21. First, we unravelled several genetic interactions between different regions of chromosome 16 and how they contribute significantly to altering the outcome of the phenotypes in brain cognition, function and structure. Then, in-depth analysis of misregulated expressed genes involved in synaptic dysfunction highlighted six biological cascades centred around DYRK1A, GSK3β, NPY, SNARE, RHOA and NPAS4. Finally, we provide a novel vision of the existing altered gene-gene crosstalk and molecular mechanisms targeting specific hubs in DS models that should become central to better understanding of DS and improving the development of therapies.

Introduction

Down syndrome (DS) is the most common genetic form of intellectual disability and was first described as a disease by John Langdon Down in 1866. A century later, genetic studies demonstrated that DS is caused by the trisomy of human chromosome 21 (Hsa21) (1). People with DS have a wide range of phenotypic and physiological features with phenotypic variability, but they always present several disabling features like intellectual disability or Alzheimer’s disease (2). The leading cause of DS is the non-disjunction of chromosome 21 (3). However, in rare cases, people with partial Hsa21 duplications have been observed with a smaller spectrum of DS features. Studying these rare conditions increased our understanding of the genotype–phenotype correlations in DS (4–10): there is no single trisomic region responsible for all DS features, rather there are several susceptibility regions when presented in three copies that contributes to DS features in people with partial duplication of Hsa21. Consequently, this can induce a wide variety of features (5,7). Nevertheless, several individuals displayed complex rearrangements such as contiguous or non-contiguous deletions or duplications, copy number variations of other regions or the duplication of genes located in the short arm of Hsa21. These rearrangements can potentially impact the phenotypic outcome of the Hsa21 duplication and add noise to the genetic dissection of clinical manifestations of human trisomy 21.

To circumvent the difficulties of studying DS in humans, several efforts have been made to generate DS mouse models (11). Indeed, there are three independent mouse regions homologous to Hsa21, carrying altogether 158 protein-coding homologous genes of the 218 protein-coding genes identified on Hsa21 (12). The largest region is found on mouse chromosome 16 (Mus musculus 16, noted as Mmu16) in a fragment of 22.42 Mb with 119 orthologous genes between Lipi and Zbtb21 (13). The most telomeric part is split into two parts. The first part is found on mouse chromosome 17 (noted as Mmu17) with 19 homologous genes in the 1.1 Mb interval between Umodl1 and Hsf2bp. Then, the second part is on mouse chromosome 10 (Mmu10) with 37 genes included in the 2.3 Mb Cstb-Prmt2 genetic interval (11,12,14). Several DS mouse models have been generated over the years, most of them carried trisomy of the largest genetic region located on Mmu16 (11,14). The Ts(1716)65Dn (noted Ts65Dn) model is the most widely used DS animal model and is quite unique, having a supplementary mini-chromosome obtained by X-ray irradiation of the male germline and containing the centromeric region of Mmu17, with genes from Psid-ps2 to Pde10a, and the 13.5 Mb telomeric fragment of Mmu16 containing genes between Mrpl39 and Zbtb21 (15–18). Several models were made by chromosomal engineering (11) and carried the segmental duplication of Mmu16. The Dp(16Lipi-Zbtb21)1Yey (noted Dp1Yey) corresponds to the duplication of the entire Mmu16 region syntenic to Hsa21 (19). The Dp(16Cbr1-Fam3b)1Rhr (noted Dp1Rhr) model carries a duplication from Cbr1 to Fam3b and demonstrates the contribution of the DS critical region (DCR) (20–23). All the DS mouse models displayed defects in behaviour and cognition which had been investigated in different laboratories with different protocols and environmental conditions, making inter-model comparison difficult (11).

To improve our knowledge of DS, we analysed seven DS mouse models that carry either a large segmental duplication, like Dp1Yey, or a transgenic line overexpressing Dyrk1a, the main driver gene of the phenotype in mouse DS models, found on Mmu16 (24–28), and specific combinations of models (see Fig. 1A). We used a unique and in-depth behaviour, morphological and transcriptomics pipeline in adult animals to dissect the contribution of the genes located on Mmu16 to DS mouse features. The behaviour pipeline relied on assessing specific hippocampal-dependent brain functions identified as altered in people with DS (29,30). Thus, we performed standardized Y-maze, Open field (OF), Novel Object Recognition (NOR), Morris Water Maze (MWM) and contextual and cue Fear Conditioning (FC) tests. All the procedures were carried out in similar environmental conditions to reduce any impact of variation (31,32). In addition, variations in specific brain regions have been observed in people with DS and mouse models (24,33–35). Neuroanatomical changes affect the whole brain volume or specific regions like the frontal region of the limbic lobe and the hippocampus in people with DS. Thus, we performed an in-depth morphological investigation of the brain by magnetic resonance imaging (MRI). Finally, whole gene expression was performed on hippocampi isolated from the six models to decipher the genes, functional pathways and biological cascades affected in the different DS mouse models.

Figure 1

Down syndrome mouse models analysed in the study (A) and Standardized behavioural profiling (B–E). (A) In the upper part of the plot the human chromosome 21 is represented, in yellow we highlighted the Hsa21 syntenic region found in mouse from Lipi to Zbtb21 (known as Zfp295 previously). The eight models analysed on this study Dp1Yey, Dp3Yah, Ts65Dn, Dp5/Dp1, Dp5Yah, Dp1Rhr, Tg(Dyrk1a), Dp5yah crossed with Tg(Dyrk1a) (noted as Dp5-Tg)) trisomic chomosomal regions were draw in comparison with the Hsa21 region. (B) Y-maze spontaneous alternation. Arm visited (A upper panel) and alternation rate (A lower panel) are presented as box plots with the median and quartiles (upper and lower 90% confidence interval are indicated by a grey box). Only the Tg(Dyrk1a) mice showed hyperactivity in this test with increased arms entries compared to the wild type (p = 0,017). Alternation rate in Dp1Yey (p = 0,002), Ts65Dn (p < 0,001), Dp1Rhr (P = 0,012), Dp5/Dp1 (p = 0,018) and Tg(Dyrk1a) (P = 0,010) mice was significantly lower than respective wild-type mice (Dp1Yey n = 10 wt and 10 Tg; Ts65Dn n = 14 wt and 14 Tg, Dp5/Dp1 n = 17 wt, 16 Dp5Yah, 15 Dp1Rhr and 17 Dp5/Dp1; Tg(Dyrk1a) n = 11 wt and 14 Tg). (C) Exploratory activity in a novel environment. Distance travelled (B upper panel) and % of distance travelled in peripheral zone recorded in the Open field arena (B lower panel). The total distance travelled was significantly higher in Ts65Dn (p = 0,022), Tg(Dyrk1a) (p = 0,008) and Dp5/Tg(Dyrk1a) (p>0,001). Moreover, the % of distance in the peripheral zone was increased in Ts65Dn (p > 0,001) mice comapred to wild type mice (Dp1Yey n = 10 wt and 10 Tg; Dp3Yah n = 15wt and 15 Tg; Ts65Dn n = 14 wt and 14 Tg, Dp5/Dp1 n = 17 wt, 16 Dp5Yah, 15 Dp1Rhr and 17 Dp5/Dp1; Tg(Dyrk1a) n = 11 wt and 14 Tg). (D) Novel Object Recognition with 24 hour (D upper panel) or 1 hour retention time (D lower panel). The results are presented as % of sniffing time (as box plots with the median and quartiles) for the novel object (NO). For 24 hours time laps, one sample t-test vs 50% (hazard) showed that Dp1Yey (p = 0,837), Dp3Yah (P = 0,173), Ts65Dn (p = 0;432), Dp1Rhr (p = 0,492), Tg(Dyrk1a) (p = 0,144) and Dp5/Tg(Dyrk1a) (P = 0,488) failed to recognize the new object. The Dp5Yah genomic fragment restored the capacity of the Dp1Rhr in the Dp5/Dp1 mice (p = 0,0157; Dp1Yey n = 10 wt and 10 Tg; Dp3Yah n = 15 wt and 15 Tg; Ts65Dn n = 14 wt and 14 Tg, Dp5/Dp1 n = 17 wt, 16 Dp5Yah, 15 Dp1Rhr and 17 Dp5/Dp1; Tg(Dyrk1a) n = 11 wt and 14 Tg). For 1 hour retention time, all the mice were able to discriminate the NO except for the Tg(Dyrk1a) (p = 0,011 preference for FO; Dp1Yey n = 10 wt and 10 Tg; Dp3Yah n = 15 wt and 15 Tg; Ts65Dn n = 5 wt and 5 Tg; Tg(Dyrk1a) n = 11 wt and 12 Tg). (E) Probe test session in Morris Water Maze. The results are presented as % of time in the target quadrant. All the mice have spent more time in the target quadrant versus non target excepted for the Ts65Dn mice (p = 0,398) (*p < 0.05, **p < 0.01, ***p < 0.001; Dp1Yey n = 9 wt and 10 Tg; Dp3Yah n = 15 wt and 15 Tg; Ts65Dn n = 10 wt and 11 Tg, Dp5/Dp1 n = 13 wt, 13 Dp5Yah, 13 Dp1Rhr and 13 Dp5/Dp1; Tg(Dyrk1a) n = 16 wt and 15 Tg).

Results

Dissecting the contribution of Mmu16 sub-regions to DS-related behavioural phenotypes in mouse models

We wanted to dissect the contribution of sub-regions located in the telomeric part of M. musculus chromosome 16 (Mmu16), homologous to Hsa21 (36), to DS-related cognitive phenotypes. First, we selected four DS mouse models: Ts65Dn, the most commonly used DS model (15), and three additional models that carry segmental duplications of well-defined sub-regions located on Mmu16, Dp1Yey (19), Dp3Yah (37) and Dp1Rhr (20). In addition, we engineered a new sub-region, Dp5Yah, corresponding to three functional copies of the genes included in the genetic interval between Cyyr1 to Clic6. This model was crossed with the Dp1Rhr sub-region to generate Dp5Yah/Dp1Rhr (noted Dp5/Dp1) compound transheterozygote carrying a trisomic gene content similar to Ts65Dn for the genes located on Mmu16. We also included a model carrying an additional copy of Dyrk1a, one of the driver genes for DS-related phenotypes (24,25), and Tg(Dyrk1a) combined with the Dp5Yah model (noted Dp5-Tg) (see Fig. 1A). We used standardized behavioural tests to study several aspects of learning and memory in mice, including the Y-maze (working memory), OF (exploration memory), NOR (recognition memory), MWM (spatial memory) and FC (associative memory) tests. For all the lines, independent cohorts of control and trisomic mouse littermates went through the pipeline at similar ages, after which the resulting data were processed using standard statistical analyses (see Supplementary Information for details). First, we assessed the potential existence of a background effect in the distribution of the measurements taken in the different tests. In our condition, the Q–Q plots with cumulative frequency were linear (see Supplementary Material, Table S1, Supplementary Material, Fig. S1) and thus, no notable difference between B6J or hybrid B6JC3B wild type controls was observed.

Mouse activity and working memory were evaluated in the Y-maze (Fig. 1B). The number of arm entries in the Y-maze showed that only the Tg(Dyrk1a) mutant line was hyperactive in this test, while a deficit in spontaneous alternation was found in Dp1Yey (38), Ts65Dn (39), Dp1Rhr and Tg(Dyrk1a) (40). Dp5/Dp1 also showed a clear deficit in the percentage of spontaneous alternation in comparison with littermate controls, while Dp5Yah trisomic animals showed normal performance. Thus, the overexpression of Dyrk1a alone was sufficient to mimic the Y-maze spontaneous alternation found in three overlapping DS models but we cannot rule out the possibility that another region is involved, as suggested by the work of Chang et al. (41).

Patterns of exploratory activity and anxiety were assessed in the OF (Fig. 1C). Ts65Dn, Tg(Dyrk1a) and Dp5-Tg presented hyperactivity with an increased travelled distance compared with wild type littermates, and support the results obtained in the Y-maze for the Tg(Dyrk1a) line. Thus, the hyperactivity phenotype is linked either to the Ts65Dn mouse model or only to the increase in Dyrk1a dosage.

The spatial reference memory was tested in the standard MWM task, in which mice have to escape from a circular pool of opaque water by localizing a hidden platform at a single fixed location using distal spatial cues. We analysed the velocity, the distance travelled by the mice to reach the platform and thigmotaxis over-training (Supplementary Material, Fig. S2). The velocity of Dp1Yey and Dp5Yah was slightly lower than that of the wild type mice (See Supplementary Material, Table S1 and Supplementary Material, Fig. S2B). As described previously (22,42–48), Ts65Dn mice displayed a longer distance travelled to find the platform during all the sessions, compared with the wild type group (Supplementary Material, Fig. S2A). Although Tg(Dyrk1a) mice were able to locate the platform, they also showed delayed acquisition compared with the control mice. Surprisingly, the Dp1Yey, Dp1Rhr, Dp5Yah and Dp3Yah mice completed this test without any difference with the wild type group. Interestingly, the Ts65Dn model, and to a lesser extent the Tg(Dyrk1a) model, presented marked thigmotaxic behaviour, spending a higher percentage of time in the peripheral zone in comparison with controls (Supplementary Material, Fig. S2C). The retention of place location was evaluated during a single probe trial (PT) with no platform available, 24 h after the last training session (Fig. 1E). All the mouse strains except Ts65Dn remembered where the platform was located after the learning sessions. Finally, to check the visual ability of the mice, we performed a visual training version of the MWM during which the platform position was indicated by a flag. All the mice were able to find the visible platform without any significant difference with controls except for Tg(Dyrk1a), which presented a small delay in session 2 (Supplementary Material, Fig. S2A).

We then evaluated non-spatial recognition memory using the NOR paradigm with a retention time of 24 h. The percentage of sniffing time for the novel object was analysed and compared with 50% (hazard). This analysis showed that Dp1Yey, Dp3Yah, Ts65Dn, Dp1Rhr and Tg(Dyrk1a) were not able to discriminate between familiar and novel objects, unlike Dp5Yah and, more surprisingly, Dp5/Dp1 (Fig. 1D). To further characterize the effect/lack of effect of Dp5Yah mutation on NOR, the Dp5Yah mouse line was crossed with the Tg(Dyrk1a) one and compared with new sets of wild type, Dp5Yah and Tg(Dyrk1a) mice. Interestingly, we found that Dp5Yah/Tg(Dyrk1a) and, as expected, Tg(Dyrk1a), displayed altered novel object discrimination while Dp5Yah spent more time exploring the novel object than the familiar one. We also assessed the short-term memory of Dp1Yey, Dp3Yah, Ts65Dn and Tg(Dyrk1a), by performing the NOR with a 1 h delay between acquisition and retention; only the Tg(Dyrk1a) mice showed a deficit. Altogether, these data suggested that several loci are involved in the cognitive phenotype: one located in the Dp3Yah region and Dyrk1a; and presumably two interacting modifier loci: one located in the Dp5Yah region and another in the Dp1Rhr region.

All the trisomic lines were also tested for associative memory in the Pavlovian FC test. All the groups showed a higher percentage of freezing during the 2 min post-shock compared with the habituation session, indicating that all the groups developed fear responses during the training session (Supplementary Material, Fig. S3). When the animals were re-exposed 24 h later to the same context, the level of freezing in all the groups was increased compared with the habituation session (PRE2 and PRE4). However, the freezing time for Ts65Dn mice was lower compared with the respective control littermates. When we assessed cued FC in a new context, all the mice presented longer immobility time with a marked difference between pre-cue and cue periods (Supplementary Material, Fig. S3). In addition, Dp1Yey and, to a lesser extent, Ts65Dn showed lower freezing during the presentation of the cues in comparison with wild type counterparts. As such, altered emotional associative memory was found in Ts65Dn and Dp1Yey but not in other DS models. Thus, depending on the behavioural traits observed, the overlapping DS models pointed to different regions with causative and modifier loci involved spread along the Mmu16 region homologous to Hsa21.

Dissecting the contribution of Mmu16 sub-regions to the DS-related brain morphological phenotypes in mouse models

DS models have been reported to show brain morphological alterations of specific regions (24). Thus, we wondered if we could detect changes in brain morphology using MRI on these different partially overlapping trisomic mice models. Data were first analysed using a volume approach and a brain region atlas. We confirmed that the brain of Tg(Dyrk1a) mice was larger (P < 0.001) (24) and the brain of Dp1Yey mice was smaller than the respective wild type. Then, we analysed different brain regions/structures taking into consideration the whole brain volume. Even with this correction, the Tg(Dyrk1a) mice were the most affected in terms of brain structure volume whereas, on the contrary, the Dp1Rhr mice did not show any significant variation compared with the wt mice. Several regions, such as the basal forebrain septum, central grey matter, the rest of the midbrain and superior colliculi were significantly larger in the Tg(Dyrk1a), Ts65Dn and Dp1Yey DS models. Moreover, the cerebellum, hypothalamus, inferior colliculi and caudate putamen were significantly different in size for the Dp1Yey and Tg(Dyrk1a) carriers compared with controls (Supplementary Material, Table S2 and Supplementary Material, Fig. S4), and a few additional areas were altered specifically in certain models (amygdala, globus pallidus, hippocampus, neocortex and thalamus for Tg(Dyrk1a); external capsule, fimbria and ventricles for Dp1Yey). Altogether, this brain morphometric analysis showed greater similarity between the Dp1Yey and Dyrk1a overexpression transgenic models, with an intermediate overlap with the Ts65Dn mouse model.

Dissecting the contribution of Mmu16 sub-regions to the DS-related transcriptome profiles in mouse models

Various studies have shown the consequences of trisomy on gene expression (49–58). Here, we took the opportunity to dissect the alteration of gene expression and functional pathways in various DS trisomic models carrying different duplications of Mmu16. We analysed Ts65Dn, Dp1Yey, Dp3Yah, Dp5/Dp1, Dp5Yah, Dp1Rhr, and we included the trisomic model for Dyrk1a alone, Tg(Dyrk1a). Considering the hippocampal formation as a hub structure involved in learning and memory, we performed gene-expression analysis in the adult hippocampus by comparing the DS models with their own littermate controls using a unique pipeline for all the models. For each DS model, we defined the expressed genes (noted as EGs) as the genes whose expression level was detected; the differentially expressed genes (noted as DEGs) as the genes whose expression level was found to be significantly altered in the trisomic model compared with the controls littermates; and then the trisomic expressed genes (TEGs) as the DEGs that are included within the duplicated chromosomal regions for each model (Supplementary Material, Table S3 and Table 1).

Table 1

Differential expression analysis results of the seven models analysed

Mouse lines	Dp1yey	Dp3Yah	Ts65Dn (Mmu16)	Dp5/Dp1	Dp5Yah	Dp1Rhr	Tg(Dyrk1a)	Ts65Dn (Mmu17)
Nb of probes detected by the array	35556
Nb of annotated probes detected (without control probes)	27359
Number of trisomic expressed genes (TEGs)	155	21	130	127	87	40	1	40
Number of differential expressed genes (DEGs)	711	826	1074	922	736	1306	850
Number of differential trisomic expressed genes TEGs	66	13	64	54	39	18	1	28
% of TEGsa	43%	62%	49%	43%	45%	45%	100%	70%
% of compensated trisomic genes detected	57%	38%	51%	57%	55%	55%	0%
Number of GAGE KEGG and GOs (CC, BP, MF) terms disregulated in the trisomic model FDR < 0.1	244	67	12	111	318	225	231
Number of GAGE KEGG and GOs (CC, BP, MF) terms upregulated in the trisomic model FDR < 0.1	207	60	3	33	4	132	222
Number of GAGE KEGG and GOs (CC, BP, MF) terms downregulated in the trisomic model FDR < 0.1	37	7	9	78	314	93	9
Number of GAGE KEGG and GOs (CC, BP, MF) terms in the trisomic model unique to each mice line FDR < 0.1	80	1	2	30	195	107	64

GO are Go functional terms involved in cellular compartment (CC), molecular function (MF) and biological processes (BP).

Although most of the genes in three copies were overexpressed in the relevant mouse model-derived hippocampi with a ratio around 1.5 (Supplementary Material, Fig. S5), from 38 to 57% of the trisomic genes showed a dosage compensation (Table 1 and Supplementary Material, Table S3). While this compensation was expected, we noticed that most of the compensated genes behaved similarly in the different trisomic contexts, although the experiments were performed independently. As such, the genes from Cldn17 to Krtap11–1, including the keratin cluster were not overexpressed when trisomic in any model (Dp1Yey, Ts65Dn, Dp5/Dp1 or Dp5Yah). This could have been due to the fact that this region seems to be under strong regulatory constraints, as on the borders two REST sites and a LaminB1 peak encompassing this region were found (UCSC browser), while Btg3 (BTG Anti-Proliferation Factor 3) and C21orf91 (Chromosome 21 open reading frame 91, also known as D16Ertd472e), and Mrpl39 (Mitochondrial Ribosomal Protein L39), Jam2 (Junctional Adhesion Molecule 2), Atp5J (ATP synthase peripheral stalk subunit F6), Gabpa (GA binding protein transcription factor subunit alpha) and App (amyloid beta precursor protein) were overexpressed in various DS models. We also found that the genes located on the trisomic segment not homologous to Hsa21 on Mmu17 in Ts65Dn hippocampi, were overexpressed (Supplementary Material, Fig. S6). This was also observed in the Ts65Dn heart in a previous study (16). Looking in detail at the homologous Hsa21 region in Mmu17, we saw two main genetic effects due to the overdosage of the Mmu16 region homologous to Hsa21. Noticeably, Cbs, coding for the Cystathionine beta-synthase, another driver gene for DS cognitive phenotypes (59), was found to be downregulated in all the models, except Dp3Yah and Tg(Dyrk1a), suggesting direct control of Cbs transcription by at least two loci, one located in Dp5Yah and another one, not due to Dyrk1a overdosage, in the Dp1Rhr trisomic region. Similarly, under-expression of the glucagon like peptide 1 receptor (Glp1r) was observed in Dp1Yey, Ts65Dn, Dp5/Dp1 and Dp5Yah. On the contrary, this gene was overexpressed in Dp1Rhr and not affected in Tg(Dyrk1a).

Here too, Dyrk1a dosage was not involved, but at least two loci controlling Glp1r expression with opposite and epistatic effects should be found respectively in the Dp5Yah and Dp1Rhr genetic intervals. Thus, a complex genetic interaction took place between different loci, controlling subsequent gene expression.

The analysis of DEGs in each model by principal of component analysis (PCA), t-SNE (Fig. 2B) or OPLS techniques (see Supplementary Information) highlighted the capabilities to separate trisomic individuals from wild type littermates (Fig. 2A). Genome-wide misregulation was found independently of the model, as DEGs were spread in all the chromosomes (Supplementary Material, Tables S3, S4 and Supplementary Material, Fig. S7), as shown previously (55), although with a stronger impact of the Dp1Rhr on the total number of DEGs detected.

Figure 2

The differential expression analysis discriminates trisomic from disomic hippocampi and identifies common dysregulared genes and pathways. (A) 3D-PCA on the DEGs for each sample allows to separate trisomic (Dp) from disomic (wt) adult hippocampi. (B) Left column: two dimensional Principal Component Analysis (2D-PCA) on the 5599 transcripts of the 4328 DEGs over all the samples identified using fcros 0.025<α<0.975. Right column: 2D-PCA on the 75 trisomic genes with a measured expression in all the models. (C) Venn Diagrams showing the overlap in gene expression between the different mouse lines represented in different colours. (D) Heatmap representation of the number and regulation sense of the pathways shared at least by two mice lines identified using the genome expression for each mice line by GAGE R package and filtered by q-value cut off < 0.1. Pathways are grouped in main meta-pathways showed on the ordinate: Synaptic related = all the synaptic related pathways excluding myelin sheath and SNARE complex formation; Transcription & epigenomics regulation; Enzymes activity; Ribosome related; Mitochondria related; Cell Structure & organelle related; Phospho-kinase related. The color key breaks represents the number of pathways within the meta-pathways 60,40,20,5,0.5. The minus or pink color represents down regulated pathways, the white color represents no pathway found in the meta-pathways and the purple or positive numbers stands for up regulated pathways.

The most overexpressed genes in terms of log2 fold change (log2FC) of expression and significance in various genetic conditions were visualized using Volcano plots (Supplementary Material, Fig. S8 and Supplementary Material, Table S4). For example, the listerin E3 ubiquitin protein ligase 1 (Ltn1) gene, coding for a major component of ribosome quality control and causing neurodegeneration in mice (60), was found overexpressed in Dp1Yey, Ts65dn, Dp5/Dp1, and Dp5Yah or Ifnar2, coding for the Interferon Alpha and Beta Receptor Subunit 2, was overexpressed as expected in models that carried three copies of this gene (Dp1Yey, Ts65Dn, Dp5/Dp1 and Dp5Yah). Instead, a more controlled gene like the neuronal acetylcholine receptor subunit alpha-3 (Chrna3) was found upregulated only in Dp1Rhr and Dp1/Dp5, certainly due to the overexpression of one gene from the Cbr1-Fam3b region but not Dyrk1a. Nevertheless, when we performed the intersection between the list of DEGs from the different models, we found only a few genes in common (Fig. 2C and Supplementary Material, Table S4).

We decided to combine the analysis of all the lines together using PCA and t-SNE and revealed a strong clustering of models that shared at least a partially duplicated region (Fig. 2B). The t-SNE analysis, based on all the 4328 DEGs detected in each mouse model added together, showed different contributions of the various DS models to the transcriptome variation (Fig. 2B, left) with two distinct groups: one encompassing four overlapping trisomies: Ts65Dn, Dp5/Dp1, Dp5, Dp1Rh and three isolated models: Dp1Yey, Dp3Yah and Tg(Dyrk1a) that were closer together, although Dp3Yah was clearly farthest from the other two. Similar distinct groups were seen when analysing the TEGs (Fig. 2B, right) and overall, the trisomic and the wild type individuals in each mouse line were nicely separated. As expected, the expression level of the TEGs and the DEGs in the different trisomic conditions were strongly correlated (Supplementary Material, Fig. S9). Interestingly, the 4328 DEGs showed a level of misregulation strongly correlated between Dp1Yey and Dp3Yah (33%), Dp5/Dp1 (50%), Dp1Rhr (40%) and Tg(Dyrk1a) (42%). Of the 75 genes detected and located on Mmu16 region homologous to Hsa21, the correlation of fold change was around 28% in the Dp1Yey and Tg(Dyrk1a) partial DS models. Thus, the correlation in gene deregulation showed that Dyrk1a overdosage is a key driver of transcriptome deregulation in the Dp1Yey and Dp1Rhr models.

Unexpectedly, the correlation of DEGs mis-expression level was lower between Ts65Dn and Dp1Yey (29%) or Dp1/Dp5 (28%). On the contrary, a large number of TEGs were misregulated in the same way between Ts65Dn and Dp1Yey (49%) and Dp1/Dp5 (52%; Supplementary Material, Fig. S9), suggesting that the other region found in three copies in the Ts65Dn over Mmu17 must affect the general DEG landscape.

Using qRT-PCR, we confirmed the mRNA overexpression of first, Dyrk1a and Sod1 genes in the DS models where they were trisomic; second, of Synaptojanin2 (Synj2) and T lymphoma invasion and metastasis inducing gene 2 (Tiam2) located on the Mmu17 centromeric region in Ts65Dn, and third, of Cholinergic Receptor Nicotinic Alpha 1 Subunit (Chrna1), a gene misregulated in the Dp1Yey, Dp5/Dp1, Dp1Rhr and Ts65Dn models (Supplementary Material, Fig. S15), and Cbs downregulated in all the models except Tg(Dyrk1a) and Dp3Yah. We also confirmed alterations of the expression of immediate early response genes Arc, FosB, Fos and Npas4 that are important for cognition.

Differential functional analysis unravels a few common altered pathways in DS models

To go further, we performed a differential functional analysis and found 12–318 misregulated pathways in the DS models (Table 1 and Supplementary Material, Table S5). Interestingly, the regulation of pathways is trisomic region-dependent, as the Dp5Yah (99%) region produced an overall downregulation whereas the Dp3Yah (89.5%) and Dp1Rhr (56%) trisomy together with the full Hsa21 syntenic model Dp1Yey (84.8%) induced preferentially an upregulation.

To facilitate understanding, we clustered the broad functional dysregulation into ten major functionality groups or meta-pathways. We found ribosomal components and mitochondrial process pathways altered in all the models, with many genes shared between models (Supplementary Material, Fig. S11). Cell structure and organelles, transcription and epigenetic regulation, synaptic meta-pathways and interferon pathway were more affected in some models than in the others (Supplementary Material, Figs S10 and S12). As such, we observed strong and connected effects in the control of transcription and epigenetic regulation, enzyme activity and cell structure, and cellular organelles involved in membrane and protein processing (endoplasmic reticulum, Golgi body, lysosome, peroxisome, etc.; Fig. 2D) in the Dp1Yey, Dp5/Dp1, Dp1Rhr and Tg(Dyrk1a) models, whereas the myelinization and 10 SNARE components, such as the Synaptosome Associated Protein genes 25 and 23 (Snap25 and Snap23), were specifically dysregulated in the Dp1Yey Dp5Yah and Tg(Dyrk1a) models.

Interestingly, we saw many shared genes between these pathways and the models, giving rise to high pathway connectivity between models (see Materials and Methods). Considering the DEGs involved in brain synaptic pathways, with the DS synaptic MinPPINet (Fig. 3A), we analysed the DS network topography and betweenness connectivity and found hubs and genes more central for network information flow. As expected from a protein–protein interaction (PPI) biological network, the likelihood of observing such connectivity in the DS network was more than one can expect by chance (P-value < 2e-16) and it showed a small world effect and scale-free topology. Using a network decomposition approach (see Supplementary Information), we highlighted six major sub-networks or biological cascades that strongly centralized six different proteins: DYRK1A, GSK3β, NPY, RHOA, SNARE and NPAS4 proteins (Fig. 3B and C). As a summary, Dyrk1a was an upregulated in Dp1Yey, Ts65Dn, Dp5Dp1, Dp1Rhr and Tg(Dyrk1a) while Npas4 was downregulated in Ts65Dn, Dp1Rhr and Tg(Dyrk1a), and Npy was upregulated in Dp1Yey, Dp5Dp1, Dp5Yah, Dp1Rhr and Tg(Dyrk1a), and downregulated in Dp3Yah. Ten genes from the SNARE complex were dysregulated in some DS partial models; from these we validated the disregulated expression of Snap25 and Snap23 by qRT-PCR (Fig. 4F). Gsk3β and Rhoa were not DEGs but these two were hubs interacting with many DEGs in the network.

Figure 3

Protein-protein interaction networks involving DEGs linked to the synaptic function. (A) STRING04 MinPPINet of genes involved in synaptic function visualized using the edge weighed spring embedded layout by betwenness index in Cytoscape. The network was built by querying STRING and selecting the PPIs with a medium confidence score (CS = 0.4) coming from all sources of evidence. The shapes of the nodes represent the following information: Shapes: i) Pallid pink ellipses: represent connecting proteins added to assure the full connectivity of the network; ii) pink octagons, represent Hsa21 syntenic genes in mouse not identified as contributing to the meta-pathway dysregulation by GAGE; iii) green inner coloured ellipses, genes identified by GAGE after q-val <0.1 cut off to be contributing even slightly, to any pathway of those found dysregulated inside the meta-pathway. If the size is similar to the octagons, they are also HSA21 syntenic genes in mouse. Additionally, the border colour represents the mouse model multi group where those genes are found altered in; iv) diamonds, genes identified by GAGE after q-val <0.1 cut off and also by FCROS as DEGs. (B) Network Structure Decomposition of the STRING04 MinPPINet. Highlighting in different colors the interactions of GSK3β, NPY, SNARE proteins, DYRK1A and RHOA respectively. In the case of NPAS4, the interactions coloured correspond up to the first level interactions. (C) The six RegPPINets were extracted from the selection of each fo the following proteins and their 2nd interactors from STRING04 MinPPINet: RHOA, DYRK1A, GSK3β, NPY, SNARE proteins and NPAS4. Then, those were further annotated with regulatory information using REACTOME (See Supplementary information). The shapes of the nodes represent the following information: Shapes: i) Pallid pink ellipses: represent connecting proteins added to assure the full connectivity of the network; ii) pink octagons, represent HSA21 syntenic genes in mouse not identified as contributing to the meta-pathway dysregulation by GAGE; iii) green inner coloured ellipses, genes identified by GAGE after q-val <0.1 cut off to be contributing even slightly, to any pathway of those found dysregulated inside the meta-pathway. If the size is similar to the octagons, they are also Hsa21 syntenic genes in mouse. Additionally, the border colour represents the mouse model multi group where those genes are found altered in; iv) diamonds, genes identified by GAGE after q-val <0.1 cut off and also by FCROS as DEGs. The edges colored represent the type of interaction annotated by following the PathPPI classification (61), and ReactomeFIViz annotations as follows i) The GErel edges indicating expression were colored in blue and repression in yellow. ii) PPrel edges indicating activation were coloured in green, inhibition in red. iii) Interactions between proteins known to be part of complexes in violet. iv) Predicted interactions were represented in grey including the PPI interactions identified by STRING DB (62) after merging both networks.

Figure 4

Evaluating RHOA and NPAS4 pathways in DS models. (A) RHOA pathway was altered in the Dp1Yey. Western blot analysis was revealed no changes in the expression of RHOA but a significant decrease of the phosphorylation of the Myosin light chain (P-MLC) in the Dp1Yey hippocampi compared to control (n =  5 per genotype). (B) Mice were housed in total darkness for 14 days and then were subsequently exposed to light for 0, 1, 3 or 7.5 h. Relative expression levels were determined, and fold change were calculated for each condition. Genotypes differences in fold change were assessed by Student t-test. (C) Only the fold change for early response genes Npas4 was up-regulated in Tg(Dyrk1a) mice compared to wt at 1 hours of light induction. (D) The late responses genes specific to inhibitory neurons Frmdp3, Slc25a36 and Igf1 were up-regulated after 7.5 hours of light induction. (E) The fold change of late responses genes specific to excitatory neurons Bdnf and Nrn1 were unchanged. (F) The fold change of late response genes shared by excitatory and inhibitory neurons Gpr3 and Nptx2 were downregulated after 3 and/or 7.5 hours of light induction. (G) Dyrk1a and RhoA showed a similar fold change along the different condition whereas Snap25 and Snap23 presented an increased enrichment for the 7.5 hours condition. Data are presented as box plots with the median and quartiles.

Overall, the network analysis of the DS synaptic RegPPINets, containing the MinPPINET and the regulatory information of the interactions, showed that DYRK1A controls 42.3% of the network nodes and 69.4% of the network seeds via second-level interactors. Hence, DYRK1A could control the DS synaptic network via PPI and regulatory interactions. Furthermore, the biological cascades centred on GSK3β, DYRK1A and RHOA are highly interconnected (Supplementary Material, Figs S13 and S14) and in fact several interactors of RHOA are connected and could somehow modulate a higher percentage (75 and 68.5%) of the nodes of the network and synaptic seeds (Supplementary Material, Table S6).

Validating the newly identified RHOA, NPAS4 and SNARE pathways in DS Models

RHOA is a small GTPase protein acting through the activation of ROCK (RhoA Kinase) and the phosphorylation of the myosin light chain (MLC). Interestingly, RhoA was not found altered in the differential expression analysis; instead, it was a connecting node introduced to obtain a full connected PPI network. To ascertain whether RHOA pathway was altered in Dp1Yey, we checked the expression levels of two proteins of the pathway by qRT-PCR and WB (Fig. 4A). We found no changes in the expression of RHOA, converging with the transcriptomic analyses, but we detected a significant decrease of MLC phosphorylation (P-MLC) in the Dp1Yey hippocampi compared with control (Supplementary Material, Fig. S16). Thus, the RHOA pathways appeared to be downregulated in the Dp1Yey DS mouse model.

In our transcriptomics and network analysis, we found that Npas4 was misregulated in the Tg(Dyrk1a), Dp1Rhr and Ts65Dn models. We verified the downregulation of Npas4 and several other immediate early genes (IEGs: Arc, Fosb and Fos) in DS mouse models (Supplementary Material, Fig. S15). It is known that these IEGs are activated when light exposure is induced after a long light deprivation period (63,64). To confirm the impact of Npas4 downregulation in Tg(Dyrk1a) mice, we performed qRT-PCR experiments to determine the specific early and late response genes altered in the visual cortex after light deprivation and de novo light exposure at three time points (1, 3 and 7.5 h) (Fig. 4B). The results showed that Npas4 was clearly induced after light deprivation following 1 h of light stimulation but the expression level of Npas4 was higher in Tg(Dyrk1a) mice compared with control (Fig. 4C). We also took the opportunity to observe the expression of late response genes specific to inhibitory neurons (Frmdp3, Slc25a36 and Igf1,  Fig. 4D) and late response genes (Grp3 and Nptx2,  Fig. 4F). We found that gene expression was altered after 7.5 h of light stimulation in Tg(Dyrk1a). Interestingly, late response genes specific to excitatory neurons (Bdnf and Nrn1,  Fig. 4E) were not affected. The Snap25, Snap23 candidate genes found in our analysis showed an altered expression after 7.5 h of light stimulation while Dyrk1a and RhoA levels were not affected (Fig. 4F).

Discussion

In this study, we explored five DS mouse models carrying three copies of the region homologous to Hsa21 found on Mmu16, to decode the DS genotype–phenotype relationships and further investigated genetic interactions between different regions. To this end, we also assessed a transgenic model overexpressing one copy of Dyrk1a, plus two combinations of models (Dp5/Dp1, Dp5-Tg, Fig. 1), using a standardized behavioural pipeline focused on hippocampus-dependent memory processes; a process found impaired in people with DS.

In this parallel comparison, we observed several loci contributing to the alteration of different brain memory and control functions. We found that the spontaneous alternation in the Y-maze was altered in most of the models (except Dp5Yah alone; Fig. 1B). The minimal common genetic part of these lines was the overexpression of Dyrk1a and the result observed for transgenic Tg(Dyrk1a). Altogether, our previous results support DYRK1A as being the main driver of working memory defect (25) although another region could be involved in controlling spontaneous alternation (41).

Similarly, DYRK1A overdosage is a major cause for increased activity in the OF in Dp5/Tg(Dyrk1a) and Tg(Dyrk1a), but not in other models. While the Dp1Rhr Ds model was not affected, it can be hypothesized that another locus, located in the Cbr1-Fam3b interval, interferes with Dyrk1a overdosage in this model. The situation should be even more complex with additional genetic interaction. Indeed, no phenotype was observed for the distance travelled in the Dp1Yey, Dp5/Dp1, Dp5Yah, Dp1Rhr models. Thus the overexpression of Dyrk1a was able on its own, or combined with Dp5Yah, to induce increased activity while some loci that were not trisomic only in the Dp5/Dp1 model were able to suppress this effect, which can be reinduced by other trisomic loci specific to the Ts65Dn (see below). Altogether, our results suggest that many different genetic influences (at least three for the distance travelled) act on different behavioural variables in DS models.

Similar analysis of the NOR test with two distinct retention times highlighted several loci. The NOR test with 24 h of retention unravelled deficiency in most of the models, except in Dp5Yah and Dp5/Dp1, suggesting that there are at least two causative and two modifier loci (Fig. 5). Taken together, our results suggest that, depending on the variable observed in the behavioural test, several genetic interactions occur to build the link of behavioural phenotyping outcome in DS mouse models with loci, spread along Mmu16, including Dyrk1a. With 1 h of retention time, the NOR test pointed only to Dyrk1a overexpression, with at least one suppressing loci in the Dp1Rhr trisomic region.

Figure 5

Genotype correlation associated to behaviour phenotype in partial trisomic DS model. Here we highlight the duplicated region carried on each model with the corresponding syntenic region in the human chromosome 21 together with the main behavioral and transcriptomics results pointing to the existence of region specific phenotypes and functional alteractions. The black lines represents the duplicated syntenic regions to human chromosome 21 on each model (represented in the yellow line). The blue lines represents the behavioral results where no alteration was found, instead the red lines identifed the tests with deficits. Over the transcriptomics meta-pathways fucntional profile summary picture, in purple is highlighted upregulation whereas in pink downregulation. The intensity of the color stands for the number of pathways included on each meta-pathway from the total number of pathways found altered on each model.

As expected, we found deficits in the Ts65Dn mice similar to those previously published in the Y-maze, the OF, NOR, MWM and contextual FC (39,65) tests. Strikingly, thigmotaxis and time in the target quadrant in the probe test of the MWM were two variables modified only in the Ts65Dn trisomic model while Dp1Yey, which carries a duplication of the complete syntenic Mmu16 region, was less affected, as described before (66,67). Remarkably, the Dp5Yah mice, with a duplication from Cyyr1 to Clic6, displayed no deficits on their own in any of the tests performed and showed the lower number of DEGs in the hippocampi. Although several genes, Sod1, Olig1, Olig2, Rcan1 and Synj1, from the region were proposed as inducing an early DS cognitive phenotype (68–71), our results indicated that the Cyrr1-Clic6 region was not sufficient to induce by itself cognitive defects in DS models. Nevertheless, no defects were found in the Dp5/Dp1 mice, contrary to Ts65Dn in the OF, or they appeared less severe in the NOR test (see below). This indicates the existence of a key modulator in the Cyyr1-Clic6 region. The major behavioural alterations found in Ts65Dn probably result from the influence of different genetic factors: first the overdosage of genes homologous to Hsa21; then the presence of the freely segregating mini-chromosome (72), and also the trisomy of about 60 Mmu17 centromeric genes, non-homologous to Hsa21 (16) and overexpressed in the hippocampus of Ts65Dn mice. In particular, the overdosage of Tiam2 and Synj2, located on the Mmu17 centromeric region could exacerbate the effect of the overexpression of their respective paralogs Tiam1 and Synj1 (71,73). Strikingly, the transcriptomic analysis showed a different global disruption of genome expression in the Ts65Dn hippocampi, compared with the other trisomic models with segmental duplications. This was emphasized by the low correlation of DEGs and deregulated pathways between Ts65Dn, Dp1Yey and Dp5/Dp1. Overall, we can hypothesize that the suppression effect seen in the Dp1Yey model compared with the Ts65Dn, is due to a suppression effect of genes overexpressed and located upstream of Mrpl39, or to an enhancing effect of genes located on the non-homologous region in the Ts65Dn minichromosome, or to the freely segregating minichromsome in the Ts65Dn. New models and further analysis will be needed to test these hypotheses.

Similar to cognition, brain morphology was affected differently in DS models. As observed in people with DS, a global decrease in brain size was observed in Dp1Yey, Ts65Dn and Tg(Dyrk1a) while Dp1Rhr showed an increase in size in many brain regions compared with the other DS models. In addition, specific changes were found in several regions including the basal forebrain septum, a predominant source of cortical cholinergic input with an early substantial loss of basal forebrain cholinergic neurons; a constant feature of Alzheimer’s disease and other deficits in spatial learning and memory (74). Besides, only Ts65Dn presented an enlargement of the ventricles, which was previously associated with a decrease of cortical neurogenesis in the brains of the Ts1Cje and Ts2Cje mouse models (33).

Comparative genome wide expression profiling in the mouse hippocampus revealed that the entire dysregulation cannot be attributed to a single gene or region. The overall effect results from a complex interplay between certain trisomic overexpressed genes and other genes spread along the genome, evidenced by the fact that the majority of DEGs were not Hsa21 genes (Supplementary Material, Table S3). Additionally, we identified 34 trisomic DEGs (TEGs) with regulatory activity (Transcription factors, chromatin modellers, etc.) as Mir99a, Usp16, Erg or Rcan1 that may be involved in the changed regulatory landscape of the models. Indeed, RCAN1 and USP16 were found to be upregulated in human brain datasets (cerebrum and cerebellum)(49,75) and USP16 was also found as a DEG upregulated in heart and adult fibroblasts while MIR99A was found upregulated in adult fibroblasts (50). Nevertheless, the expression of DEGs was strongly correlated and conserved in Mmu16 based DS models. This is similar to the behavioural results obtained where related phenotypes were found in models carrying correlated partial duplications. Unexpectedly, Dp1Yey DEGs correlation was closer to Tg(Dyrk1a) than to Ts65Dn (42% against 25%) and there was a negative correlation between Dp3Yah & Dp5Yah (22%) and Tg(Dyrk1a) & Dp5Yah (13%). These correlations point to different gene dysregulations in these models and to the existence of epistasis with several regulatory trisomic genes countering the effect of genes in other trisomic regions.

We carried out an in depth functional annotation analysis to characterize 10 major meta-pathways with ribosome and mitochondrial functions, transcription and epigenomic regulation, and the synapse function categories highly affected. We also found a strong upregulation of genes involved in the Interferon beta pathway (Supplementary Material, Fig. S13B) as some interferon receptors were found upregulated in Mmu16 DS models. As such Ifnar2 and Il10rb were found upregulated in all the mice lines (except Tg(Dyrk1a)), pointing to a potentially critical role in interferon pathway dysregulation. The same phenomenon was observed with other genes like Irgn1, Ifit1, Ifit2 or Ndufa13. This upregulation of the interferon beta pathways was previously reported in the Ts1Cje mouse model (76,77) and linked to a possible increase of activity of the JAK-STAT signalling pathway, recorded here by the up regulation of Stat1.

The expression of genes involved in long-term synaptic potentiation (LTP) and synaptic plasticity were decreased in Dp1Yey, Dp5Yah, Dp1Rhr, Dp5/Dp1 models, respectively, corroborating previous reports on different DS mouse models and in vitro studies. The only upregulated pathways were the myelin sheath and SNARE, both found in Dp1Yey and Tg(Dyrk1a) models. Interestingly, models carrying the Dp1Rhr region duplication showed dysregulation mainly in synapse transmission, plasticity and LTP, while models carrying the Dp5Yah region duplication showed dysregulation associated with genes involved in stemness and differentiation. Together, the models with both Dp5Yah and Dp1Rhr duplicated regions were involved in post-synapse modulation and transmission. Thus, there seems both region-specific effects and region-dependent effects.

Moreover, the high intermodel inter-pathway connectivity approach showed six major sub-network biological cascades controlled by DYRK1A, GSK3β, NPY, SNARE complex, RHOA and NPAS4, which play a crucial role in the brain function. DYRK1A is a well-recognized driver of DS phenotypes and the target of several therapeutic approaches (25,38), which also interacts with GSK3β and NPY in DS models (25,78–82). Thus, we were pleased to detect these three pathways and surprised to notice how closely interconnected those sub-networks were. We also described new pathways with SNARE, RHOA and NPAS4 in models based on Mmu16. Interestingly, SNARE complex proteins were also found modified in a DS model for the region homologous to Hsa21 on Mmu17 (59). RHOA is a member of the RHO GTPase involved in several intellectual disabilities that affect dendritic structure in adult neurons (83–86), a phenotype also described in certain DS models (87–89) or linked to DYRK1A (90,91). Misregulation of Npas4 and IEGs was found in various DS mouse models and can induce abnormal regulation when activated in specific biological processes, such as cognition in the Ts65Dn model (92). Furthermore, the network analyses highlighted NPAS4 as a potential modulator of synaptic dysfunction via well-connected interactors. NPAS4 could affect the main altered biological cascades in addition to the GABA and NMDA receptors involved in the modulation of the excitatory/inhibitory balance of the brain (63).

The betweenness centrality index value is used to measure the interconnectivity of the network and showed that RHOA, DYRK1A, GSK3β, and their interactors were more closely knitted together and populated the central part of the network while SNARE, NPAS4 and NPY with their first- and second-layer interactors were more in the periphery of the network. This strong interconnectivity is of interest for two reasons: it makes the full network highly sensitive to targeted attacks against these proteins while, on the contrary, the network is robust against such attacks if they do not target several proteins simultaneously, for example during a drug trial. Thus, studying further closely connected altered genes and understanding their interactions could provide novel insights into the possible molecular mechanisms explaining why so many compounds, including DYRK1A specific kinase inhibitors, are capable of restoring learning and memory in DS models (11,38,93,94). Additionally, these nodes show a large number of connections. Indeed, using the betweenness index, these nodes can be seen to orchestrate the network and their interactors occupy its centre, illustrating their extreme importance for its stability in terms of network theory. Moreover, similar to the observation in DS people, where the affectation/severity of gene dysfunction varies from one individual to another, we propose here that different DS mouse models show different changes in the six signalling interconnected cascades affecting brain dysfunction, leading to similar behavioural phenotypes. Our observation may support the developmental instability hypothesis, which postulates that the non-specific triplication of a relatively small number of genes causes genetic imbalance with a significant impact on global gene expression. This hypothesis is in agreement with the study of rare individuals, carrying Hsa21 duplications, who displayed intellectual disabilities (5,7) and should be taken into account when therapeutic assays are planned. We suggest that preclinical observation in one partial trisomic mouse model should be replicated in more genetically complex models to test potential genetic influences (38,59,95). This is probably the limit of the model, since although behaviour and memory mechanisms are common between mice and humans, the complexity of the model is lower. Conducting the same studies in more complex animal models carrying all the trisomic genes homologous to Has21, would definitely permit better deciphering of genes having an impact on cognitive behaviour.

Taking advantage of DS mouse models, we investigated behaviour and cognition, brain morphology and hippocampal gene expression in a standardized and controlled manner. Our results with the partial duplication of the Mmu16 region homologous to Hsa21 are in agreement with human genetic analysis (4–10) and showed how multiple genetic interactions between different regions of chromosome 21 contribute towards altering the outcome of the behavioural, morphological and molecular/pathway phenotypes. We are now faced with the challenge of carefully dissecting all these genetic interactions. Nonetheless, we found that overlapping DS models show convergence in the biological cascades altered, observed via building PPI and regulatory networks, and centred on six main hubs: DYRK1A, GSK3β, NPY, SNARE, RHOA and NPAS4. We propose to name these hubs the centre of the DS biological cascade. Some of them have already been described as altered in certain DS models, and we validated two additional ones, RHOA and NPAS4. Thus, we have built a novel vision of existing altered gene–gene crosstalk and molecular mechanisms, with six specific highly interconnected DS hubs in mouse models. They may well prove essential in improving our understanding of DS neurobiology and making progress in therapy development.

Material and Methods

Mouse lines

The duplications of different Mmu 16 regions (Dp(16Lipi-Zbtb21)1Yey (or Dp1Yey), the Dp(16Hspa13-App)3Yah (noted Dp3Yah), the Dp(16Cbr1-Fam3b)1Rhr (or Dp1Rhr)) and BAC transgenic mice for Dyrk1a (Tg(Dyrk1a)) models were described previously (19,20,24,37). The genetic background of the DS lines carrying each duplication was pushed towards the C57BL/6J (B6J) genetic background for more than seven generations of backcrossing. The only exception was the trisomic Ts65Dn (Ts(1716)65Dn) mice, initially obtained from the Jax, which were kept on a F1 B6C3B genetic background (with the C3B line as a C3H/HeH congenic line for the BALB/c allele at the Pde6b locus (96)). The Dp(16Cyyr1-Clic6)5Yah (noted Dp5Yah) was generated by the in vivo TAMERE technology inserting loxP sites in App and Runx1 (see Supplementary information). In the Dp3Yah and Dp5Yah models, only two complete copies of App and Runx1 genes were expressed. The Dp5Yah line was crossed with the Dp1Rhr line in order to generate Dp5Yah/Dp1Rhr (also noted Dp5/Dp1) compound transheterozygote carrying a similar trisomic Mmu16 gene content to that of the Ts65Dn. Indeed, 15 Hsa21 homologous genes (Mrpl39, Jam2, Atp5j, Gabpa, App, Cyyr1, Runx1, Setd4, Mx2, Tmprss2, Ripk4, Prdm15, C2cd2 and Zbpb21) are not in three copies in Dp5/Dp1 compared with Ts65Dn. In addition 46 protein-coding genes located on the Mmu17 centromeric region in the Ts65Dn minichromosome (13,16) are not trisomic in the Dp5/Dp1. The Dp5Yah model was also combined with Tg(Dyrk1a) by crossing Dp5Yah/+ and Tg(Dyrk1a)/0 animals and generating the four genotypes (Dp1Yah, Dp5Yah, Tg(Dyrk1a) and [Dp5Yah; Tg(Dyrk1a)] noted here Dp5-Tg), to test specific interaction between Dp5Yah and Dyrk1a overdosage.

All the lines were maintained under specific pathogen-free conditions and were treated in compliance with the animal welfare policies of the French Ministry of Agriculture (law 87 848), and the phenotyping procedures were approved by our local ethical committee (Com’Eth, no. 17, APAFIS no. 2012-069).

Behaviour pipeline

A series of behavioural experiments were conducted in male mice with an age-range starting at 2.5 up to 7 months for the last test, as described in the Supplementary information. The tests were administered in the following order: Y-maze, OF, NOR (24 h), MWM and FC (contextual and cue). Behavioural experimenters were blinded to the genetic status of the animals. Separate groups of animals were composed for each line (as indicated in the Supplementary Material, Table S1). Several mouse models found defective for the NOR performed with 24 h of retention memory were also tested after 1 h of retention. The Dp5Yah crossed with Tg(Dyrk1a) was tested for Y-maze and NOR at 24 h. All the standard operating procedures for behavioural phenotyping have been already described (95,97–100) and are detailed in the Supplementary Information.

Magnetic resonance imaging

A dedicated cohort of animals at the age 102 ±7 days was anesthetized and perfused with 30 ml of room temperature 1× phosphate buffer saline (PBS) complemented with 10% (% w/v) heparin and 2 mm of ProHance Gadoteridol (Bracco Imaging, Courcouronnes, France) followed by 30 ml of 4% PFA complemented with 2 mm of ProHance Gadoteridol. Then, the brain structure was dissected and kept in PFA 4% 2 mm ProHance overnight at 4°C. The next day, each specimen was transferred into 1× PBS 2 mm ProHance until imaging.

Just prior to imaging, the brains were removed from the fixative and placed into plastic tubes (internal diameter 1 cm, volume 13 ml) filled with a proton-free susceptibility-matching fluid (Fluorinert® FC-770, Sigma-Aldrich, St. Louis, MO). Images of excised brains were acquired on a 7T BioSpec animal MRI system (Bruker Biospin MRI GmbH, Ettlingen, Germany). Images were reconstructed using ParaVision 6.0. An actively decoupled quadrature-mode mouse brain surface coil was used for signal reception and a 72-mm birdcage coil was used for transmission, both supplied by Bruker. The first protocol consisted of a 3D T2-weighted rapid-acquisition with relaxation enhancement. The second imaging protocol consisted of a 3D T2*-weighted Fast Low Angle (FLASH) sequence. The image matrix for both sequences was 195 × 140 × 90 over a field of view 19.5 × 14.0 × 9.0 mm3 yielding an isotropic resolution of 100 μm and treated and analysed for anatomical parameters, as detailed in the Supplementary Information.

Gene expression assay

Hippocampuses were isolated from DS trisomic models and their littermate controls (N = 5 per group) and flash frozen in liquid nitrogen. Total RNA was prepared using the RNA extraction kit (Qiagen, Venlo, The Netherlands) according to the manufacturer’s instructions. Sample quality was checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Gene expression analysis was carried out using GeneChip® Mouse Gene 1.0 ST arrays (Affymetrix, Santa Clara, CA, USA). All the procedures and analyses are detailed in the Supplementary Information.

Availability of data and materials

All the mouse lines are available in the Jax or the EMMA/Infrafrontier repository. Raw microarray data and re-analysed data have been deposited in GEO (Accession No. GSE149470).

Bioinformatic analysis

The gene expression profile of the mouse hippocampi isolated from Dp1Yey, Dp3Yah, Ts65Dn, Dp5/Dp1, Dp5Yah, Dp1Rhr and Tg(Dyrk1a) trisomic mouse models was analysed with a specific bioinformatics pipeline and controlled for quality prior to and after data pre-processing and normalization (see Supplementary Information in the detailed Material and Methods section). The DEGs were identified using a method based on fold-change (FC) rank-ordering statistics (FCROS) (101). In the FCROS method, k pairs of test/control samples are used to compute FC. For each pair of test/control samples, the FCs obtained for all genes are ranked in increasing order. Resulting ranks are associated with genes. Then, the k ranks of each gene are used to calculate a statistic and resulting probability (f-value) used to identify the DEGs after fixing the error level at 5% false discovery rate (FDR).

We performed the functional differential analysis using GAGE (102) and grouped all the pathways into 10 functional categories (noted meta-pathways). Functional intermodel meta-pathway connectivity was studied by identifying the genes shared between pathways and models inside the same meta-pathway. Then, to assess gene connectivity, we built a minimum fully connected PPI network (noted MinPPINet) of genes known to be involved in synaptic function as they were associated with synaptic pathways via GO (103) and KEGG databases (104). Furthermore, regulatory information was added to build the final RegPPINet. We used the betweenness centrality analysis to identify hubs, keys for maintaining the network communication flow. The relevance of the connecting nodes was further predicted by the machine learning algorithm Quack (105). Finally, we computed 100 000 random networks with a similar degree, to assess if the likelihood of observing such connectivity in the DS network was more than one can expect by chance using statnet and sna R packages (https://cran.r-project.org/web/packages/statnet/index.html;  https://cran.r-project.org/web/packages/sna/index.html). The full list of R packages used can be found in the Supplementary Material, Table S7.

Western blot

The expression levels of the RHOA protein and MLC phosphorylation by the MLC Kinase part of the RHOA pathway were analysed using western blot in five animals Dp1Yey and five control (wt) littermates (see Fig. 4G, Supplementary Material, Fig. S16, and Supplementary Information). We used the following primary antibodies: anti-RHOA (2117, Cell Signalling, USA, 1:1000), anti-pMLC (Thr18/Ser19 #3674, Cell signalling, Boston, MA, USA, 1:1000) and mouse monoclonal Anti-β-Actin-Peroxidase antibody (A3854 Sigma, 1:150 000); and HRP conjugated Goat anti-Rabbit IgG secondary antibody (A16096, Invitrogen, France). Protein signals were visualized with Amersham™ Imager 600 and were quantified using ImageJ and statistical analysis using Sigma Plot. The relative amount of RHOA and p-MLC proteins was calculated as the ratio of the signal detected for each protein of interest compared with the β-actin signal detected and normalized by the mean signal of the wt samples.

Visual stimulation

Mice raised in a standard light cycle were housed in constant darkness for 2 weeks. Then, animals in the light-exposed condition group, were consecutively exposed to light for 0, 1, 3 and 7.5 h before being sacrificed. The animals belonging to the dark-housed condition group were sacrificed in the dark. After euthanasia, their eyes were enucleated before visual cortex dissection in the light and flash frozen in liquid nitrogen. cDNA and quantitative PCR were performed as indicated in the Supplementary Information. The Ct values were transformed to quantities by using the comparative Ct method. Hence, all data were expressed relative to the expression of the most expressed gene. These relative expression levels, were normalized with Genorm by keeping the more stable reference genes (106). To calculate fold-induction, the relative quantity of gene expression at each time point was divided by the mean of the relative level of gene expression of dark-housed mice for the corresponding genotype. The mean and standard error were calculated at each time point from these fold-induction values.

Statistical analysis

All data are expressed as mean group value ± standard error of the mean (SEM) or as box plots with the median and quartiles. For each data set, we analysed if the data were normally distributed using the Shapiro–Wilk test and Quantile–Quantile plots (Supplementary Material, Fig. S2) and the homogeneity of variances by the Brown–Forsy test. Differences between groups were inferred by one-way ANOVA (OF) and ANOVA for repeated measures, or we performed the Kruskal–Wallis non-parametric test in the case of datasets where the assumptions of normality or homogeneity of variances were not fulfilled. The post-hoc tests (Fisher’s LSD Method) were conducted only if the F parameter in ANOVA achieved a level of 0.05. All the behavioural analysis results can be found in the Supplementary Material, Table S1. For the MRI data, intergroup comparisons on region-based data were conducted on the normalized volumes (i.e. ratio between the volume of the structure and the whole brain volume) of each segmented structure using the Student’s t-test while correcting by multiple testing setting up an FDR correction.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.