Microbiota–Muscle/–Immune Interactions in Rhesus Macaque under Simulated Weightlessness Revealed by Integrated Multi-Omics Analysis

Long-term exposure to microgravity during spaceight has adverse effects on human health including muscle atrophy, impaired immune function, and alterations in gut microbiome prole. Gut microorganisms inuence a wide range of host biological processes, but their interactions with skeletal muscle and the immune system under microgravity are not known.

the reduced abundance of butyrate-producing colon bacteria Eubacterium, Roseburia and their crossfeeding bacteria Bi dobacteria may contribute to both the impaired immune function and muscle atrophy caused by HDBR.

Conclusions
We rst reported the HDBR-associated changes in gut microbiota composition, metabolomics of skeletal muscle and transcripts of PBMCs in non-human primate. Particularly, we revealed the underlying microbiota-muscle and microbiota-immune interactions during simulated microgravity, implicating that modulation of gut microbiota may represent a new strategy in enhancing crewmembers' health and safety during long-term space expeditions.

Background
The harsh environment of space poses a signi cant physiologic challenge to humans and has hindered the progress of deep space exploration. Astronauts experience a variety of pathophysiologic changes during space ight such as space motion sickness [1], spatial disorientation [2], osteoporosis [3], muscle atrophy [4], and impaired cardiac [5,6], cognitive [6,7], and immune [8] functions. The operation of the International Space Station (ISS) further extended our knowledge of the effects of microgravity on human health. For example, many astronauts on the ISS reported ocular issues including optic disc edema, globe attening, hyperopic shifts, and cotton-wool spots [9]. The National Aeronautics and Space Administration (NASA) twins study reported increased subfoveal choroidal thickness and peripapillary total retinal thickness [6]. Recently, a case of venous thromboembolism during space ight was also con rmed by radiologists on earth [10]. Generally, the physiological and functional alterations caused by space ight have been gradually well recognized.
Over the past two decades, a variety of high-throughput omics technologies such as methyl-CpG-binding domain protein sequencing [11], RNA sequencing [12], mass spectrometry (MS) analysis and metabolomics [13], have been used to investigate the impact of microgravity on human health. However, each omics technologies emphasize the role of molecules of their corresponding omic layer, but miss the complementary effects and interactions between omic layers [14]. Thus, an integrative approach to multiomics data analysis is needed to fully elucidate microgravity-induced physiological changes.
The gastrointestinal tract harbors complex communities of microbes that play important roles in maintaining human health, including energy extraction, vitamin biosynthesis, protection against pathogens, and development of the innate and adaptive immune systems [15,16]. Although changes in the gut microbiome pro le in response to microgravity have been reported in both human and rodents [6,17], how these relate to other microgravity-induced pathophysiologic alterations is unclear. To address this point, we carried out an integrative multi-omics analysis covering metagenomics of fecal samples, metabolomics of skeletal muscle and transcripts of peripheral blood mononuclear cells (PBMCs) from rhesus macaques subjected to -6° head-down tilted bed rest (HDBR). The multi-omics dataset generated in this study can serve as a resource for future investigations on the effects of microgravity on human health and may guide the development of potential countermeasures for future space ights.

Animal experiments
A total of 15 healthy male rhesus macaques aged 4-6 years and weighing 4-8 kg were provided by the Beijing Institute of Xie'erxin Biology Resource (Beijing, China). The animals were acclimatized for 3 months at the Laboratory Animal Center of China Astronaut Research and Training Center prior to being used in experiments. Ultimately, 5 fully acclimatized rhesus macaques were included in the study and subjected to 42 days of -6° HDBR as previously described [18]. After the treatment, the animals were individually housed in stainless steel mesh cages where they were allowed to recover for 32 days. The animals had free access to food and water and their general health condition was closely monitored for the duration of the study.

Muscle biopsy
Muscle tissue samples were obtained by biopsy as previously described [19]. Brie y, after general anesthesia with iso urane, the pre-HDBR muscle tissue sample (Pre-3) was obtained from the left soleus of rhesus macaques using an open biopsy technique. At the end of HDBR (H + 42), a biopsy was performed on the right soleus; and 32 days after HDBR, when the pre-HDBR biopsy site in the left soleus was fully healed, the muscle was biopsied at a different site (R + 32) (Additional le 1: Figure S1). All muscle samples were immediately frozen in liquid nitrogen after biopsy until use.

