Gut microbiome partially mediates and coordinates the effects of genetics on anxiety-like behavior in Collaborative Cross mice

Growing evidence suggests that the gut microbiome (GM) plays a critical role in health and disease. However, the contribution of GM to psychiatric disorders, especially anxiety, remains unclear. We used the Collaborative Cross (CC) mouse population-based model to identify anxiety associated host genetic and GM factors. Anxiety-like behavior of 445 mice across 30 CC strains was measured using the light/dark box assay and documented by video. A custom tracking system was developed to quantify seven anxiety-related phenotypes based on video. Mice were assigned to a low or high anxiety group by consensus clustering using seven anxiety-related phenotypes. Genome-wide association analysis (GWAS) identified 141 genes (264 SNPs) significantly enriched for anxiety and depression related functions. In the same CC cohort, we measured GM composition and identified five families that differ between high and low anxiety mice. Anxiety level was predicted with 79% accuracy and an AUC of 0.81. Mediation analyses revealed that the genetic contribution to anxiety was partially mediated by the GM. Our findings indicate that GM partially mediates and coordinates the effects of genetics on anxiety.


Results
Computational pipeline for mouse behavior characterization. We used the light-dark box to assess the variation in anxiety-like behavior in 445 mice (228 female; 217 male) representing 30 CC strains with video recording of the light compartment. The number of mice for each strain ranged from 7 to 25 (Table S1). A computational pipeline was developed to quantify anxiety levels by tracking mouse location ( Fig. 1A-D). Seven anxiety-related phenotypes were extracted from each video file: number of full and partial transitions between light and dark, speed, distance traveled and total time spent in the light compartment, average time spent in light for each transition and latency to first transition into the dark compartment (Table 1) Anxiety related behavior phenotypes vary across CC strains. The travelling distance in the light compartment ranged from 1649mm (CC039) to 16311mm (CC004), the number of full transitions ranged from 2 (CC039) to 15 (CC042), the number of partial transitions ranged from 8 (CC011) to 30 (CC040), the total time spent in the light compartment ranged from 28.58s (CC039) to 159.84s (CC004), the average speed in the light compartment ranged from 39.54mm/s (CC010) to 95.24mm/s (CC004), the latency time for the first transition ranged from 6.44s (CC030) to 86.12s (CC032), and the average time in light ranged from 6.44s (CC030) to 86.18s (CC032) (Fig. 1E, Table S1). No significant difference was observed between male and female mice from any of the CC strains (adjusted p>0.05; Fig. S2). We observed significant variation in these seven anxiety-related phenotypes across different CC strains, strongly suggesting that host genetics influences anxiety.
Anxiety level assessment based on seven anxiety related phenotypes. In order to classify mice into different anxiety level, we combined the anxiety-related phenotypes of our mouse cohort and employed consensus clustering to obtain anxiety-related subgroups using different numbers of clusters (K = 2, 3 and 4; Fig. 2A). Consensus clustering visualizes the consistency by which each sample is assigned to a specific cluster. The cumulative distribution function (CDF) and change in the area under the curve for CDF at different values of K suggests maximum stability at K=4 (Fig. 2B,C). However, visual inspection of the consensus matrices, showed that dividing the mouse cohort into two subgroups (K= 2) corresponding to low and high anxiety levels resulted in the most consistent matrix ( Fig. 2A). Based on the result of the consensus clustering, we have divided the mouse cohort into two groups: high anxiety (HA) and low anxiety (LA). The seven anxiety-related phenotypes were significantly different between the two groups (Mann-Whitney U-test; FDR < 0.01; Fig. 2D), and multivariate logistic regression analysis indicates that five out of seven anxiety-related phenotypes independently contribute to the anxiety classification (Table S2; p<0.05). All mice from strains CC013, CC004 and CC008 were classified as LA, whereas all mice from strains CC010, CC038 and CC051 were classified as HA (Fig. S3). For each of the remaining strains, individual mice were assigned to either LA or HA suggesting incomplete penetrance of the anxiety-like phenotype. Data from all mice was included in downstream analyses.   CC003  CC002  CC004  CC016  CC019  CC021  CC009  CC028  CC041  CC040  CC052  CC038  CC042  CC051  CC055  CC039  CC057  CC059  CC061  CC033  CC032  CC036  CC037  CC008  CC011  CC010  CC013  CC026  CC030   12 CC001  CC003  CC002  CC004  CC016  CC019  CC021  CC009  CC028  CC041  CC040  CC052  CC038  CC042  CC051  CC055  CC039  CC057  CC059  CC061  CC033  CC032  CC036  CC037  CC008  CC011  CC010  CC013  CC026 3A; Tables S3 and S4). Gene Ontology analysis revealed that 141 genes were significantly enriched in biological processes related to neuronal function including synapse assembly, neuron fate specification and presynaptic endocytosis (p<0.05; Fig. 3B). Our screen identified 62 genes known to be associated with anxiety, behavioral alterations and neurodevelopment of which 40 genes show expression in the brain based on in situ hybridization data from the Mouse Brain Atlas (Allen Brain Atlas) ( Table S4). For instance, allele variants of NTRK3, PPP2R2B, and ESR1 were associated with anxiety in humans 32-34 , knock-out mice of Cacna1h, Rapgef2, Clstn2, and Tnr exhibited abnormal anxiety-like behaviors [35][36][37][38] , and Isl1, Abl2, Dlgap1 and Csmd2 knock-out mice showed alterations in neurodevelopment and behavior [39][40][41][42] . In addition to these 62 known genes, our screen identified 79 genes not previously associated with anxiety, which includes 38 genes that show expression in the brain based on in situ hybridization data from the Mouse Brain Atlas (Allen Brain Atlas) ( Table S4). The spatial gene expression data suggests that these 38 genes may play a role in anxiety.

