Changes in gut viral and bacterial species correlate with altered 1,2- 2 diacylglyceride levels and structure in the prefrontal cortex in a non-human 3 primate model of depression.

31 Background: Major depressive disorder (MDD) is a debilitating mental disease, 32 but its underlying molecular mechanisms remain obscure. Gut microbiome can 33 modulate brain function and behaviors through the microbiota-gut–brain (MGB) 34 axis in depression. Our previously established non-human primate model of 35 naturally occurring depression-like behaviors, which is characterized by MGB 36 axis disturbances, can be used to interrogate how a disturbed gut ecosystem 37 may modulate the MDD onset. To better clarify the molecular interrelationships 38 and downstream functional consequences in the MGB axis on MDD pathology, 39 here, gut metagenomics were used to characterize how gut virus and bacterial 40 species, and associated metabolites, change in depressive monkey model. 41 Results: We identified a panel of 33 gut virus and 14 bacterial species that 42 could discriminate the depression-like (DL) from control M. fascicularis . In 43 addition, using lipidomic analyses of central and peripheral samples obtained 44 from these animals, we found that the DL macaque were characterized by 45 alterations in the relative abundance, carbon-chain length, and unsaturation 46 degree of 1,2-diacylglyceride (DG) in the prefrontal cortex (PFC), in a brain 47 region-specific manner. In addition, lipid-reaction analysis identified more active 48 and inactive lipid pathways in PFC than in amygdala or hippocampus, with DG 49 being a key nodal player in these lipid pathways. Significantly, co-occurrence 50 network analysis showed that altered gut viral and bacterial species, and their 51 interaction may be relevant to the onset of negative emotions behaviors by 52 modulating the DG levels in PFC in the depressive macaque. 53 Conclusions: Our findings suggest that altered gut virus and bacteria as well 54 as DG levels and structure in the PFC are hallmarks of the DL macaque, thus 55 providing a new framework for understanding the gut microbiome's role in 56 depression. 57