Transcriptome analysis
Total mRNA of PBMCs was ampli ed using oligo (dT) primers and sequenced by Complete Genomics (San Jose, CA, USA). At least 20 million reads were generated for each sample, and SOAPnuke was used to lter out those of low quality [20]. HISAT2 was used to map sequence reads to the reference genome and RNA-Seq by Expectation Maximization was used to calculate fragments per kilobase million (FPKM) [21]. Spearman correlations between all samples were calculated based on FPKM and were found to be very high (> 0.95) and biological replicates clustered well together. We used DESeq2 package in R to identify differentially expressed genes (DEGs), and performed Gene Ontology (GO) analysis for functional annotation [22].

Metabolome analysis
Metabolite detection was carried out by ultrahigh-performance liquid chromatography (UPLC)-MS. Quality control samples were included to evaluate data quality, and low-quality samples were removed. Feature alignment, picking, and identi cation were performed using Progenesis QI software (Nonlinear Dynamics, Newcastle, UK). MetaX software was used for data cleaning and statistical analysis [23]. Signi cantly altered features (P < 0.05, fold change < 1/1.2 or > 1.2, and variable importance in projection [VIP] > 1) were identi ed by combining uni-and multivariate analyses and annotated using Progenesis QI and the Human Metabolome Database (HMDB) v.3.6 and Kyoto Encyclopedia of Genes and Genomes (KEGG) database (www.genome.jp/kegg/) [24].
Clusters of co-abundant metabolites in muscle tissue were identi ed by weighted gene correlation network analysis (WGCNA) in R package [25]. A soft threshold of β = 9 for muscle features was selected by scale-free topology analysis for signed, weighted features co-abundance correlation network construction. A dynamic hybrid tree-cutting algorithm with deepSplit = 4 and a minimum cluster size of 30 was used for cluster identi cation. If the biweight mid-correlation between 2 clusters' eigenvectors was > 0.8, they were considered similar and were merged. Muscle feature clusters were labeled as M1-M6.

