An integrated metagenomics and metabolomics approach implicates the microbiota-gut-brain axis in the pathogenesis of Huntington's disease

BACKGROUND
Huntington's disease (HD) is an autosomal dominant neurodegenerative disorder with onset and severity of symptoms influenced by various environmental factors. Recent discoveries have highlighted the importance of the gastrointestinal microbiome in mediating the gut-brain-axis bidirectional communication via circulating factors. Using shotgun sequencing, we investigated the gut microbiome composition in the R6/1 transgenic mouse model of HD from 4 to 12 weeks of age (early adolescent through to adult stages). Targeted metabolomics was also performed on the blood plasma of these mice (n = 9 per group) at 12 weeks of age to investigate potential effects of gut dysbiosis on the plasma metabolome profile.


RESULTS
Modelled time profiles of each species, KEGG Orthologs and bacterial genes, revealed heightened volatility in the R6/1 mice, indicating potential early effects of HD mutation in the gut. In addition to gut dysbiosis in R6/1 mice at 12 weeks of age, gut microbiome function was perturbed. In particular, the butanoate metabolism pathway was elevated, suggesting increased production of the protective SCFA, butyrate, in the gut. No significant alterations were found in the plasma butyrate and propionate levels in the R6/1 mice at 12 weeks of age. The statistical integration of the metagenomics and metabolomics unraveled several Bacteroides species that were negatively correlated with ATP and pipecolic acid in the plasma.


CONCLUSIONS
The present study revealed the instability of the HD gut microbiome during the pre-motor symptomatic stage of the disease which may have dire consequences on the host's health. Perturbation of the HD gut microbiome function prior to significant cognitive and motor dysfunction suggest the potential role of the gut in modulating the pathogenesis of HD, potentially via specific altered plasma metabolites which mediate gut-brain signaling.


Introduction
Huntington's disease (HD) is a progressive neurodegenerative disorder caused by a trinucleotide expansion (CAG) in the huntingtin gene, which is expressed ubiquitously throughout the brain and peripheral tissues (Moffitt et al., 2009). Despite being a genetic disease, the onset and severity of the symptoms is affected by a combination of environmental factors including diet, physical activity and stress (Gubert et al., 2020). There is increasing evidence implicating the gut microbiome as a factor influencing neurodegenerative diseases including Alzheimer's disease and Parkinson's disease (Sampson et al., 2016;Scheperjans et al., 2015;Vogt et al., 2017). Recently, two independent studies in mouse models and a study of clinical patients have shown evidence of gut dysbiosis in HD (Kong et al., 2018;Radulescu et al., 2019;Wasser et al., 2020). This was not entirely surprising, as gut dysbiosis is tightly associated with gastrointestinal motility and HD patients have bowel complications (Kobal et al., 2018). However, compositional differences in the gut microbiota may not have significant impact on its metabolic capacity due to the overlap in the functional capacity of different taxa (Louca et al., 2018). As these studies were performed using 16S rRNA amplicon sequencing, the functional implications of the HD gut dysbiosis could only be predicted.
One important function of the gut microbiome is the production of metabolites for the maintenance of the host's metabolic homeostasis as well as the shaping of the gut-brain-axis (Forsythe et al., 2014). A significant part of the metabolites circulating in mammalian blood are derived from the intestinal microbial community, some of which are important in regulating brain health (Matsumoto et al., 2013;Wikoff et al., 2009). In particular, short-chain fatty acids (SCFA) are notable microbial-derived metabolites from the fermentation of dietary fiber in the gut (Cummings et al., 1987). Manipulation of SCFA levels in the plasma via oral ingestion have shown improvements in the severity of symptoms in various metabolic and neurodegenerative diseases (Kidd and Schneider, 2010;Laurent et al., 2013;Kilgore et al., 2010).
Additionally, a stable configuration exists for an individual's gut microbiome which is critical for optimal health. This healthy, stable configuration consists of a highly personalized suite of gut microbiota or functional gene profiles that is consistently present during adulthood and naturally undergoes substantial changes during the extreme ends of life (Faith et al., 2013;Caporaso et al., 2011;Costello et al., 2009;Claesson et al., 2011;Yatsunenko et al., 2012;Turnbaugh et al., 2009;Qin et al., 2010). Environmental factors such as diet or antibiotic intake drive normal day-to-day variability in the gut microbiome (Lozupone et al., 2012). However, a stable gut microbiome typically restores its ecology and returns to a functional equilibrium state when these environmental pressures are mitigated (Lozupone et al., 2012). On the other hand, a volatile (unstable) gut microbiome lacks the ability to resist change in composition, restore its original composition, or recover it is initial function despite compositional changes in the face of disturbances (Lozupone et al., 2012;Moya and Ferrer, 2016). Failure to do so may result in an irreversible dysbiotic state, deviated from the normal healthy gut composition or function. The majority of previous studies on the stability/volatility of the microbiome had focused on healthy controls, but less is known about the effects of disease states (Faith et al., 2013;Caporaso et al., 2011;Costello et al., 2009;Moya and Ferrer, 2016). Whilst there are reports of increased volatility in patients with inflammatory diseases of the gastrointestinal tract, to our knowledge, there are no studies investigating the volatility of the gut microbiome in neurological diseases (Halfvarson et al., 2017;Martinez et al., 2008). Given that host genetics are one of the factors influencing gut microbiome composition, we hypothesized that the mutant HD protein in the gut could exert constant pressure on the gut ecology and result in a volatile gut microbiome.
Hence, in this study, we sought to provide holistic insights into HD gut dysbiosis from young to adult stages of life, prior to onset of overt motor deficits modelling clinical symptoms. The R6/1 transgenic mouse model of HD was used in this study. The R6/1 line is a mouse model with a slower progression of symptoms and better recapitulates the adult onset of HD, representing the majority of the cases, compared to the R6/ 2 line which is a better model for juvenile HD (Mangiarini et al., 1996). Previous findings of gut dysbiosis in the R6/1 line using 16S rRNA amplicon sequencing of fecal DNA was reported at 12 weeks of age (Kong et al., 2018). This present study aimed to extend those findings by determining whether the HD gut microbiome is significantly perturbed at earlier time points and characterizing the functional metagenomics profile in the HD gut. Additionally, we aimed to determine the impact of the HD gene mutation on the volatility of the gut microbiome. Finally, to characterize the plasma metabolome profile of these mice at 12 weeks of age, targeted metabolomics was performed. We aimed to unravel prospective relationships between the gut microbes and plasma metabolome which could be involved in the gut-microbiome-brain axis and potentially modulate other HD-related symptoms. In this study, the gut microbiome of the R6/1 mice was profiled at five timepoints: 4, 6, 8, 10 and 12 weeks of age, using shotgun sequencing and metagenomic analyses. Targeted metabolomics was performed on the plasma of these mice at 12 weeks of age, followed by multi-omics integration analysis on the metagenomics and plasma metabolome data, to determine any potential effects of gut dysbiosis on the plasma metabolome profile.
Overall, this exploratory study indicated that an alteration of gut microbial composition and function in pre-motor symptomatic R6/1 mice at 12 weeks of age but not at any time points before that. However, the HD gene mutation had more profound effects on the stable state of the gut microbiome throughout the course of the experiment, showing heightened volatility in the gut microbial species and functional capacity compared to their wild-type littermates. Analysis of plasma metabolomic profile revealed a signature of metabolites consisting of ATP and some muscle biomarkers that could distinguish the HD plasma from the controls. Integration of the 12-week timepoint bacterial species, functional metagenomics and the plasma metabolomics data provided some potential links between the gut microbiome and plasma metabolomics which can be validated in future studies.