7 associated with the reads and their mated reads were removed. Metagenomics 166 data were assembled using MEGAHIT [28](https://github.com/voutcn/megahit), 167 contigs with the length being or over 300 bp were selected as the final 168 assembling result. All predicted genes with a 95 % sequence identity (90% 169 coverage) were clustered using CD-HIT [29] (http://www.bioinformatics.org/cd-170 hit/). Reads after quality control were mapped to the representative sequences 171 with 95% identity using SOAPaligner [30] (http://soap.genomics.org.cn/). Here, 172 prior to construction of gene sets, we initially removed the host sequences with 173 high homology to the macaca genome from microbial metagenomes. Based on 174 the NCBI NR database, we annotated gene sets for bacteria, fungi, viruses, 175 protozoa, and archaea using Diamond (Version 0.8.35). Based on a unified 176 database, each gene is assigned to the highest scoring taxonomy, 177 which facilitates simultaneous assessment of these microbial species in the gut 178 ecosystem of depressive macaca [31]. The differential bacterial species and gut 179 viruses between the two groups were identified using Linear discriminant 180 analysis (LDA) Effect Size (LEfSe) with an LDA score >2.0 [32]. 181

Lipid identification and analysis 182
Details of the lipidomic methods were based on LC-MS/MS, similar to our 183 previously published protocols [8,26]. Briefly, brain tissue samples were 184 in the pathway, as follows: 221 As shown in references [34,35], also follows a normal distribution. To 223 determine if pathway A was significantly more active, we again chose the 224 significance level (P) to be 0.05, corresponding to a > 1.645 for reaction to 225 be considered as significantly active. Again, these Z scores were multiplied with 226 -1 for the reactions which were significantly more active in HC as opposed to 227 DL. For visualization, a negative Z-score indicates inactive pathway in DL 228 relative to HC, and a positive Z-score indicates active pathway in DL relative to 229

Statistical analyses 231
Statistical analyses were carried out using SPSS (version 21) and R studio 232 (version 3.5.2, 2018), unless otherwise described. The discriminating lipid 233 species between the two groups were considered significant at the threshold of 234 VIP > 1.0 and P-values < 0.05. The levels of all lipid subclasses were analyzed 235 by unpaired two-tailed Student's t-tests. Error bars represent the standard error 236 of the mean (SEM) in all cases. LEfSe analysis was used to identify the different 237 gut bacteria and viruses by estimating the effect of the abundance of each 238 species (LDA > 2 and P value < 0.05). In established lipid pathways, the 239 reactions were considered as significantly active (or inactive) for their calculated 240 Z score > 1.645 (or < -1.645), corresponding to the P-value < 0.05. Correlation 241 data analysis was performed using the spearman's regression model with the 242 significance threshold of P < 0.05. The investigators were not blinded to the 243 group classification while analyzing the data. 244

Gut virus differences between DL and HC macaques 246
The obtained metagenomic sequencing was mapped to the known viral 247 genomes NCBI NR database. From the filtered sequence data, we annotated 248 1311 viral species based on the NR database. Using LEfSe analysis, we 249 identified 33 differential gut viral species that could discriminate between the 250 DL and HC groups (Figure 1a,     Brain function relies on the homeostasis of cellular membranes, and its 299 perturbation might partially explain the neuronal deficits found in depression. 300 The carbon chain length and the degree of unsaturation of fatty-acyls are two 301 key elements of lipid architecture that impact membranes' functional 302 biophysical properties. Our lipidomic analyses demonstrated not only 303 differences in relative lipid composition and levels between HC and DL 304 macaques, but also differences in lipids' biochemical structures with potential 305 structure-function consequences on neural cell activity. Among the altered lipid 306 subclasses, DG showed profound structural alterations, with multi-level 307 alterations of carbon chain length and degree of unsaturation. Compared with 308 HC, the DL macaques had altered levels of DGs with carbon chain lengths of 309 34C, 36C, and 43C in PFC, and 38C in AMY and HIP (Figure 3d). The 310 differential carbon chain length DGs were primarily in the range of 34C to 40C, 311 a range which comprised the majority of DG species. We also analyzed the 312 degree of unsaturation for each subclass and DG. Compared with HCs, the 313 degree of unsaturation of DG was significantly increased in 2, 3, and 10 double 314 bones in PFC, but there was no difference in unsaturation degree of DG in AMY 315 or HIP (Figure 3e). Taken together, DG structure and degree of unsaturation 316 was the most different in the PFC, in HC versus DL macaques. 317 In the periphery (plasma), the structural changes in DG species were 318 distinctly different than seen in brain. In plasma, the typical DG carbon chain 319 length was longer, mainly 62C to 67C (Figure 3d). The degree of unsaturation 320 was lower, mainly 0 to 3 (Figure 3e), compared with 2 to 6 carbon bones of DG 321 in brain. There were no difference in carbon chain length or unsaturation degree 322 in plasma DG. 323 Given the altered abundance and structure of brain DG, next we employed 325 lipid-reaction analyses to explore whether specific lipid pathways were 326 dysregulated in key depression-relevant brain regions in the DL macaque. activity was significantly changed across more than one brain region examined 341 (Figure 4b, 4c). This most dysregulation of lipid pathway activity occurred in 342 pathways and found some interesting lipid shape changes in these altered lipid 354 pathways in the PFC of the DL macaques. Based on the size ratio of the 355 hydrophilic head to the hydrophobic tail, glycerophospholipids and glycerolipids 356 can be classified into cone (small head, two or more hydrophobic tails; including 357 DG, PA, and TG), cylinder (large head and two hydrophobic tails; including PC 358 and PI), or inverted cone (large head and single hydrophobic tail; including LPI 359 and LPC) shapes. As shown in Figure 4f, the altered lipid pathways tend to 360 promote (red arrow) lipid shape evolution from inverted cone lipids to cylinder 361 and cone lipids, and to inhibit (blue arrows) the opposite evolution. DG lay at 362 the intersection of these evolutions of lipid molecular geometry, and was 363 involved in multiple processes. 364

DG levels in DL macaques versus HCs. 366
To explore the potential interactions of these microbiome and molecular 367 changes along the MGB axis of DL macaques, we constructed co-occurrence 368 networks of altered gut viruses, bacteria, and DG in the PFC of DL versus HC 369 macaques. Using an edge-weighted spring-embedded layout, the network was 370 visualized and the nodes were spontaneously mutually attractive or exclusive 371 based on the coefficient between nodes (Figure 5). Overall, co-occurrence 372 analysis showed that gut viruses and bacteria formed strong and broad co- as well as a potential factor to inhibit plasma acetate levels and intrarenal RAS also found some Streptococcus enriched in healthy controls, suggesting that 435 gut Streptococci may play protective roles in the MGB axis. 436 Brain is particularly enriched in lipids, with a diverse lipid composition 437 compared to other tissues [44]. Changes in the composition and structure of 438 lipids in the brain profoundly affect neurodevelopment and signal transduction 439 in perception and emotional behavior, which may lead to depression and Nonetheless, there are some limitations of our study, which can also 521 provide direction for future research. First, due to the low reproductive rate and 522 morbidity, as well as ethical considerations, the sample sizes were relatively 523 small, thus the reliability of the association reported may be impacted. Second, 524 the effects on key viruses, bacteria, and brain lipids that we reported in this 525 study need further longitudinal independent validations with larger samples. 526 The key viruses and bacteria need further isolation and culturing from macaque 527 fecal samples, and more independent verifications in multiple animal models, 528 such as fecal microbiota transplantation (FMT) or microbial agents. Third, 529 depression as a complex mental disorder, and depression pathology is not 530 isolated to only the PFC, HIP, and AMY regions studied here (although those 531 are key regions of depression neuropathology). Thus further studies with 532 broader brain regions and technologies such as FMRI, calcium imaging, and 533 optogenetics are needed. Fourth, MGB crosstalk in depression involves 534 multiple mechanisms and metabolic pathways beyond those which we focused 535 on in this study. The vagus nerve, hypothalamic-pituitary-adrenal axis, and 536 neuroimmune mechanisms, are all worthy of further study. 537

Conclusions 538
Taken together, using metagenomic data, we found the altered gut virome, 539 especially bacteriophages, plays a role in the onset of depression. Through 540 multiomics approaches, we have presented evidence that DL macaques were 541 characterized by disturbances of gut-virus, bacteria, and DGs in the PFC.

646
Analysis of lipid pathway activity showed altered predicted lipid fluxes occurred in brain regions.

647
Calculated Z-scores were chosen to indicate the active/inactive pathways in DL relative to HC,

Supplementary information 678
Addition file 1: Figure S1. Discriminating lipids in brain regions and plasma between HC and 679 DL groups.

680
Addition file 2: Figure S2. Pathway activity of lipid showed changed reactions among brain 681 regions in DL group relative to HC group.

682
Addition file 3: Table S1. Discriminatory gut virus between DL and HC groups.