Associations between microbiome and anxiety-like behavior.
To investigate the association between specific microbes and anxiety-like behavior, we collected fecal samples from 30 CC strains for 16S ribosomal RNA profiling. Sequence reads were mapped to 5761 OTUs corresponding to 71 bacterial families 30 . We identified five families with significantly different abundance levels between high anxiety and low anxiety groups (Mann-Whitney U-test; FDR < 0.1; Fig. 4A). We observed higher relative abundances of Ruminococcaceae, Clostridiaceae, Clostridiales_family unknown, and lower relative abundances of Bacteroidaceae and Bac-teroidales_S24-7 in HA mice compared to LA mice (Fig. 4A). Logistic regression confirmed that these families were significantly correlated with anxiety (FDR < 0.01; Fig. 4B). Abundance levels of Ruminococcaceae and Clostridiales_family unknown were positively correlated with anxiety level, while the abundance of Bacteroidaceae, Clostridiaceae and Bacteroidales_S24-7 were negatively correlated with anxiety level ( Fig. 4B; FDR < 0.05). At the OTU level we found that the abundance level of 327 OTUs was significantly different between the high anxiety and low anxiety groups (Mann-Whitney U-test; FDR < 0.05; Table S5). For example, we observed decreased abundance of bacteria from the genus Bacteroides (family Bacteroidaceae) and increased abundance of several genera from the family Ruminococcaeae including Ruminiclostridium, Ruminococcaceae_UCG-014 and Oscillibacter in high anxiety compared to low anxiety mice.
To determine the importance for each of the five families related to anxiety, we performed random forest classification analysis and show that Bacteroidaceae and Ruminococcaceae contribute most to the anxiety phenotype (Fig. 4C). Furthermore, random forest classification based on the five families related to anxiety level (i.e., HA and LA) resulted in a predictive accuracy of 79% with an AUC around 0.81 (Fig. 4D), where the predictive power is statistically identical with the one derived from all microbiome features (Fig. 4E). Our results demonstrate an association between microbiome features and anxiety-like behavior.
Gut microbiome partially mediates the effects of genetics on anxiety-like behavior. We then performed mediation analysis to investigate whether genetic variants indirectly contribute to anxiety-like behavior by controlling the abundance of the five families associated with anxiety. Mediation analysis is a statistical model to determine whether the relationship between two variables (genetic variants and anxiety) is mediated through a third variable (gut microbiome). We first identified 12,368 genetic variants significantly (Mann-Whitney U-test; p<1E-09) associated with any of the five microbial families. We then selected only genetic variants significantly associated with both the abundance of any of the five microbial families and anxiety-like behavior. Mediation analysis revealed that Ruminococcaceae, Clostridiaceae, Bacteroidaceae and Clostridiales (family unknown) function as mediators for the effect of 17 SNPs within 7 genetic loci on anxiety (Fig. 5, Table S6; FDR <0.05). These analyses indicated that the effect of genetic variations on anxiety is at least partially mediated by the gut microbiota, and therefore suggested that the gut-brain axis plays important roles in anxiety disorders.
Anxiety related genes show significant overlap with human GWAS for psychiatric conditions. To discover candidate genes with potential impact on human anxiety, we compared anxiety-related mouse genes (141 genes) to a compiled list of human genes associated with seven psychiatric conditions that   (Table S7). We observed significant overlap between the mouse anxiety associated genes and genes associated with attention deficit hyperactivity disorder, depression, feeling nervous, and neuroticism ( Fig. 6; p<0.05). We found that 25 out of our 141 anxiety-related mouse genes were associated with one or more human GWAS for psychiatric conditions. Interestingly, STAG1 and SORCS3 were both found in four phenotypes. STAG1 was associated with autism spectrum disorder, feeling nervous, feeling worry and neuroticism, whereas SORCS3 was associated with Alzheimer, depression, feeling nervous and neuroticism. We thus conclude that the candidate mouse genes identified in this anxiety associated study exhibit significant relevance with human psychiatric conditions.

