Mouse models of immune dysfunction: their neuroanatomical differences reflect their anxiety-behavioural phenotype

Extensive evidence supports the role of the immune system in modulating brain function and behaviour. However, past studies have revealed striking heterogeneity in behavioural phenotypes produced from immune system dysfunction. Using magnetic resonance imaging, we studied the neuroanatomical differences among 11 distinct genetically modified mouse lines (n = 371), each deficient in a different element of the immune system. We found a significant and heterogeneous effect of immune dysfunction on the brains of both male and female mice. However, by imaging the whole brain and using Bayesian hierarchical modelling, we were able to identify patterns within the heterogeneous phenotype. Certain structures—such as the corpus callosum, midbrain, and thalamus—were more likely to be affected by immune dysfunction. A notable brain–behaviour relationship was identified with neuroanatomy endophenotypes across mouse models clustering according to anxiety-like behaviour phenotypes reported in literature, such as altered volume in brains regions associated with promoting fear response (e.g., the lateral septum and cerebellum). Interestingly, genes with preferential spatial expression in the most commonly affected regions are also associated with multiple sclerosis and other immune-mediated diseases. In total, our data suggest that the immune system modulates anxiety behaviour through well-established brain networks.


INTRODUCTION
The immune system plays an important role in brain-behaviour interactions. During pregnancy, for example, maternal illness increases the risk of nervous system disorders in offspring [1][2][3][4][5] and maternal IL-6 cytokine concentration is associated with altered neurodevelopment in offspring [6][7][8]. In rodents, in utero exposure to maternal immune activation (MIA) affects brain structure and behaviour [9][10][11]. Immune system effects on the brain and, ultimately, behaviour are heterogenous and multifaceted. To better understand the complex pathophysiology responsible for immune dysregulation impacting brain development, the role of specific immune components has been investigated using targeted experiments in mice. For example, anxiety-like behaviour is modulated by the production of various cytokines [12][13][14], adaptive immune response [15][16][17][18], nitric oxide synthesis [19], and mast cell activity [20].
While the role of specific immune components is an area of active investigation, existing studies generally examine brain regions or immune components in isolation, thereby sacrificing the ability to observe broader patterns that may explain the observed heterogeneity. Structural magnetic resonance imaging (MRI) is a useful methodology for studying relationships between genetics, brain regions, and behaviour with several studies demonstrating a strong link between the volume of brain structures and organism behaviour [21,22]. As structural MRI has the benefit of whole-brain coverage without the constraint to identify a region-of-interest (ROI) ahead of time, it allows the identification of novel, unanticipated phenotypes in mouse models of human disease. For example, we previously showed using ex vivo MRI that mice lacking all functional T cells demonstrated a loss of sexual dimorphism in many brain regions [15]. We corroborated this finding with a loss of sexual differentiation in activity-related behaviours, highlighting the importance of T cells for the development of sex differences in neuroanatomy and behaviour. However, such studies remain rare and the brain-wide effects of alterations in other elements of the immune system are unknown. High-throughput data acquisition and analysis of structural MRI allows for the investigation of similarities and differences across a range of immune system mutations.
Thus, here we used whole-brain MRI to explore the effect of immune system mutations on the brains of male and female mice. We imaged the brains of 11 genetically engineered mouse lines, each deficient in different elements of the immune system, using a standard pipeline for specimen preparation, data acquisition, and image registration. Brain structure was highly susceptible to immune dysfunction and we found regional variations in this susceptibility. The most susceptible regions tended to express genes associated with immune-mediated diseases. We also identified that strains with similar anxiety-behavioural phenotypes displayed analogous neuroanatomical effects. Lastly, we identified candidate brain regions and networks that may be involved in modulating the anxiety behaviours seen in mouse models of immune dysfunction.