Animal husbandry
Hemizygous R6/1 transgenic male mice were crossed with F1 (CBA x C57Bl/6) females to generate male wild-type (WT) and R6/1 littermates, herein termed HD mice (n = 9 per group). The mice were co-housed according to genotype due to the coprophagic nature of mice, and they were distributed among 4 cages per group, with 2-3 mice per cage. The animals had ad libitum access to water and chow in a temperate (22 • C) and humidity (45%) controlled room, with a 12-h light / dark cycle (lights on at 0700 h). Cage changes were performed weekly and body weight was monitored. All experiments were approved by the Florey Institute of Neuroscience and Mental Health Animal Ethics Committee and performed in accordance with research guidelines of the National Health and Medical Research Council.

Fecal samples and plasma collection
Fecal samples were collected at 4, 6, 8, 10 and 12 weeks of age as previously described (Kong et al., 2018). Briefly, individual mice were placed in clean cages for 5-10 min. Fresh pellets were collected and immediately frozen with liquid nitrogen prior to storage at − 80 • C until use. Blood was collected at 12 weeks of age via terminal cardiac puncture in EDTA tubes. The plasma was obtained after centrifugation at 5000g for 10 mins.

Fecal DNA extraction and shotgun sequencing
The 90 fecal samples were randomized prior to DNA extraction to ensure the samples were equally represented across sequencing runs so as to minimize batch effects. DNA was extracted from 50 to 200 mg of raw sample with preliminary step of bead beating using 0.1 diameter glass beads (BioSpec Products #11079101) on the Powerlyser 24 homogenizer (Mo-Bio #13155). Each sample was added to a bead tube filled with 750 μL of Bead Solution (Qiagen #12855-100-BS) and 60 μL of solution C1 was added to each sample tube and vortexed to mix. The tubes were heated at 65 • C for 10 min. The sample was then bead beat for 5 mins at 2000 rpm, and then centrifuged for 1 min at 10000g. The resulting lysate was transferred to a new collection tube. Extraction was as per Qiagen DNeasy Powersoil Kit (#12888-100). Final elution volume was 50 μL.
Libraries were prepared according to the manufacturer's protocol using Nextera XT Library Preparation Kit (Illumina #FC-131-1096). The only alterations to the protocol as outlined was the reduction of the total reaction volume for processing in 96 well plate format. Library preparation and bead clean-up was run on the Mantis Liquid Handler (Formulatrix) and Epmotion (Eppendorf #5075000301) automated platform. These programs cover "Tagment Genomic DNA" to "Amplify DNA" in the protocol (Mantis-Nextera XT library prep protocol) and "Clean Up Libraries" in the protocol (Epmotion -Library Clean Up protocol). On Completion of the library prep protocol, each library was quantified and QC was performed using the Quant-iT™ dsDNA HS Assay Kit (Invitrogen) and Agilent D1000 HS tapes (#5067-5582) on the TapeStation 4200 (Agilent #G2991AA) as per the manufacturer's protocol.
Nextera XT libraries were pooled at equimolar amounts of 2 nM per library to create a sequencing pool. The library pool was quantified in triplicates using the Qubit™ dsDNA HS Assay Kit (Invitrogen). Library QC is performed using the Agilent D1000 HS tapes (#5067-5582) on the TapeStation 4200 (Agilent #G2991AA) as per the manufacturer's protocol. The library was prepared for sequencing on the NextSeq500 (Illumina) using NextSeq500/550 High Output v2 2 × 150 bp paired end chemistry in the Australian Centre for Ecogenomics (University of Queensland, Australia) according to the manufacturer's protocol. The reads obtained were between 6 and 8 million reads for each sample.

Taxonomic assignment and microbial structure analysis
The reads were de-replicated and the host contaminated reads were removed for the subsequent analysis. The contaminant-free dataset has been made publicly available at NCBI under SRA accession PRJNA613182 and the R code for the following analysis have been uploaded onto https://github.com/gkong1/longitudinal_metagenomics _analysis. The quality of the reads was checked using FASTQC. To identify the taxonomic composition, the reads were classified using Kaiju v1.6.3 against the NCBI BLAST nr + euk database, under Greedy mode with MEM 11 parameter to maximize the classification rate (Menzel et al., 2016).
For alpha diversity analysis, reads were rarefied to 25,000 reads to calculate alpha diversity metrics including Observed and Inverse Simpson metrics using the Phyloseq v1.28 R package (McMurdie and Holmes, 2013). Two-way ANOVA was performed to determine the effect of genotype and age in any of these metrics while Tukey's HSD test was performed to determine the significance of the difference at each timepoint.
The data was normalized to their relative abundance and centered log ratio (CLR) transformation was applied to the data to mitigate the effect of compositional data before subsequent analysis. This includes analysis of cage and sequencing batch effects, time course analysis, differential abundance analysis, multivariate Partial Least Square Discriminant analysis (PLS-DA) and the sparse version (sPLS-DA). Only the features which passed the relative abundance threshold of 0.001% were included in the following analysis, resulting in a table of 2041 species. The normality of the CLR-transformed data was assessed before testing for differential abundance in the Phylum and Family level using two-way ANOVA.
To determine at which time point the data was discriminative of genotype, PLS-DA was performed on the CLR-transformed data at each time point. The significance was assessed using 999 permutations on the classification error rate (CER) obtained from leave-one-out cross-validation (LOO CV) method using MVA.test() from "RVAideMemoire" R Package, as described in (Westerhuis et al., 2008). To identify the microbial signature distinguishing the two genotypes, sPLS-DA which is implemented in mixOmics v6.8.5 R package, was performed (Le Cao et al., 2016). The default parameters were used to select a signature of 30 species and the LOO CV method was used to evaluate the performance of the method in terms of CER. As sPLS-DA is a supervised method, it is prone to overfitting which may lead to spurious results. Therefore, the model's performance needs to be assessed using crossvalidation to determine how the results could potentially be generalized to new, unseen data.

Metagenomic reads functional analysis
Metagenomic reads were annotated using MG-RAST Version v4.0.3 (https://www.mg-rast.org) (Glass et al., 2010).The Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used to explore the microbial functions in the samples. Similarity search between protein reads and the SEED/KEGG databases was conducted by using an E-value cutoff of 10-5 minimum identity cutoff of 60%, and minimum alignment length of 15aa. The annotated reads were sorted into Level 3 and Level 4 subsystems before applying a relative abundance threshold of 0.1%, resulting in a total of 333 genes and 245 KEGG orthologs (KOs) for the subsequent analysis. Similar to the taxonomic data analysis, PLS-DA was performed on the CLR-transformed data followed by CER (LOO CV) permutation tests to determine at which time point the data is discriminative of genotype, as described above. Variable selection was performed using sPLS-DA analysis, identifying a signature of 20 genes and 6 KOs which could distinguish between the two genotypes. The performance of the sPLS-DA model was assessed using the LOO CV method.

Identifying cage and sequencing run effects
To assess potential cage and sequencing run effects, we used linear mixed models on the CLR-transformed dataset to test for cage and run effects with cage or sequencing run as random effects (significance level 0.05). The PCA sample plots was also examined across all time points, and at each time point to assess potential age effects.

Time course analysis
We used the timeomics R pipeline from (Bodein et al., 2019) to perform the following analysis, focused on assessing relative stability or volatility over time. To cluster bacterial features, such as species/KO/ gene profiles across time, we modelled each feature across the samples of the same genotype using Linear Mixed Effect Model Splines (LMMS) (Straube et al., 2015). LMMS fits one of the best possible models for each species: a linear model which assumes that the response is not affected by individual variation (Model 1), a model of nonlinear response patterns (Model 2), a model fitted with subject-specific random intercept (Model 3) or a model fitted with both subject specific intercept and slope (Model 4). If little to no change is observed across time, the species is assigned to Model 0. This method accounts for between and within individual variability and irregular time sampling whilst avoiding underor over-smoothing.
Straight line modelling (Model 1) may occur when the interindividual variability is very high. An additional filtering step based on Breusch-Pagan test and mean square error were used to assess homoskedasticity of residuals to remove profiles that were considered noisy (Bodein et al., 2019). Temporal profiles of species/KO/gene profiles were clustered using PCA, as proposed by (Bodein et al., 2019). Each species/KO/gene makes a positive, negative or zero contribution to each component in the PCA and was assigned to clusters based on its maximum contribution to one of the components. Following clustering by PCA, sparse PCA was applied to obtain a smaller cluster signature.
As correlations can be spurious for compositional microbiome data, we used the proportionality distance between species identified as relevant in our method, as a measure of association (Bodein et al., 2019). A small value indicates that in proportion, the pair of profiles are strongly associated (0 represents the strongest). Median distances within clusters were calculated and compared to all the median distance between all species (Supp. Table 1). Closer number to 0 indicates stronger association between the members within each model.

Metabolomics sequencing 2.8.1. Quenching and extraction of polar metabolites
Polar metabolites were extracted in mouse plasma from different conditions using acetonitrile: methanol: milli-Q water (2,2:1 v/v/v, 180 μL) containing internal standards ( 13 C 10 , 15 N 5 -AMP, 13 C, 15 N-UMP, 13 C 6 -Sorbitol and 13 C 5 , 15 N-Valine) at 2 μM concentration per sample and vortexed. The homogenate was incubated for 10 mins at 4 • C on the thermomixer to ensure lysis of all cell membranes and then centrifuged to maximum speed 14,000 rpm at 4 • C to remove the cell debris. Around 150 μL of the clear supernatant was transferred into a HPLC vials with inserts for LC/MS analysis.
The mass spectrometry analysis was performed on an Agilent 6545B series quadrupole time-of-flight mass spectrometer (QTOF MS) (Agilent Technologies). The LC flow was directed to an electrospray ionisation source (ESI), where metabolite ionisation in negative mode was performed with a capillary voltage of 2500 V, a drying gas (N 2 ) pressure of 20 psi with a gas flow rate of 10.0 L/min, a gas temperature in the capillary of 300 • C and fragmentor and skimmer cap voltages of 125 V and 45 V, respectively. Data was collected in centroid mode with a mass range of 60-1200 m/z and an acquisition rate of 1.5 spectra/s in all-ion fragmentor (AIF) mode, which included three collision energies (0, 10 and 20 V). Prior to analysis, mass calibration was performed for negative mode to 0.5 ppm accuracy of the m/z value. Internal mass calibration was performed using Agilent ESI-TOF Reference Mass Solution containing purine (119.036320) and hexakis ( 1 H, 1 H, 3 H-tetrafluoropropoxy) phosphazine (981.99509) (API-TOF Reference Mass Solution Kit, Agilent technologies), which was continuously infused into the ESI source at a flow rate of 200 μL/min.
Targeted data matrices were generated using MassHunter Quantitative Analysis software (version B.09.00, Agilent Technologies) with metabolite identification (Metabolomics Standard Initiative (MSI) level 1) based on the retention time and molecular masses matching an authentic standard (Sumner et al., 2007).

Metabolomics analysis
The missing values (2 in total) were imputed and the metabolomic dataset was median-normalized. Similarly to the metagenomics analysis, PLS-DA was performed along with CER (LOO CV) permutation tests to assess the discriminative capacity of the metabolome data according to genotype. To select metabolite signature that could distinguish between the two groups, sPLS-DA was performed, obtaining a signature consisting of 20 metabolites, and the performance of the model was evaluated using LOO CV method.

Multi-omics data integration
We applied the DIABLO multi-omics integration method from the mixOmics v6.8.5 R package to integrate the shotgun sequencing data, both bacterial species and genes, and the plasma metabolomics data at 12 weeks of age, as well as to discriminate between the WT mice and HD littermates (Rohart et al., 2017). The CLR-transformed shotgun sequencing data and the median-normalized metabolomics data were used for this supervised integration method, allowing variable selection. A design matrix of 0.1 was used to focus primarily on the discrimination between the genotype groups leading to a selection of 30 species, 20 genes and 20 metabolites, as detailed in (Rohart et al., 2017). The resulting network was analysed further using Cytoscape v3.7.2 (Smoot et al., 2010).

SCFA extraction, derivatization and LCMS analysis
SCFAs were extracted and analysed by LCMS in accordance with the previously described protocol (Han et al., 2015). Briefly, 380 μL of acetonitrile:water (1:1, v/v) containing 4 μM of 4-methylvaleric acid as internal standard was added to the Eppendorf tubes holding 20 μL mouse plasma and the tubes were vortexed for 30 s before incubating for 10 min on a thermomixer (Eppendorf) at 10 • C with 950 rpm. Insoluble material was removed by centrifugation for 5 min at 4 • C at 16,200 g (Beckman Coulter Microfuge 22R refrigerated microcentrifuge) and the supernatant transferred to a fresh tube. SCFA standards were prepared from stock solutions (1 mM) in acetonitrile. Calibrators (0.1 to 50 μM) were made in acetonitrile:water (1:1, v/v) from serial dilutions of the standards. A small portion (40 μL) of the supernatant and calibrators was derivatized using 3-nitrophenylhydrazine (NPH) and 1-Ethyl-3-(3-dimethylaminopropyl) carbodiimide(EDC), which converts SCFA to their 3nitrophylhydrazones. The reaction was stopped using 20 μL of 200 mM quinic acid dissolved in acetonitrile:water (1,1, v/v). The samples were then spiked with 20 μL of 10 μM isotope-labelled internal standard mix prepared in accordance with the previously described protocol (Han et al., 2015). Finally, the samples (1 μL) were then injected into an Agilent 6490 triple quadruple mass spectrometer using the conditions previously described (Han et al., 2015).
Welch's t-test was performed to test whether the concentration of acetate, propionate and butyrate differed between WT and HD mice plasma. Significance threshold was set at 0.05.

Results
We present longitudinal fecal microbiome data for WT and HD mice from early adolescence (4 weeks) to the adult stage of life (12 weeks) to determine the onset of gut dysbiosis in HD. The 90 fecal samples yielded an average of 2 gigabases per sample (after trimming low-quality bases and removing sequences aligned to the mouse genome). A total of 20,928 species were detected, and those above a relative abundance of 0.001% were included for further analysis. The remaining 2041 species were not affected by cage or sequencing run effects.

Table 1
Number of features which are: "stable" in WT but "variable" in HD, "variable" in WT but "stable" in HD, "stable" in both WT and HD (Overlapping Stable), and "variable" in both WT and HD (Overlapping Variable). "Stable" group represents those features that are categorized into Model 0, whilst "variable" represents features categorized into the remaining the Models.

Subtle differences in bacterial species composition at 12 weeks of age
At the phylum level, the gut microbiome in both groups were predominantly members of the Bacteroidetes and Firmicutes phyla, followed by members of the Proteobacteria, Actinobacteria, and Verrucomicrobia phyla in much lower abundances (Fig. 1A-D). The phylum composition was examined at 5 different time points from the shotgun metagenomics sequencing data (Fig. 1A-D). No significant differences were observed in the any of the phyla when comparing HD with WT mice at any of the time points. At the family level, the Lachnospiraceae were the most abundant with an average abundance of 11.6%, followed by similar abundances of Bacteroidaceae (8.9%), Porphyromonadaceae (7.5%), Prevotellaceae (6.2%) and Clostridiaceae (2.4%) Similarly, no significant differences were observed in any of the bacterial families when Fig. 1. Overview of the gut microbiome in the HD mice and their WT littermates from 4 to 12 weeks of age. Boxplots of the relative abundances of the major phyla present in the mouse gut microbiome: (A) Bacteroidetes, (B) Firmicutes, (C) Proteobacteria and (D) Actinobacteria. No significant differences were observed in any of the phyla between the HD and WT mice at any of the timepoints measured. The alpha diversity indices, including the (E) number of bacterial species observed and (F) Inverse Simpson index, a measure of evenness, was accounted for from 4 to 12 weeks of age in these mice. There were significant effects of Age in both the Observed (p < 0.001) and Inverse Simpson index (p = 0.004) but only a significant genotype effect in the Inverse Simpson index (p = 0.017). A significant Age*Genotype interaction effect was found only in the Observed index (p = 0.019). Tukey's HSD post-hoc test showed no significant differences between the two genotypes at any timepoints in any of these measures. (n = 9 per genotype per time point).
comparing HD and WT mice at any of the time points (data not shown).
To investigate the microbial diversity within individual samples over time, the number of bacterial species observed (Observed) and Inverse Simpson index were computed. Two-way repeated ANOVA revealed no significant effects of age or genotype on the number of bacterial species observed (Fig. 1E). There was a significant effect of age (p = 0.006) but not genotype based on the Inverse Simpson index (Fig. 1F). Tukey's HSD post-hoc test showed no significant differences between the two genotypes at any timepoints in any of these measures.
To determine at which time point the data was discriminative of genotype, PLS-DA was applied to the dataset of each time point and permutation tests based on CER were performed to conclude on the significance of such discrimination. No significant effect was observed from Week-4 to Week-10 (p week4 = 0.352, p week6 = 0.708, p week8 = 0.701, p week10 = 0.428). The time point Week-12 was however, found to be significant (p = 0.017) to discriminate the two genotypes. Week-12 corresponds to the timepoint prior to onset of overt motor symptoms in HD mice. SPLS-DA was performed to obtain a signature consisting of 30 bacterial species that could distinguish between the HD and WT gut microbiome (CER = 0.27) at this time point.

The HD gut microbiome function is perturbed at 12 weeks of age
One advantage of whole genome shotgun sequencing is that it allows the profiling of microbial genes to interrogate the function of the gut microbiome, thus enabling the establishment of specific functional effects due to gut microbial composition alteration in the HD mouse. A total of 333 genes were above the cut-off relative abundance of 0.1%, Fig. 2. At time point 12 weeks, sPLS-DA analysis was performed on the KO pathways data to determine the signature capable of discriminating the HD mice from their WT littermates. The leave-one-out cross validation method shows a CER of 0.17 for the model, permutation test based on CER confirmed a significance of this signature (p = 0.01). The relative abundances of the top KO pathways selected are shown above. Decreased levels of galactose metabolism and benzoate degradation, as well as increased levels of sulfur metabolism, lysine degradation, glutathione metabolism and butanoate metabolism, were observed in the gut microbiome of HD mice compared to their WT littermates. (n = 9 per genotype). corresponding to 245 KEGG orthologs (KOs). No effects of cage and sequencing run were found on the bacterial genes and KOs, hence, the subsequent analysis were performed without further filtering.
For the KO analysis, the most abundant pathways were aminoacyl-tRNA biosynthesis with an average relative abundance of 6.7%, followed by alanine, aspartate and glutamate metabolism (5.2%) and ABC transporters and purine metabolism (4.6%). Permutation tests based on CER of PLS-DA revealed a significant discriminative effect only at Week-12 but not at any earlier time points, although a trend towards significance was observed at Week-8 (p week4 = 0.474, p week6 = 0.549, p week8 = 0.055, p week10 = 0.271, p week12 = 0.024). The KO pathways selected by sPLS-DA at Week-12 identified include galactose metabolism and benzoate degradation, both showing decreased levels in the HD gut microbiome compared to their WT littermates (CER 0.17, Fig. 2). Additionally, sulfur metabolism, lysine degradation, glutathione metabolism and butanoate metabolism were identified by sPLS-DA, showing increased levels in the HD gut microbiome compared to their WT littermates at 12 weeks of age (Fig. 2).

Modelling of bacterial species and functional trajectories revealed significant volatility in the HD gut
The gut microbiome is known for its stability, typically consisting of a set of 'core' or 'stable' entities with minimal variation over time within a healthy individual (Faith et al., 2013;Costello et al., 2009;Turnbaugh et al., 2009). As HD is a genetic disease and the expression of mutant HD protein in the gut could potentially affect the stability of gut microbial structure as well as the function from early adolescence (4 weeks) to the adult stage of life (12 weeks).
LMMS models categorized each species, genes or KO pathways into one of the four different models denoted as 0, indicating a relatively flat profile, 1 indicating either increasing or decreasing over time, 2 and 3 indicating complex curves across time (Fig. 4). Features that fit Model 0 are categorized as 'stable' while those that fit in the remaining models (Model 1, 2 and 3) are categorized as 'variable'. After filtering of the noisy profiles, models for 1708 species, 1228 genes and 124 KOs were able to be fitted for both HD and WT mice.
Out of the 1674 species, 1198 genes and 174 KO pathways which are 'stable' in the WT gut microbiome, a large proportion were 'variable' in the HD gut microbiome (Table 1). Namely, fatty acid metabolism, fatty acid biosynthesis, tryptophan metabolism and propanoate metabolism are among the stable pathways in WT that are variable in HD. The majority of the features that are stable in the HD gut microbiome were also stable in their WT counterparts. A small number of 'variable' features in WT remain 'stable' in the HD gut microbiome. The full lists of classification results of bacterial species, genes and KOs into 'stable' or 'variable' are provided in Supp. Table 2.

Targeted metabolomic profiling revealed no major differences in HD plasma 12 weeks of age
To further investigate the potential effects of gut dysbiosis in HD on the circulatory system, which is one of the routes of bidirectional communication between the brain and the gut microbiome, targeted metabolomics was performed. Due to the volume of blood required and the progressive nature of the gut dysbiosis in HD, only the plasma metabolites at 12 weeks of age in these mice were examined, following terminal cardiac blood collection. In total, the platform identified 221 metabolites, which belong to the following broad categories: amino acids, carbohydrates, lipids, nucleotides, peptides and xenobiotics. Although the PLS-DA CER permutation test revealed no significant discriminative effect on the plasma metabolites (p = 0.816), sPLS-DA identified 20 metabolites with some ability to distinguish the samples based on their genotype (CER 0.22, Supp. Fig. 2). The top metabolites selected included adenosine triphosphate (ATP), 3-Methylhistidine, urocanic acid, carnosine, threonic acid, homocitrulline, orotic acid, ADP, p-Aminobenzoic acid, 2-methylbutyrlglycine, trigonelline, alphahydroxyisobutyric acid, propionic acid, butyric acid, pipecolic acid, 2hydroxybutyric acid, isobutyrylglycine, 2-hydroxy-3-methylbutyric acid and ribitol (Supp. Fig. 2B).

Plasma SCFA were not significantly altered in HD mice
The gut function analysis revealed that butanoate metabolism was affected in the HD gut, hence, we sought to determine a possible change of SCFAs in the blood of the HD mice. Even though SCFA data was detected by the targeted metabolomics, the majority of the SCFAs in the circulatory system would have been metabolized by the liver and a In WT mice, the temporal profile of species, KOs and genes remained largely stable (Model 0: 98%, 97.6%, 100% respectively). In the HD mice, the majority of the temporal profiles of species, KOs and genes were categorized into Model 1 (Model 0: 23.5%, 41%, 39.5% respectively, Model 1: 74.9%, 57.3%, 58.1% respectively). (B) Concentration of short-chain fatty acid assay in WT and HD mice validated independently at 12 weeks of age. Welch's t-test revealed that the decrease in the concentration of propionate and butyrate in the HD plasma was not statistically significant when compared to their WT littermates. (n = 9 per genotype). minimal amount would have passed through to the bloodstream. Therefore, the targeted metabolomics approach may not accurately determine the concentration of these SCFAs. Hence, an independent SCFA assay was conducted to validate the SCFA levels in the plasma. Decreases in the plasma concentration of propionate (p = 0.11) and butyrate (p = 0.11) were observed, although statistically insignificant, in the HD mice when compared to their WT littermates. No differences in the concentration of acetate were observed (Fig. 4B).

Multi-omics integration revealed previously unknown relationships between the gut microbiome and plasma metabolites
Integrative analysis with DIABLO identified a signature of 30 bacterial species, 20 genes and metabolites that were highly associated (correlation >0.7 after variable selection, Fig. 5A). Further visualization in Cytoscape revealed that many microbes were found to be highly correlated with butyrate, homocitrulline, ATP, L-asparagine, 3-methylhistidine, orotic acid, and isobutyrylglycine.
The subsequent analysis was narrowed down to the metabolites with known functions in the nervous system, including ATP, butyrate and pipecolic acid, and the associated bacterial species which have been classified. Closer inspection of each network revealed many overlaps in the bacterial species between the three networks. Several Bacteroides species including B.pyogenes, B.ihuae, B.oleiciplenus and B.timonensis were positively correlated with butyrate but negatively correlated with ATP and pipecolic acid (Fig. 5B, Fig. 6). Similarly, Prevotella ruminicola and Odoribacter laneus positively correlated with butyrate levels and negatively correlated with ATP and pipecolic acid (Fig. 6).
Blautia producta was uniquely identified to be negatively correlated with butyrate. Prevotella scopos was identified to be negatively correlated with ATP levels and no associations were found with butyrate and pipecolic acid.

Discussion
The adult gut microbiome is noted for its ecological stability, consisting of a relatively unchanging 'core' which is resistant to changes in its composition, restore to its initial composition, or is able to recover its initial function despite compositional changes when facing disturbances due to external pressure (Caporaso et al., 2011;Claesson et al., 2011;Lozupone et al., 2012;Moya and Ferrer, 2016). This stability is considered critical for host health as it allows the gut microbiome to react to disturbances and failure to do so results in an irreversible dysbiotic state (Lozupone et al., 2012). Heightened volatility has been reported in inflammatory diseases of the gastrointestinal tract (Halfvarson et al., 2017;Martinez et al., 2008). To our knowledge, there are currently no studies investigating volatility of the gut microbiome in neurological diseases and the present study provides the first evidence of this kind. Previous studies investigating the stability/volatility of the gut microbial composition have adopted a Venn diagram-like approach in which taxa (or genes) that were consistently shared between individuals over time are considered as members of the 'core' microbiome (Faith et al., 2013;Caporaso et al., 2011;Costello et al., 2009). Another study has averaged both weighted and unweighted UniFrac distance across all individuals at each time point to estimate deviation of the diseased gut microbiome composition from the healthy state (Claesson et al., 2011). These aforementioned approaches in establishing a stable and unchanging core microbiome have some limitations. For example, the Venn diagram-like approach does not account for the abundances of individual taxa (or genes) at each time point. In addition, weighted distance-based analysis favors abundant features, thereby resulting in patterns driven by the most abundant features. We took an unbiased analytical approach in the current study by taking the temporal profile of each feature and categorizing into 'core/stable' features (whose abundance remains relatively unchanged throughout the course of experimental timepoint) or 'volatile/variable' features (whose abundance either increases and/or decreases over time).
Concurrent with previous findings, the WT gut microbiome has a stable microbial composition and function post-weaning (Schloss et al., 2012).The HD mice displayed a high intra-temporal volatility in the microbiome structure based on LMMS splines analysis, and even more strikingly, the functional attributes of the microbiome also fluctuated compared to the WT littermates. A large proportion of the species/ genes/KO pathways were 'stable' in WT but 'variable' in HD. Approximately a quarter of the bacterial species and about 40% of the microbiome function were 'stable' in both WT and HD mice, suggesting a set of microbial species and microbiome function that is unperturbed by the HD gene mutation in the gut. The factors which may confer stability of the gut microbiome includes diversity of bacterial species and functional response (Lozupone et al., 2012). Communities rich in diverse species are less susceptible to foreign colonization since the occupying species are specialized in utilizing potentially limiting resources more efficiently (Lozupone et al., 2012). However in the present study, there are no significant differences in the alpha diversity of gut microbial species or function (data not shown) between the of HD and WT littermates. Other possible explanations would include unexplored factors which are somehow affected by the HD gene mutation, resulting in the observed volatility. Nevertheless, the observed heightened volatility suggests that the HD gut microbiome lacks resilience to the disturbances of its stable state which may lead to altered gut-microbiome-brain axis during growth, resulting in widespread detrimental effects.
One of the aims of the current study was to extend the previous finding of gut dysbiosis in the R6/1 mouse model at 12 weeks of age, which utilized 16S rRNA amplicon sequencing of fecal DNA (Kong et al., 2018). This study aimed to extend those findings by determining whether the HD gut microbiome is perturbed at earlier time points as well as to characterize the functional metagenomics profile in the HD gut. The permutation test based on the CER of PLS-DA revealed a significant discriminative effect of genotype on the gut microbial composition and function at Week-12 only and not at any earlier time points. At this time point, the bacterial species signature selected by sPLS-DA that can distinguish the two genotype has a relatively high error rate, being able to correctly predict 66.6% of the samples. Given the progressive nature of the disease, subtle changes were anticipated at the onset of gut dysbiosis during the early pre-motor symptomatic stage of HD.
A perturbation in the functional capacity of the HD gut microbiome was observed at this time point as well. The sPLS-DA model showed higher accuracy when examining the genes (73%) and KO pathway (83%) signatures. Galactose metabolism and sulfur metabolism were the top KO pathways selected by the sPLS-DA. Increase in sulfur metabolism indicates an elevated hydrogen sulfide production in the HD gut which could adversely influence gut health and overall microbiome composition (Carbonero et al., 2012). Another notable pathway selected is butanoate metabolism, suggesting an increase in butyrate production in the HD gut microbiome at 12 weeks of age. In the gut, butyrate is the preferred source of energy by the colonocytes with a plethora of roles in promoting gut integrity through tight junction regulation, cell proliferation and stimulation of mucin production by goblet cells (Morrison and Preston, 2016) as well as anti-inflammatory effects (den Besten et al., 2013). However, elevated butyrate concentration in the gut lumen is cytotoxic and may induce severe epithelial cell apoptosis and disrupt the intestinal barrier, and is thereby toxic to the host (Peng et al., 2007). Additionally, pathogenic butyrate producers typically utilize lysine and other sources rather than pyruvate for butyrate production while releasing toxic byproducts into the gut lumen (Vital et al., 2014;Anand et al., 2016;Barker et al., 1982). The increase in lysine degradation coupled with elevated butanoate metabolism capacity indicates a bloom of pathogenic butyrate producers in the HD gut.
Notably, the HD blood plasma does not reflect this increase in butanoate metabolism in the gut. Validation of the plasma SCFA showed only a decreasing trend of butyrate concentration in the HD plasma at 12 weeks of age. Since excess butyrate not taken up by the colonocytes bacterial species (purple), genes (green), and metabolites (light blue), and lines indicate measure of association (correlation), either positive or negative. (B) Further visualization of the network from DIABLO using Cytoscape, focusing on butyrate and their relationship with classified bacterial species and genes which may contribute to their abundance. The size of the nodes indicate the number of interacting partners within the network. Many Bacteroides species, as well as Odoribacter laneus and Prevotella ruminicola, were found to correlate positively with butyrate. Blautia producta was found to correlate negatively with butyrate. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) will be taken up by the liver through the hepatic portal vein, it is possible that butyrate absorption by the liver or colonocytes is impaired, which would have system-wide effects throughout the host. (van der Beek et al., 2015). Nevertheless, changes in the HD gut microbiome function at the pre-motor symptomatic stage could have widespread downstream consequences in the modulation of other HD-related symptoms including alterations of metabolic homeostasis which are evident in HD (Liot et al., 2017). Future studies could focus on profiling the fecal metabolome as a functional readout of the gut microbiome to validate these exploratory findings (Zierer et al., 2018).
Additionally, targeted plasma metabolomics revealed no major alterations in HD plasma according to the CER based permutation tests. Regardless, there are some notable metabolites selected by sPLS-DA which may be useful in understanding HD pathogenesis. For example, the top metabolite marker selected was ATP, with an increase of over 4-fold in the HD plasma, and ADP increased to a lesser extent. This was opposite to expectations since mitochondrial dysfunction is prominent in HD, resulting in impaired ATP production, suggesting that the gut microbiome could contribute to the bulk of the changes found in the plasma (Cui et al., 2006;Jenkins et al., 1993). Indeed, bacteria could release ATP via mechanosensitive channels or stimulate ATP production by providing SCFA as food source for the colonocytes (Firmansyah et al., 1989;Donohoe et al., 2011). Additionally, the present metabolome data indicated decreased levels of carnosine and elevated levels of 3-methylhistidine (3MH) in the HD plasma, both biomarkers of muscle mass. 3MH is a metabolic product formed by methylation of histidine residues in the skeletal muscle and the only source of 3MH in the plasma is from the breakdown of these tissues (Asatoor and Armstrong, 1967;Long et al., 1975). On the other hand, carnosine is a dipeptide synthesized in the skeletal muscle and, to a lesser extent, in the brain and it is known for Fig. 6. Network analysis from DIABLO using Cytoscape, focusing on (A) ATP and (B) Pipecolic acid, and their relationship with classified bacterial species and genes which may contribute to their abundance. The size of the nodes indicate the number of interacting partners within the network. There were many overlaps between the members of the two networks shown. Many Bacteroides species, as well as Odoribacter laneus and Prevotella ruminicola, correlated negatively with ATP and pipecolic acid. In addition, Prevotella scopos also correlated negatively with ATP.
its role in exercise performance as well as brain function (Suzuki et al., 2002;Baguet et al., 2010;Schön et al., 2019). Reduced levels of carnosine were reported in individuals indiv with autism spectrum disorder or Alzheimer's disease, and dietary supplementation led to some improvements in clinical symptoms (Schön et al., 2019;Ming et al., 2012;Fonteh et al., 2007). The decrease in plasma carnosine levels may reflect diminished carnosine synthesis due to skeletal muscle atrophy, as indicated by 3MH levels, in these pre-motor symptomatic HD mice, concurrent with previous reports of HD-related muscle wasting (Zielonka et al., 2014). Future validation of these plasma metabolites would be required, and if confirmed, dietary supplementation of carnosine may be beneficial in improving HD symptomatology.
We sought to identify potential relationships between gut microbial and plasma metabolites involved in communication within the gutmicrobiome-brain axis . Hence, in the multi-omics integration, focus was placed on ATP and pipecolic acid, in addition to butyrate, since there is evidence supporting their role in CNS and they can be microbially derived (Matsumoto et al., 2013;Firmansyah et al., 1989;Donohoe et al., 2011;Ivanova et al., 2006;Takahama et al., 1986;Matsumoto et al., 2003;Edwards et al., 1992;Proietti et al., 2019). Matsumoto and colleagues have provided evidence for the contribution of the gut microbiome to the cerebral levels of pipecolic acid and ATP using germfree mice, however it is unclear which specific bacterial species contributed to those observations (Matsumoto et al., 2013). DIABLO analysis highlighted several key species overlapping between the three networks. Bacteroides species including B.pyogenes, B.timonensis, B.oleiciplenus as well as Odoribacter laneus and Prevotella ruminicola, were found to be highly correlated with all three metabolites. Odoribacter laneus is a common constituent of the human fecal microbiome although no functions, to date, have been tied to this particular species (Nagai et al., 2010). Prevotella ruminicola has known functions in SCFA production and was found to be reduced in our fecal samples (Strobel, 1992). On the other hand, negative correlations were found between Blautia producta and plasma butyrate levels. One report had shown that supplementation of Blautia producta together with other potential commensal species in germ-free rats increased the fecal concentration of butyrate (Becker et al., 2011). It is possible that the distal gut's functional output is not entirely reflected in the plasma despite a high correlation between the two datasets (correlation coefficient of 0.76). Nevertheless, this approach could help in untangling the complicated cross-talk between the host and the microbiome as well as identifying potential candidates for a probiotic intervention as part of future treatments for HD. Future studies would be required to determine the mechanistic link between the gut microbial species, its functional capacity and the plasma metabolites to better understand microbiome-gutbrain axis in HD in order to develop novel therapeutic interventions.
A recent study reported gut dysbiosis in clinical HD patients, comparing healthy controls with gene-positive patients, however no difference in the gut microbial structure was found when comparing between pre-symptomatic and symptomatic patients (Wasser et al., 2020). The current study revealed a small effect of gut dysbiosis at the pre-motor symptomatic stage which would be difficult to detect in a clinical setting when compared to animal models reared in a controlled environment, thereby minimizing the environmental effects on gut microbiome composition. The present study also has the limitation of having a small sample size, constrained due to the current costs of shotgun sequencing. Regardless, the preliminary findings here have confirmed and extended previous findings of gut dysbiosis in HD as well as uncovered a novel characteristic, which is the instability, of the HD gut microbiome. Given that power analyses are mainly restricted to parametric statistical methods, which were not applied in the present study to identify metagenomic signatures, future studies with larger sample size will enable to ascertain the validity and reproducibility of our findings.

Conclusions
The present study suggests the occurrence of gut dysbiosis, and more importantly, perturbation in the gut microbiome function in pre-motor symptomatic HD mice when compared to WT littermates. Additionally, heightened volatility in the HD gut microbiome indicates a lack of resilience to the disturbances of its stable state which may have widespread detrimental effects. The multi-omics integration analysis discovered novel relationships between the gut bacteria and plasma metabolome, suggesting the potential role of gut in modulating the pathogenesis of HD via specific altered plasma metabolites which mediate gut-brain signaling. The preliminary findings of the current study have generated new hypotheses to guide future directions in aiding the understanding of peripheral and central pathology in HD, and may pave the way for future development of novel therapeutic approaches.

Ethics approval
All experiments were approved by the Florey Institute of Neuroscience and Mental Health Animal Ethics Committee and performed in accordance with research guidelines of the National Health and Medical Research Council.

Funding
The shotgun sequencing was funded by a Computational Biology Research Initiative (CBRI) Grant to KALC and AJH, from the University of Melbourne. KALC is an NHMRC Career Development Fellow (GNT1159458). TR is an NHMRC-ARC Dementia Research Fellow. AJH is an NHMRC Principal Research Fellow. None of the funding bodies had any role in the design of the study or collection, analysis, and interpretation of data, or in writing the manuscript.