Discussion
In our study, a systematic genome and metagenome analysis on anxiety-like behavior was performed on 445 mice across 30 genetically defined CC strains to identify the effects of host genetics and gut microbiota and their interaction on anxiety. Our findings have the potential to provide insights on the mechanisms of host-microbe interactions related to anxiety. We further demonstrate that genetic effects on anxiety are partially mediated through modulating the abundance of specific gut microbes, suggesting links between host genetics and anxiety via intestinal health.  www.nature.com/scientificreports/ more accurately define anxiety level subgroups. Here, we integrated seven anxiety-related behavior phenotypes and utilized consensus clustering for the automatic identification of anxiety level subgroups. Using the CC mouse population model with diverse and reproducible genetic backgrounds, we identified 264 SNPs corresponding to 141 known that are significantly associated with anxiety-like behavior. A number of genes identified in our study have been linked to anxiety-like phenotypes in transgenic animal models and confirm that population-based studies using CC mice are a powerful and unbiased approach to identify candidate genes associated with specific phenotypes. For example, Ntrk3, Tnr, Cacna1h, Clstn2, and Rapgef2 were found to be involved in pathological processes of anxiety. The mRNA expression of Neurotrophic Receptor Tyrosine Kinase 3 (Ntrk3) decreased in central amygdala nucleus of young primates with high anxious temperament 46 . Overexpression of NTRK3's endogenous ligand (Ntf3) in the dorsal amygdala resulted in reduced anxious temperament and altered function in the anxious temperament neural circuit, implicating the role of neurotrophin-3/NTRK3 signaling in mediating primate anxiety 47 . Tnr knockout mice displayed decreased motivation to explore and an increased anxiety, which was more easily influenced by environmental factors 38 . Anxiety analysis in the BXD recombinant inbred mouse population also identified Tnr and subsequent systems genetic analysis showed that Tnr was co-expressed with genes related to psychiatric disorders 48 . Overexpression of Cacna1h induced anxiety and genetic ablation of the Cacna1h gene results in an anxiety-like phenotype in mice, suggesting normal Cac-na1h state is crucial to anxiety and both activation and inhibition of these channels in stressful condition may produce anxiety 35 . Clstn2 knockout mice displayed high exploration and hyperactivity affecting anxiety parameters 37 . Rapgef2 knockout mice exhibited hyperlocomotion phenotypes and decreased anxiety-like behavior 36 .
In addition to genes already associated with anxiety, our study identified a number of genes not previously associated with anxiety. A role for Abl2, Csmd2, Dlgap1, and Isl1 in neurodevelopment and psychiatric and neurodegenerative diseases have been reported in previous studies and our association analysis suggests a role in anxiety. Activation of Abl2/Arg kinase can alleviate corticosteroid-induced dendrite loss and behavioral deficiencies whereas Arg knockout mice exhibit synapse and dendrite loss and behavioral deficiencies 40 . Knock-down of Csmd2 results in reduced filopodia density in immature developing neurons and reduced dendritic spine density and dendrite complexity, implicating its association with certain psychiatric disorders 42 . Dlgap1 knockout mice exhibit post-synaptic density (PSD) disruption and reduced sociability, consistent with reports of Dlgap1 variants in schizophrenia and autism spectrum disorder (ASD) 41 . Conditional deletion of Isl1 using a Six3-cre transgene results in an early and persistent defect in cholinergic neuron differentiation and these dysfunctions have been implicated in various psychiatric and neurodegenerative diseases 39 . The comparison of our anxietyrelated genes with human GWAS for neurological conditions identified Stag1 and Sorcs3 associated with four conditions including "Alzheimer's disease", "autism spectrum disorder, "depression", "feeling nervous", "feeling worry", and "neuroticism". The roles of Stag1 and Sorcs3 in anxiety were not reported before, but the human GWAS data and our mice GWAS data imply that there could be a connection between these genes and anxiety.
The gut-brain axis plays important roles in neuropsychiatric disorders. Bidirectional interactions between the central nervous system and gut microbiota are maintained by different pathways including: direct activation of neuronal pathways, microbial metabolism of nutrients and production of circulating mediators, and immune activation and circulating inflammatory mediators. Our results showed that abundance of Ruminococcaceae was significantly higher in HA than in LA mice, and was positively correlated with the level of anxiety. Kang et al also showed that abundance of Ruminococcaceae correlated negatively with "percent of time in light" (i.e. lower Ruminococcaceae levels correlated with less anxiety) 49 . This study suggests that specific gut microbes could be used as a biomarker for anxiety or cognition and perhaps even targeted for therapy 49 . In a separate study, the levels of depression, anxiety, and eating disorder psychopathology at an inpatient admission were associated with the composition and diversity of the intestinal microbiota 50 . Significant changes in the composition of the intestinal microbiota were seen in patients with anorexia nervosa during re-nourishment, particularly among genera falling in the family Ruminococcaceae. Furthermore, in a study of psychological distress in patients with irritable bowel syndrome identified that patients with anxiety were characterized by elevated Bacteroidaceae 51 . However, not all literature is consistent in terms of abundance levels of Ruminococcaceae and anxiety. For example, anxiety and depression were associated with decreases in OTUs belong to the family Ruminococcaceae 52 . Also, abundance levels of Ruminococcaceae_UCG-014 correlated negatively with anxiety severity and positively with anxiety reduction 53 . To date, there appears to be no known mechanism regarding the role of Ruminococcaceae in anxiety, and further studies are needed to explore this. Studies in mice and humans have shown that the abundance level of pathogenic bacteria can increase anxiety-like behavior 54,55 , which could be due to their ability to produce exotoxins and promote conditions favoring inflammation. Our mediation analysis indicated that the effect of seven genetic loci could affect anxiety by altering abundance levels of Ruminococcaceae, Bacteroidaceae and Clostridiaceae. For example, our study suggests that Sorcs3 regulates anxiety by modulating abundance of Bacteroidaceae. Interestingly, genetic variants in Sorcs3 have previously been associated with abundance levels of Peptoniphilus in the nasopharynx 56 . These findings suggest that the association of SORCS3 with Alzheimer, depression, feeling nervous and neuroticism could be mediated by changes in the host microbiome.
The pathogenesis of anxiety disorders is complex, and involves intricate interactions between biological factors, environmental influences and psychological mechanisms. Even though we have found genetic and microbiome influences on anxiety disorders the genetic underpinnings of anxiety remains poorly understood. Together, our results suggest a complex genetic-microbiome interplay in the modulation of mouse anxiety. This study lays the foundation for future research to evaluate treatments for anxiety taking into account both host genome and microbiome. Light/dark test and mouse behavior characterization. We performed a light/dark (LD) behavioral test commonly used to study anxiety-like behavior in mice 59 . The light/dark test is based on rodents' natural aversion to bright areas and their spontaneous exploratory behavior. The device consists of two compartments each a polyvinylchloride box (20 cm × 40 cm × 40 cm) covered with plexiglas, one of the boxes is black and the other is transparent and illuminated. A small opening (7 cm ×7 cm) connects the two compartments. The light/white compartment was illuminated by a light from the ceiling (~350 lx at the floor of light/dark apparatus), while the black/dark compartment was not illuminated. Individual mice were placed into the light/white compartment, and then allowed to explore the enclosure freely for 5 minutes. All testing was conducted during the animal's light cycle. A video camera recorded the light compartment of the apparatus. Between consecutive tests, the instrument was cleaned with water. After the LD test, we constructed a computational pipeline to systematically evaluate mouse behavior following three steps: (1) mouse tracking, (2) behavior profiling and (3) feature extraction (Fig. 1A-E). The first two steps resulted in the mouse behavior profile, which characterizes the dynamic mouse size in the light chamber across time and provides a visualization tool for phenotype definition, extraction and quality control. Finally, seven anxiety-related phenotypes were quantified (Table 1).