METHODS Mice
The study consisted of 371 mice (Table 1) from 11 different genetically engineered strains (Jackson Labs, Bar Harbor, Maine, USA)-which encapsulate the different components of the immune system (Supplementary Fig. 1). Sample size was chosen to detect 4% hippocampal volume (Type I error = 0.05, power = 80%) [23]. All mice were housed in individually ventilated cages (10 air changes/hour at 20-22°C), and provided with ad libitum water and diet#2918 (Envigo, Indianapolis, Indiana, USA). All breeding colonies were maintained using homozygous mutants with two exceptions: Cxcr2 colony was maintained using heterozygous mutants due to the fragility of homozygotes [24], and Rag2 mice were ordered directly from Jackson Labs. Genotypes were determined from tail biopsies using real-time PCR (Transnetyx, Cordova, Tennessee, USA) or fur colour in the case of the Kit strain. The control group consisted of C57Bl/6J mice and Cxcr2 wild-type littermates [25]. The two wild-type control strains were pooled for analyses as there were no significant neuroanatomical differences between them [26]. On P65±3 (mean±max range), mice were sacrificed for fixation (physical measurements reported in Supplementary Table 1 and Supplementary Fig. 2). Fixed brain samples were prepared for MRI using a previously described fixation protocol [27]. Briefly, mice were perfused with a solution of phosphatebuffered saline (Wisent Incorporated, Saint-Jean-Baptiste, Quebec, Canada), gadoteridol (Bracco Diagnostics Incorporated, Monroe Township, New Jersey, USA), and heparin (Sandoz Canada Incorporated, Mississauga, Ontario, Canada) followed by gadoteridol and 4% paraformaldehyde (PFA) (Electron Microscopy Sciences, Hatfield, Pennsylvania, USA). After decapitation, brains (kept within the skull) were stored at 4°C in a solution of PFA, gadoteridol, and sodium trinitride (Fisher Scientific, Mississauga, Ontario, Canada). Brains were imaged 7-8 weeks postmortem [28].
All animal experiments were approved by The Centre for Phenogenomics (TCP) Animal Care Committee (AUP 17-0175H) in accordance with recommendations of the Canadian Council on Animal Care, the requirements under the Animals for Research Act, RSO 1980, and the TCP Committee Policies and Guidelines. As subjects were assigned to groups based on sex and strain, experimenters were not blinded to subject and subjects were not randomised.

Frequentist statistics
For each structure, volume was fit using a linear model with predictors of sex, strain, and their interactions; and partial F-test was used to assess significance. Multiple comparisons were corrected using false-discovery rate (FDR) [40]. Post hoc simulations to assess power were also performed ( Supplementary  Fig. 6). Variances were similar between groups (Supplementary Table 4).

Bayesian statistics
Structure volumes were standardised using Z-score. An anatomical hierarchy [41] was used to reduce the number of structures to 95 bilateral brain regions (Supplementary Figs. 7 and 8). A Bayesian hierarchical model (BHM) [42] was used to explore the effects of strain on neuroanatomical volumes of all 95 structures. It contained global and structure-specific predictors of sex, strain, and their interaction, and individual-specific intercept (priors and healthy model diagnostics given in Supplementary  Table 2). Similar results were seen when refitting the model with different priors for the correlation matrices. For each structure, the posterior distribution was used to compute the effect-size (d) distribution.
The posterior also provided the distribution of neuroanatomy phenotype as beta-values for each structure and strain (referenced to wild-type). Hellinger distance was used to assess the dissimilarity of neuroanatomy phenotype between each pair of strains. Strains were annotated with anxiety phenotype based on existing literature (increased, decreased, or unchanged; Supplementary Table 3). Since Cxcr2 strain did not have an anxiety phenotype reported in the literature, it was assumed to be unchanged. To assess whether strains with similar anxiety phenotype have lower Hellinger distances, permutation testing was performed (10 5 iterations). Similar results were seen after excluding the Cxcr2 strain and when using a bootstrap-based networking procedure [21]. To determine the effect of anxiety for each structure, the posterior distribution was used to calculate the anxiety effect size (η 2 ) distribution. Hellinger distance, η 2 , and d statistics are illustrated in Supplementary Fig. 9.