Metagenome analysis
Raw reads sequenced on the Illumina HiSeq 2000 platform (Expression Analysis, San Diego, CA, USA) at BGI were ltered to remove adapter contamination, low-quality reads, and host genomic DNA (rhesus macaque, assembly Mmul_8.0.1, National Center for Biotechnology Information [NCBI]). The remaining high-quality reads were assembled using metaSPAdes v.3.10.1 [26]. Open reading frames (ORFs) in contigs of each sample were determined using GeneMark v.2.7 [27]. The nonredundant gene set of all ORFs was clustered using CD-HIT v.4.5.7 based on the nucleotide sequence with thresholds of 95% identity and 90% coverage [28]. Taxonomic annotation of gene sets was carried out using CARMA3 based on BLASTP alignment with bacteria and archaea genes from the NCBI-NR database [29]. The gene set was also annotated against the KEGG v.59 database using BLAST v.2.2.23.
Gene abundance was calculated based on SOAP2 alignment and species abundance and functional pro les were summarized from their respective genes [30]. The differential alpha diversity, species, and KEGG orthologs (KOs) at different time points were calculated with the Kruskal-Wallis test using R software (https://www.r-project.org/). Spearman's correlation coe cient was used to evaluate the relationship between different genera and KOs or metabolic features according to their abundance.

HDBR-associated changes in gut microbiome pro le
To investigate the effect of HDBR on the gut microbiome, we carried out metagenome sequencing of 35 gut samples of rhesus macaque collected at 7 time points throughout the whole experiment (Additional le 1: Figure S1). A total of 286.86 Gb of raw data were generated; after removing low-quality and host reads, 275.99 Gb of high-quality data were obtained (average of 7.89 Gb per sample) (Additional le 2: Table S1). We constructed a reference library of the rhesus macaque gut metagenome using all of the samples (Additional le 3: Table S2); 63.10%, 21.00%, and 2.43% of the genes were annotated at the phylum, genus, and species levels, respectively, and 45.51% were annotated to 6631 KOs. Firmicutes and Bacteroidetes were the predominant phyla while Prevotella and Clostridium were the predominant genera (Fig. 1A, B). Genes related to metabolism constituted the largest proportion of total genes in the KEGG pathway analysis (Fig. 1C). 95.58% of KO functions and 49.48% of genera were shared between rhesus macaque and human gut gene catalogs (Additional le 4: Figure S2).
To identify microbial taxa affected by space ight, we compared pre-HDBR (control) and HDBR as well as HDBR and post-HDBR (recovery) gut microbiome pro les with the Kruskal-Wallis test. A number of bacterial taxa were signi cantly different between control and HDBR samples (1 class, 5 orders, 11 families, 55 genera, and 122 species), and some of these showed the reverse trend during recovery period (1 class, 1/5 orders, 3/11 families, 24/55 genera, and 72/122 species) (Additional le 5: Table S3). For example, the abundance of Acinetobacter and Lactococcus, two genera involved in the regulation of in ammation and protection against infections, was decreased by HDBR but returned to the baseline (pre-HDBR) level after 17 days of recovery ( Fig. 2A). However, the abundance of some genera was not restored even after 28 days of recovery; for example, the level of Bi dobacterium, a bene cial gut bacterium, decreased continuously throughout the experiment ( Fig. 2A).
We also identi ed 537 KOs with signi cantly altered abundance assigned to 139 genera, of which 27 were signi cantly altered throughout the experiment (P < 0.05). A correlation analysis showed that 44.76% of pairwise correlations between KOs and speci c genera were positive while 3.38% were negative (R 2 > 0.3, P < 0.05) (Fig. 2B). For example, a higher abundance of Myroides and Acinetobacter was associated with increased representation of K00121, K00151, K00276, K00451, and K01555, which are involved in tyrosine metabolism, indicating that changes in the abundance of speci c genera are associated with altered gut microbial function during HDBR.

Metabolic pro le associated with HDBR in skeletal muscle
To investigate the metabolic signature associated with HDBR in greater detail, we carried out a UPLC-MS-based metabolomic analysis to quantify metabolites in skeletal muscle tissue of rhesus macaque before, during, and after HDBR (Additional le 1: Figure S1 and Additional le 6: Table S4). We identi ed 356, 287, and 100 differentially abundant features (DAFs) in the HDBR vs control, HDBR vs recovery, and recovery vs control comparisons, respectively (P < 0.05, fold change > 1.2, VIP > 1) (Additional le 7: Figure  S3A). We focused on the 154 DAFs that were signi cantly altered in HDBR relative to the control and recovery samples (P < 0.05) (Additional le 7: Figure S3B). In the KEGG pathway analysis, the top 5 most affected pathways were tyrosine metabolism; biosynthesis of amino acids; protein digestion and absorption; alanine, aspartate, and glutamate metabolism; and tryptophan, consistent with accelerated protein degradation occurring in HDBR-induced muscle atrophy (Additional le 8: Table S5).
In addition, the interaction of these DAFs was evaluated by WGCNA. We classi ed 552 DAFs in muscle into 6 modules (Fig. 3A). DAFs involved in tyrosine metabolism (mcc00350) were mainly distributed into 3 modules (M1, M5, and M6) (Fig. 3B, C). DAFs in M1 showed co-abundance in muscle that the levels of all of these metabolites were decreased during HDBR but were reversed in the recovery group (Fig. 3D, E). Consistent with the known metabolic and functional changes in skeletal muscle under microgravity, the DAFs in M1 contained some critical amino acids involved in gluconeogenesis (eg, L-alanine [KEGG compound C00041, HMDB00161]) and fatty acid transport (eg, L-carnitine [C00318, HMDB00062]) (Fig. 3F). Notably, we identi ed epinephrine (C00788, HMDB00068) as the hub metabolite of M1 (Fig. 3G). Pyridoxamine (C00534, HMDB01431), a form of vitamin B6 that participates in free radical scavenging, was identi ed as the hub metabolite of M6 (Fig. 3G).

Transcriptional signature associated with HDBR in PBMCs
We collected PBMCs from 5 rhesus macaques throughout the experiment to evaluate the effect of HDBR on gene expression in immune cells (Additional le 1: Figure S1). A total of 779.09 Mb reads were obtained at 6 time points (Additional le 9: Table S6). By comparing every two time points (fold change > 2 and P < 0.05) and removing duplicate genes, we detected 65 DEGs (Additional le 10: Figure S4A) that were signi cantly enriched in 44 biological processes (P < 0.05); 41 were related to immune regulation including leukocyte activation (GO:0002694), T cell differentiation (GO:0045580), and interleukin (IL)-2 production (GO:0032663) (Additional le 10: Figure S4B and Additional le 11: Table S7). DEGs associated with these biological processes were mainly downregulated by HDBR, which is consistent with previous ndings that immune function was impaired by simulated or actual microgravity [8,31].

Integrated analysis of -omics data
To clarify the functional implications of microgravity-associated changes in gut microbiome pro le, we performed a multi-omic analysis of combined metagenomic, transcriptomic, and metabolomic data from fecal samples, PBMCs, and muscle tissue, respectively. As shown in Fig. 4 By integrating metagenomic data from fecal samples and metabolomic data from atrophied muscle, we found that 25 of the 27 differentially represented genera were closely correlated with 174 of 372 differentially expressed metabolites in muscle (Fig. 5), with Oligella, Sporosarcina, Citrobacter, Weissella, and Myroide showing the highest correlations. Notably, the abundance of Oligella was positively correlated with up to 45 metabolites and negatively correlated with 12 metabolites; and that of Myroide was negatively correlated with leucodopachrome (C05604, HMDB04067) and dopaquinone (C00822, HMDB01229), which are involved in tyrosine metabolism (mcc00350). Some bacterial genera were closely correlated with crucial metabolites identi ed in HDBR-induced muscle atrophy. For example, the abundance of Oligella was positively correlated whereas that of Citrobacter was negatively correlated with l-alanine; Lactococcus showed a positive correlation with L-carnitine; and Providencia was positively correlated with p-cresol (C01468, HMDB01858), a metabolite of tyrosine (Fig. 5).
Microbiota-derived short-chain fatty acids (SCFAs) are thought to mediate interactions between gut bacteria and other tissues [32,33]. For example, butyrate-a major SCFA produced by these microorganisms-prevents excessive in ammation by stimulating the function of M2 macrophages and regulatory T cells and inhibiting neutrophil in ltration [34]; moreover, administration of butyrate was shown to increase muscle mass and cross-sectional area in aged mice [35]. The abundance of butyrateproducing colon bacteria such as Eubacterium and Roseburia spp. and their cross-feeding bacteria Bi dobacterium was reduced during HDBR, suggesting that a lower level of butyrate may contribute at least in part to the impaired immune function and muscle atrophy caused by HDBR (Additional le 12: Figure S5). 3-Hydroxyphenylacetate can be transformed into 3,4-dihydroxybenzeneacetic acid (3,4DPHAA; C01161, HMDB01336)-a metabolite involved in tyrosine metabolism in skeletal muscle-by the enzyme 4hydroxyphenylacetate 3-monooxygenase (EC: 1.14.14.9), which is produced by Providencia rettgeri. 3,4DPHAA level was increased during HDBR and returned to the baseline level during the recovery phase, re ecting the changes in P. rettgeri abundance in the gut and suggesting that this bacterium in uences tyrosine metabolism in HDBR-induced muscle atrophy via 3,4DPHAA (Additional le 12: Figure S5). Besides, tryptophan 2,3-dioxygenase (EC: 1.13.11.11) produces N′-formylkynurenine (C02700, HMDB60485); which is transformed into formylanthranilic acid (C05653, HMDB04089) involved in tryptophan metabolism by kynureninase (EC: 3.7.1.3). Both tryptophan 2,3-dioxygenase and kynureninase are produced by gut microbes, such as Myroides and Comamonas, that showed signi cant changes in abundance during HDBR (Additional le 12: Figure S5). Collectively, these ndings indicate that changes in microbiota community composition in uence amino acid metabolism in skeletal muscle as well as immune function in rhesus macaque exposed to HDBR.

Discussion
It is well established that long-term space ight results in a broad spectrum of deleterious effects on human health. For example, the homeostasis of intestinal microbiota had been reported to be disturbed both in humans and rodents exposed to microgravity [6,17], however, the former human study involved only a single set of twins, while the ndings of the latter were not fully applicable to humans because of interspecies physiologic differences. Given the limited availability of biological samples from astronauts, more animal studies-especially in nonhuman primates-are needed to clarify the impact of microgravity on physiologic functions in mammals. In this study, we performed a longitudinal analysis of changes in the gut microbiome in response to HDBR. Consistent with the ndings from the NASA twins study [36,37], microgravity induced a decrease in the abundance of bene cial gut bacteria such as Bi dobacterium and an increase in that of opportunistic pathogens such as Escherichia coli; the lower abundance of Bi dobacterium during HDBR was also observed in fecal samples and continued to decrease during the 28-day recovery period. However, HDBR does not fully simulate the environmental conditions encountered during space ight and has relatively modest effects on the gut microbiome. HDBR caused no changes at the phylum level, and only 1 class showed differential abundance in fecal samples of rhesus macaque. In contrast, a 1-year space ight altered the abundance of 3 phyla and 3 classes of bacteria in the fecal microbiome of astronauts. Surprisingly, HDBR in rhesus macaque produced more differentially abundant families (11 vs 8), genera (55 vs 13), and species (122 vs 36) than did microgravity in humans. Given that microbiota community composition is strongly in uenced by diverse factors such as diet, lifestyle, and medication use, we speculate that the speci c diet and nutritional supplementation for astronauts may alleviate the negative effects of microgravity on gut microbiome composition, especially from family level to species level.
In addition to nutrient absorption, intestinal microorganisms maintain human health through bidirectional interactions with host biological processes including the immune system [38,39]. Studies in pathogenfree or gnotobiotic animals have demonstrated that the intestinal immune system in uences the compartmentalization of commensal microbiota and microbial community composition; reciprocally, gut bacteria affect the development of organized lymphoid structures and intestinal and systemic immune cell function [38,39]. In this study, HDBR reduced the abundance of butyrate-producing colon bacteria including Eubacterium, Roseburia, and Bi dobacterium. The SCFA butyrate is believed to mediate communication between commensal bacteria and the immune system [33]. It suppresses the production of proin ammatory cytokines released by M1 macrophages and neutrophils while activating Tregs, stimulating the production of the anti-in ammatory cytokine IL-10, and increasing the expression of interferon (IFN)-γ and granzyme B in cytotoxic T lymphocytes and Tc17 cells [34,40]. As simulated or actual microgravity can weaken immune function and increase susceptibility to infection as a result of reduced lymphocyte proliferation and IFN-γ production [8], the reduced abundance of butyrate-producing microbiota may contribute to the immune dysfunction caused by HDBR. Besides Bi dobacterium, we also found that Klebsiella and Kluyvera were closely associated with immune dysfunction. Klebsiella species causes infections such as pneumonia, urinary tract infections (UTIs), bloodstream infections, and sepsis [41]; and Kluyvera is a potentially pathogen that infects hosts under various conditions [42]. It remains to be determined how changes in the abundance of these two genera in uence host susceptibility to infection.
The relationship between intestinal microbiota community composition and skeletal muscle has recently been investigated [43,44]. Reduced muscle mass was shown to be closely correlated with a speci c gut microbiome signature [43][44][45]. For example, muscle atrophy in ghrelin-null mice was accompanied by selective depletion of butyrate-producing bacteria such as Clostridium XIVa and Roseburia [44], and there is evidence for an association between frailty and reduced abundance of butyrate-producing bacteria [46]. The administration of probiotics such as Lactobacillus reuteri and Faecalibacterium has been shown to alleviate the loss of muscle mass [47][48][49]. Butyrate treatment was also found to increase muscle mass and cross-sectional area in aged mice [35]. Based on these previous ndings, it is reasonable to come to the conclusion that muscle atrophy caused by HDBR can be partly attributed to the reduced abundance of butyrate-producing bacteria. Five other bacterial genera (Oligella, Sporosarcina, Citrobacter, Weissella, and Myroide) that have been mainly implicated in UTIs were closely related to abnormal metabolism of amino acids in HDBR-induced muscle atrophy, although the mechanistic basis for this association is unclear. We speculate that intestinal mucosal barrier disruption and enhanced intestinal permeability under microgravity results in the release of LPS and proin ammatory cytokines by these pathogens, which may promote muscle atrophy [36]. Supporting this possibility, bacteria-derived indoxyl sulfate was reported to enhance the expression of myostatin and atrogin-1 in atrophic skeletal muscle [50].

Conclusions
In summary, we described microbiota-immune/-muscle interactions in rhesus macaque under simulated microgravity using a multi-omics approach. Despite the phylogenetic proximity, our results may not be completely recapitulated in human. Nevertheless, our research has extended our understanding of the effects of gut microbiota on host physiology under microgravity and provided us with a new strategy to enhance astronauts' adaptation in space, for instance, increasing the abundance of butyrate-producing microorganisms in the intestine during space missions through dietary supplementation.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.