Anxiety level characterization.
Consensus clustering (ConsensusClusterPlus package v1.50 in R) was performed using all mice in our study to determine behavior subgroups that correlate with different anxiety levels based on seven anxiety-related phenotypes, where 80% of the samples were bootstrapped 100 times. Euclidian distance was adopted for similarity measurement, and k-mean clustering was used as the clustering algorithm. After consensus clustering, the number of clusters was determined by the consistency of clusters, and the cumulative distribution information among different numbers of clusters. Finally, based on K=2 clusters, we assigned each mouse to either high anxiety (HA) or low anxiety (LA) groups (Table S1). Data from all mice was used for genetic association, microbiome association and mediation analysis.
Fecal sample collection and microbiome analysis. Fecal samples were collected from individual cages 16 h after a cage change as reported previously 14, 28 24 hours prior to the light/dark test. At least four independent cages were sampled for each CC strain. After collection, fecal samples were stored at -80°C until microbial analysis. We extracted genomic DNA from the homogenized fecal samples using the PowerSoil DNA Isolation Kit (http://www.mobio .com/) according to the manufacturer's instructions. PCR amplification of the V4 region of the 16S rRNA gene was performed using modern primers 60 . Amplicons were sequenced on an Illumina MiSeq using paired, 250 base-pair reads, according to the manufacturer's instructions and analyzed as previously described 30 . Animals were group housed whenever possible and individual mice were assigned the microbiome profile from their respective cage.
QTL analysis of anxiety. Genotype data of 134,593 SNPs was obtained from the UNC Systems Genetics Core website (http://csbio .unc.edu/CCsta tus/index .py). SNPs were filtered based on minor allele frequency ≥3 out of the 30 CC strains, leaving 70,273 SNPs. At each SNP, anxiety level (i.e., sub-group assignment) for all CC mice were assigned to their respective alleles. The Chi-square test was used to test the significance of associations between anxiety level (high anxiety or low anxiety based on consensus clustering) and allele classes at each SNP. SNPs with a p-value less than 1.00E-13 were selected. Putative candidate genes were defined as those genes (gencode.vM741) containing a significant SNP within the boundaries of the gene sequence. Gene Ontology biological annotations of putative genes were determined using ClueGO and visualized in Cytoscape. Human GWAS data was downloaded on 09/13/2019 from https ://www.ebi.ac.uk/gwas/. The significance of overlap between mouse and human candidate genes was calculated based on the hypergeometric distribution 61 .
Association between microbiome and anxiety levels. Mann-Whitney U Test (Matlab 2012b; Statistics Toolbox version 8.1) was utilized to identify microbiome features that have significantly different abundance (FDR < 0.1) between anxiety levels. Pre-selected microbiome features were correlated with anxiety levels using logistic regression (R 3.6.0; stats Package 3.6.1) (FDR < 0.01). To evaluate the effectiveness of significant microbiome features for the prediction of anxiety levels, random forest classification model (Matlab Random Forest Package version 4.6-14) was used with 100 cross-validation iterations, where, in each iteration, the number of trees was experimentally optimized to be 1000, 90% of samples were randomly selected for training and the rest 10% samples were used for testing.