Preferential gene expression and disease enrichment analysis
Using previously published methods [43,44] and the Allen Brain Institute (ABI) gene expression atlas [41], we explored the relationship between brain regions with altered volume and gene expression. Determined by the BHM, the top 25 brain structures with the highest absolute effect size (i.e., all structures with |d| > 1) constituted the ROI for this analysis. Preferential expression for each gene was evaluated using fold-change: defined as mean expression in ROI divided by mean expression in the brain.
For disease enrichment and Gene Ontology (GO) enrichment analysis, all genes were part of the background set, while the target set was the top 4000 genes with the highest preferential expression. Similar results were seen with target sets of the top 3000 and 5000 genes. For disease enrichment analysis, mouse genes were annotated with human diseases using two databases: NCBI [45] (mapping mouse genes to homologous human genes) and DisGeNET [46] (annotating human genes with associated diseases). GO enrichment analysis was conducted using GOrilla [47]. Significance was assessed using hypergeometric tests with FDR correction. Differences in probability density were assessed for significance using the Kolmogorov-Smirnov test. ABI developmental gene expression atlas [48] was used to investigate gene expression patterns through development. Expression was clustered using k-means with k = 4 determined by elbow-method [49].

RESULTS
Immune dysfunction affects male and female neuroanatomy A significant effect of strain was observed on the volumes of nearly all brain structures (Fig. 1A), and this effect was similar for male and female mice ( Supplementary Fig. 10). While the effect of strain was highly significant, the magnitude and direction of these effects showed a great deal of heterogeneity ( Fig. 1B for females, Supplementary Fig. 11 for males, Supplementary Table 4 for complete list). For example, mutants lacking cytokines IL-6, IL-10, and IL-18 had a similar phenotype of larger cerebellum compared to wild-type. However, they differed in phenotypes of the left frontal association cortex with IL-10 being smaller than wild-type, and IL-6 and IL-18 being larger.
To find important patterns through these heterogeneous effects, a BHM was fitted to the neuroanatomy data. Consistent with the frequentist analysis (Fig. 1A), all brain structures had a high probability of having moderate effect-size magnitudes (95% credible interval exceeded 0.50 for all structures). However, certain structures were particularly sensitive to strain and had a high chance of having large effect sizes ( Fig. 2A). For example, the midbrain (Fig. 2B), corpus callosum (Fig. 2C), dorsal striatum (Fig. 2D), and thalamus ( Fig. 2E) Fig. 1 Immune system mutations have a highly heterogeneous effect on mouse brain anatomy. A Nearly all brain structures showed a significant effect of strain evaluated using F-statistics from ANOVA. B The directional effect in females of the various mutant strains relative to the wild-type strains is visualised using t-statistics and shows a heterogeneous neuroanatomical phenotype. Regions larger or smaller in mutants relative to wild-type are given maroon-pink and blue-turquoise colours, respectively, if effects are <5% FDR. Saturated colours represent effects <0.01% FDR.
probabilities of having large effect-sizes across the 11 strains evaluated (more examples in Supplementary Fig. 12).
Neuroanatomical differences among strains cluster by anxiety phenotype Given that immune strains had heterogeneous neuroanatomical phenotypes, we investigated whether clustering these phenotypes could reveal important underlying biology. Because prior literature has shown that immune dysregulation causes heterogeneous anxiety-like behaviours (Supplementary Table 3), we assessed whether there was a relationship between a strain's neuroanatomical endophenotype and anxiety behaviour reported in literature. We found that strains with similar anxiety-behavioural phenotypes had similar neuroanatomical endophenotypes (Fig. 3A  inset). This association is particularly driven by strains with increased anxiety behaviours, as visualised in a network plot (Fig. 3A). We then assessed which brain structures in particular drive this association with anxiety behaviour. Nearly every structure had at least a moderate effect-size of anxiety (Fig. 3B); the culmen and lateral septum were plotted as representative structures ( Fig. 3C; more examples in Supplementary Fig. 13). As the Cxcr2 mutant strain did not have a documented anxiety phenotype, we assumed this strain had no anxiety-like behaviours for this analysis; repeating this analysis excluding this strain resulted in similarly strong grouping by anxiety phenotype.
Affected neuroanatomy has preferential spatial expression of genes associated with immune-mediated disease We sought to link the neuroanatomical findings with immune dysfunction by comparison with genome-wide gene expression patterns from the ABI database [41]. For this purpose, an ROI was defined from the MRI results by selecting the set of 25 structures with the highest median effect-size magnitude (shown in Fig. 4C first row). Spatial gene expression analysis [43,44] was used to identify genes that were preferentially expressed within this ROI using a fold-change measure (average expression in ROI divided by average expression in the brain). Gene ontologies of preferentially expressed genes were enriched in many biological processes associated with white matter: ensheathment of neurons (GO:0007272, p < 10 −6 , FDR < 10 −3 ), axon ensheathment (GO:0008366, p < 10 −6 , FDR < 10 −3 ), and myelination (GO:0042552, p < 10 −6 , FDR < 10 −2 ). Using the DisGeNet [46] and NCBI [45] databases, we then identified diseases in which these preferentially expressed genes have been implicated. Many of the disease terms with significant enrichment have known or suspected immune-mediated pathophysiology, such as Parkinson's disease (p < 10 −8 , FDR < 10 −4 ), Alzheimer's disease (p < 10 −8 , FDR < 10 −4 ), and multiple sclerosis (MS) (p < 10 −4 , FDR = 0.016).
As there is strong evidence that MS is an immune-mediated disease [50], we decided to explore it further. Consistent with the previous enrichment analysis, genes associated with MS tended to have higher expression in brain regions sensitive to immune dysfunction (Fig. 4A). We also clustered the expression patterns of MS-associated genes in the ROI through the course of postnatal development [48] and found four main clusters of genes (Fig. 4B). Two clusters had low expression in adulthood with expression peaks in early life: Cluster 1 had a peak at~P4 and Cluster 2 peaked later~P14. The other two clusters had high expression in adulthood: Cluster 3 and 4 had high expression starting from~P14 and~P56, respectively. All gene expression analysis results are provided in Supplementary Tables 5-8. For each pair of strains, the dissimilarity of endophenotypes was assessed using Hellinger distance and visualised using a network (thicker edges imply greater similarity). A Strains with increased anxiety behaviours (red nodes) had similar endophenotypes and clustered together in the network. A similar pattern was seen for unchanged anxiety behaviours (dark purple), but not decreased anxiety behaviour (blue). Cxcr2 anxiety phenotype is not known and assumed unchanged (light purple). The inset plot shows that pairs of strains (represented as dots) with similar anxiety phenotype had similar neuroanatomy endophenotypes (p < 0.01 from permutation testing). B The effect size (η 2 ) for the anxiety grouping was computed for each structure. The green colour bar represents at least medium effect sizes and saturates for large effect sizes. C The culmen and lateral septum were chosen as representative examples to illustrate large effect-size for anxiety grouping. The 95% credible interval of predicted volume for each strain (grey bars) and anxiety phenotype (coloured bars) are shown. Fig. 4 Brain regions susceptible to immune system mutations have a preferential spatial expression of genes involved with multiple sclerosis (MS). The top 25 brain structures with the highest effect-size magnitudes are shown in C (first row) and constitute the region-ofinterest (ROI). For all genes in the mouse genome, preferential spatial expression was assessed using a fold-change measure (i.e., gene expression signal in ROI relative to the whole brain). A Genes associated with MS (solid line) had significantly higher expression in the ROI compared to the genome (two-sided KS test D = 0.05, p < 10 −3 ). B MS genes showed four different clusters of temporal expression within the ROI. The shaded region represents the 95% confidence interval. C A representative example from each cluster was chosen to visualise gene expression signals within the ROI over the course of neurodevelopment. Each example was closest to its respective cluster's centroid.

DISCUSSION
Using neuroanatomical imaging with ex vivo MRI, we confirm heterogeneous and widespread effects of the immune system on the brain using a selection of immune-related genetically engineered mouse strains. To identify patterns underlying this heterogeneity, we studied neuroanatomical phenotypes on three scales: analysing each strain individually, analysing groups of strains unified by a common anxiety-like behavioural phenotype, and aggregating across all strains. When analysing strains individually, we found results consistent with other studies. For example, we did not find altered morphology of the dorsal hippocampus in Rag1 mutants, which is consistent with histology studies [51]. We found thalamic changes in mast cell-deficient Kit mutants, consistent with mast cell migration into the thalamus during development [52]. Phenotypes also showed some degree of lateralisation that has previously been seen in autism models [53,54]. MRI-detectable volume changes are poorly specific to evaluation of biological determinants: depending on the context, studies have implicated glutamate concentration [55], neurogenesis [56], cellular composition [57], and axonal/dendritic processes [58,59]. By identifying phenotypes in each strain at mesoscale resolution, we hope targeted investigations can be performed in future studies to elucidate cellular mechanisms.
Analysis of the heterogeneous structural findings for association with behavioural phenotypes revealed important insights and suggests that neuroanatomical endophenotypes cluster significantly with anxiety-like behavioural phenotypes. Previous use of mouse models has revealed several neural circuits associated with producing anxiety behaviours [60,61], but it is not known which of these neural circuits are responsible for driving anxiety behaviours in models of immune dysfunction. We found that alterations in the lateral septum and hypothalamus cluster with anxiety; and activation of this circuit can promote persistent anxiety [62]. We also found that the cerebellum and midbrain were correlated with anxiety. Studies have shown the important role of these structures in fear-and anxiety-related brain networks [63,64]. It is unclear whether volume changes are associated with increased or decreased activation of a neural circuit [65]; however, these structures are promising candidates for future studies investigating anxiety behaviours in immune dysfunction. Interestingly, anxiety disorders are a notable comorbid condition in neurological disorders associated with immune dysfunction. For example, anxiety disorders are three times more prevalent in MS patients than in the general population [66]. The prevalence of anxiety disorders in Parkinson's is 31% [67].
Immune dysfunction plays an important role in many neurological disorders, but the impact of the immune system is difficult to study because the aetiology is highly polygenic. For example, MS is associated with over 1000 genetic risk variants [46], but many of which are associated with immune function [68]. Although requiring further validation, our data provide a potential new avenue for examining the basis of these interactions. By studying the neuroanatomy of several mouse models, our methodology may be useful in finding patterns underlying polygenic immune dysfunction. The corpus callosum, thalamus, striatum, and midbrain were amongst the most-affected structures across all mouse models, which reflects neurological findings in autoimmune disorders [69][70][71][72]. Furthermore, several genes associated with myelination are expressed in these regions, which is consistent with molecular studies conducted on MIA model primates [73]. Preferentially expressed genes were also homologous to human genes associated with MS. This finding suggests that expression patterns of MS genes could be followed through neurodevelopment to identify when their expression peaks, pointing to possible developmental timewindows for further study. Our data support the hypothesis that, while there are many genetic variants associated with immune system dysregulation, certain brain regions may be more sensitive to alterations in the immune system. Thus, further study of mouse models with targeted immune system mutations could uncover potential insights regarding neurological symptoms in polygenic immune-mediated disorders.
There are limitations in our work that are important to discuss. Our study focuses solely on static loss-of-function mutations in isolated immune system components, which is a simplification used throughout the existing literature [12][13][14][15][16]. However, it is important to consider that triggering immune system dysfunction in different windows of neurodevelopment can have differing effects [11]. Thus, our work cannot identify when in the maturation process anatomical phenotypes emerge. To capture these dynamic effects, our ex vivo imaging study design could be extended to have a longitudinal in vivo imaging component. Although in vivo MRI generally has inferior performance compared to ex vivo MRI studies of the same animals [74], the additional information could help identify important periods in neurodevelopment where certain immune system components have the greatest effects. Our analysis of behaviour was limited to anxiety because it was best characterised by currently available literature. Social behaviour has also been shown to be reduced by MIA [75] that may suggest social avoidance as a behaviour to reduce susceptibility to pathogen infection. Initiatives like the International Mouse Phenotyping Consortium could be tremendously useful as they would have consistent behavioural assays for every single-gene mutation in the mouse [26]. This would allow future studies to not only investigate the association between immune mutations and anxiety, but social behaviours as well.
In summary, we imaged the neuroanatomy of 11 different mouse mutants, each deficient in particular components of the immune system. By using a consistent high-throughput protocol for data acquisition and analysis, we characterised the diverse effects of immune dysregulation on the brain and identified patterns underlying this heterogeneity. The observed neuroanatomical differences in the mouse models of immune dysfunction are clustered by their anxiety phenotype, recapitulating known brain networks implicated in modulating anxious behaviours. Exploratory analysis of gene expression shows that the brain regions that are most affected by immune dysfunction have a preferential spatial expression of genes associated with MS and other immune-mediated diseases.