All surgical and postoperative care procedures were performed in accordance with the Animal Use Committee (AUC) guidelines of the Universidad Católica de la Santísima Concepción, Concepción, Chile. Twelve female and thirteen male C57BL/6 mice (8 weeks old) were purchased from the ISP (Chile) and randomly allocated to either the sham or DCM group. We divided 12 females into sham (n = 6) or DCM (n = 6) groups, whereas 13 males were divided into sham (n = 6) and DCM (n = 7) groups. The female and male DCM and sham groups were housed with 2–3 mice per cage. Fecal samples were randomly selected from 11 females and males for bacterial 16S rRNA analysis (Figure 1A) from different cages to avoid all sequenced mice from each group coming from the same cage. Tissue samples from 25 female and male mice were isolated and analyzed using flow cytometry (Figure 1A). The mice were anesthetized with 2% isoflurane, and a midline incision at the C4 level was performed. An aromatic polyether material was implanted underneath the C5–6 laminae to induce chronic and progressive compression of the cervical spinal cord as previously described (Vidal et al., 2017; Ulndreaj et al., 2022). A sham group was included for comparison, in which the aromatic polyether material was inserted under the C5-6 laminae for 30 s, and then removed. Animals within each cage received the same treatment (sham or DCM surgery). All animals received the same feeding diet (PROLAB RMH 3000) and water ad libitum, to avoid changes in microbiota composition owing to variations in feeding behavioral (Lee et al., 2022). All animals were housed in a facility with a 12:12 light/dark cycle.

Fecal samples were collected at 6, 9, and 12 weeks after DCM induction. The mice were removed from their home cage and placed in an empty cage with a cage filter top. Samples were collected over a 30-min period in a sterile tube, and immediately frozen, and stored at −80°C until further processing (Kigerl et al., 2016). Fecal samples were collected in the morning, to avoid confounding differences associated with circadian patterns. Fecal output was measured as previously described. Briefly, mice were individually placed into an empty cage, and fecal pellets were counted over a 15-min period (Sampson et al., 2016). Additionally, the gut length was measured at the end of the experiment.

The ladder walk consisted of two plexiglass plates (69 × 9 cm) connected by 10 rungs, each 1 cm apart as previously described (Ulndreaj et al., 2022). It was used to assess locomotor deficits in females and males. Two trials were video-recorded for each animal, and the foot fault steps out of all steps during each trial were counted. The average (of two trials) percent (%) foot fault steps were plotted for each animal (Ulndreaj et al., 2022).

Frozen fecal samples were shipped on dry ice to Microbiome Insights Inc. (Vancouver, Canada) for 16S rRNA sequencing. DNA extraction was performed using the MoBio PowerMag Soil DNA Isolation Bead Plate on a KingFisher robot. PCR-based amplicons were generated from the extracted DNA using primers targeting the V4 region of bacterial 16S SSU rRNA (Caporaso et al., 2012). Amplicons were further sequenced using an Illumina MiSeq machine in paired-end mode with 250 base-pair read length (Kozich et al., 2013). Thirty three female and male sequencing datasets were, respectively, processed using the nf-core/ampliseq framework (Straub et al., 2020), implemented on the NextFlow platform (Ewels et al., 2020). After read quality control and filtering using the Phred score, Amplicon Sequence Variants (ASVs) were obtained using the DADA2 pipeline (Callahan et al., 2016). Ribosomal ASVs were predicted using Barrnap.¹ To distinguish between known and novel taxa and remove unwanted taxa, filtered 16S rRNA sequences were extracted and classified with QIIME2 (Hall and Beiko, 2018) by a combination of a naïve-bayes algorithm and sequence alignment based on BLAST+ (Camacho et al., 2009) and VSEARCH (Rognes et al., 2016), respectively. The classified sequences were taxonomically assigned using the reference database Silva 138.1 prokaryotic SSU database (Pfeiffer et al., 2014). Finally, 744 high-quality ASV sequences were obtained, abundance tables were constructed, and the final results were obtained using MultiQC (Ewels et al., 2016). To obtain Bray-Curtis ordination plots, alpha-diversity estimations via Shannon entropy, and microbial abundance bar plots, respectively, DADA2 output tables were imported into R language and processed using the phyloseq R package (McMurdie and Holmes, 2013).

The 744 high-quality ASV sequences were aligned using the MAFFT multiple sequence alignment program v7.490 using the –reorder flag (Katoh et al., 2002, 2009). A phylogenetic tree was then constructed using FastTree software v2.1 using the MAFFT output alignments FASTA format as input along with the -nt flag. Edition and plot of the phylogenetic tree were generated using the Interactive Tree of Life server (iTOL), collapsing all clades whose average branch length distance was below 0.0002 (Letunic and Bork, 2007).

Differentially abundant ASVs between sham and DCM conditions were obtained by reading the DADA2_table.tsv into the R software and employing the R edgeRun package to interrogate ASV differentially abundant between sham and DCM conditions across sexes (Dimont et al., 2015). We considered significantly differentially abundant ASVs with a fold change over 0.5 and a False Discovery Rate (FDR) less than 0.05.

We aimed to match every differentially expressed ASV from the 16S rRNA sequencing experiment to a candidate microbial genome using the iMGMC mouse-gut dereplicated metagenome-associated genomes collection (iMGMC-hqMAGs-dereplicated_genomes.tar.gz) containing 132,958 contigs distributed across 830 high quality microbial genomes (N50 = 40,316, 2,396 Mb) (Lesker et al., 2020). A BLAST database was constructed containing these genomes and employing BLASTn with a cutoff of 95% to match an ASV to a given genome (Camacho et al., 2009). The BLASTn results were parsed using a script that relies on the biopython library (Cock et al., 2009) and SeqKit tool (Shen et al., 2016)² and microbial MAGs matching ASVs were parsed and extracted from the complete iMGMC fasta collection. Using this procedure, we matched 17 MAGs with male ASVs and seven MAGs with female ASVs.

The identified metagenomes-assembled genomes (MAGs) from male and females ASVs were merged as single metagenomes, respectively (hereafter, male and female metagenome), and each metagenome was processed using the Anvi’o platform for metagenomics (hope, v7.1) (Eren et al., 2015). Per the metagenome, a separated contig database was constructed, and then we profiled k-mer frequencies, functional annotation of open reading frames (ORFs), and GC content of every contig, respectively. We initially binned contigs exactly as in the MAGs collection and further refined these bins by manual human-guided binning to obtain >90% completion and < 10% redundancy values per bin. Binning refinement is not necessary for the female metagenome, given that both the completion and redundancy values for the female metagenome had already met the defined criteria during the MAGs creation process in Anvi’o, a process that involved contig separation and clustering based on tetranucleotide frequency and phylogenetic differentiation, the binning refinement was not necessary. In contrast, the male metagenome contained a higher number of MAGs than the female, suggesting greater microbial diversity and complexity in the male sample. As a result, the male metagenome required manual bin refinement. After refining the bins in the male metagenome, we identified the taxonomy of each bin according to the closest match with genomes from the GTDB database, consisting of 317,542 genomes (database release 07-RS207, 8th April 2022) (Chaumeil et al., 2022). Metaproteomes were resolved using prodigal v2.6.3 (Hyatt et al., 2010). We assessed MAG redundancy and completion using the CheckM tool (Parks et al., 2015). In each MAG, genes and proteins were annotated from each metagenome using NCBI’s Clusters of Orthologous Groups (COGs) database (Tatusov et al., 2000).

In both the female and male metagenomes, we employed the anvi’o function “anvi-get-sequences-for-gene-calls” with the flag --export-gff3 to export a single gff file belonging to each metagenome. Each gff file was parsed by the MAG using a modified bash script available at: https://groups.google.com/g/anvio/c/DWO8fDQ_g7M. In this manner, in both female and male metagenomes, we obtained gene calls per MAG in gff format, which were placed in a given folder and served as input to construct a pangenome matrix via Roary, using the following command: roary -e --mafft -p 30 (Page et al., 2015). As described in the Roary README, each pangenome matrix was plotted using the roary_plots.py script using as inputs the MAG phylogenetic tree (accessory_binary_genes.fa.newick) and gene_presence_absence.csv tabular file, both outputs of Roary.

To identify metabolic genes along metagenomes, we re-annotated all proteins from male and female metagenome databases using the KEGG ortholog assignment KOfam method (Aramaki et al., 2020), including the use of the KEGG BRITE Database (Kanehisa et al., 2017). Enriched metabolic modules were obtained by determining the presence of a given metabolic module in MAGs associated with sham or DCM and quantifying the distribution of this module across all MAGs belonging to a given metagenome. A module was considered to be enriched in sham or DCM-associated MAGs if the q-value of the Rao score test was below 0.05, as described here (Shaiber et al., 2020). The latter function is fully implemented in the anvi’o function “anvi-compute-metabolic-enrichment.” Statistically significant modules were plotted as Sankey plots using the R package network3D.³ The entire metabolic reconstruction was benchmarked using the gapseq program (Zimmermann et al., 2021). Contigs from all MAGs were annotated using a curated reaction database composed of the MetaCyc (Caspi et al., 2016), KEGG (Kanehisa et al., 2019) and ModelSEED (Henry et al., 2010) databases. Prediction of metabolic pathways and gap-filling per MAG was performed by running the gapseq_find.sh script with the following flags: -v 0 -b 200 -p all -t auto. To model the influx and efflux of metabolites, we chose as gapfill-medium a formulation made of glucose and a minimal amount of acetate, with no amino acids, into a BacArena modeling framework, as described here: https://bacarena.github.io/ (Bauer et al., 2017). The chosen parameters for the simulations were as following: arena size of 50 × 50 grid cells, random population of two grid cells as starting culture, removal of acetate from the initial substrate list allow the community to grow with bacterial acetate from cross-feeding, and performing simulations for 14-time steps, respectively. Parsing of all reconstructed pathways, including text-processing, was executed through the Jupyter Notebook using libraries such as NumPy (Harris et al., 2020) and pandas.⁴

Blood samples and cells from the spleen and colonic lamina propria, were collected and analyzed by flow cytometry 12 weeks after DCM induction, as previously described (Ulndreaj et al., 2016; Vidal et al., 2017; Ugalde et al., 2021). Red blood cells from the spleen and blood samples were lysed in red blood cell lysis buffer ACK (Thermo Fisher) for 5 min and viable cells were stained with a viability dye (Fixable viability dye Red 780, TonboBio) for 20 min. The following antibodies were used to distinguish granulocytes, monocytes, dendritic cells, and T lymphocytes: Ly6C (clone HK1.4) conjugated to PerCP/Cyanine5.5, Ly6G (clone 1A8) conjugated to APC, CD11b (clone M1/70) conjugated to FITC, CD3 (clone 17A2) conjugated to PerCP/Cyanine5.5, CD8 (clone 53–5.8) conjugated to FITC, CD4 (clone GK1.5) conjugated to APC, CD11c (clone N418) conjugated to FITC, and CD25 (clone 3C7) conjugated to PE (BioLegend, San Diego, United States). Matching isotype controls or fluorescence minus one (FMO) was used to gate the positive cells. Data were acquired using a BD FACS Canto II (BD Biosciences, San Jose, CA, United States) and analyzed using FlowJo X 10 (BD Biosciences, Ashland, United States). Samples that did not meet the minimum number of events (100) were excluded from analysis.

The animals were transcardially perfused with phosphate buffered saline (PBS). The spinal cords were dissected out (0.3 cm rostral and 0.3 cm caudal from the compression epicenter), and post-fixed with 4% paraformaldehyde (PFA) in PBS, and cryoprotected in 30% sucrose/PBS for 48 h. Coronal sections (20 μm thick) were prepared and blocked (5% normal serum and 0.3% triton X-100 in PBS) for 1 h at room temperature (RT). Incubation with the primary antibody GFAP (1,1,000, C9205, Cy3-conjugated, Sigma) and 4′,6-diamidino-2-phenylindole (DAPI, 1:200, Sigma), was performed at RT for 1 h. Sections were systematically sampled every 580 μm over 3,420 μm (13 sections per animal). Reconstruction was performed using ImageJ software bundled with 64-bit Java 8 (Schneider et al., 2012) and the related Fiji distribution (Schindelin et al., 2012) using built-in plugins, including the MosaicJ external plugin.

All statistical analyses were performed using the GraphPad Prism 5.01 software (GraphPad Software, Inc.). Data sets were analyzed for normal distribution using the D’Agostino-Pearson normality test. GI function was analyzed using two-way analysis of variance (ANOVA) for repeated measurements with Tukey’s multiple comparison test. Flow cytometry analyses and ladder walk test results were analyzed using the nonparametric Mann–Whitney U test. Data are presented as mean ± standard error of the mean (SEM). Differences were considered statistically significant at p ≤ 0.05.

To characterize the systemic changes in the blood during DCM, we divided the animals into four groups that received either sham or DCM surgery (Figure 1A). We used a previously characterized DCM mouse model that resembles the human pathology (Vidal et al., 2017; Laliberte et al., 2021; Ulndreaj et al., 2022). At 12 weeks post DCM-induction there is a 47% compression ratio in the mouse DCM model (Vidal et al., 2017), loss of motor neurons (Satkunendrarajah et al., 2018; Ulndreaj et al., 2022), and development of motor impairment (Vidal et al., 2017; Laliberte et al., 2021; Ulndreaj et al., 2022). Representative images of reconstructed spinal cords of sham and DCM groups after DCM confirmed spinal compression versus sham (Figure 1B, upper and lower panels), as well as increased astrogliosis (Figure 1B, upper and lower panels, small inserts), and locomotor deterioration (Figure 1C) resembling previously published findings from the DCM mouse model (Vidal et al., 2017; Desimone et al., 2021).

At the systemic level, T cells in DCM females presented a 1.7-fold increase compared to the sham group (Figure 2A; **p = 0.0022) and a 6-fold increase in the frequency of CD4⁺ T cells (Figure 2A; *p = 0.0411). In males, no significant changes were detected at the systemic level in the T cells (Figure 2B). Non-significant changes between the DCM and sham groups were observed in the CD8⁺ T cell subset, dendritic cells, granulocytes, monocytes, and monocyte subsets in both females and males (Figures 2A,B, respectively).

To assess gastrointestinal function and constipation in our DCM mouse model, we analyzed the total output of fecal pellets and the total colon length in the sham and DCM groups. We did not observe significant differences neither in the time course of fecal output over a 15-min time period or the colon length between the sham and DCM groups in females (Figure 3A) or males (Figure 3B). The gut wall is normally composed of a layer of epithelial cells (enterocytes), secretory cells, and the GALT. The latter mainly uses lymphocytes to carry on immune responses (Mowat and Agace, 2014) and communicates with the CNS by interacting with nerve fibers, secreting metabolites, and immune cells (Mowat and Agace, 2014). We previously reported that systemic T cell levels vary during DCM progression and following surgical decompression (Vidal et al., 2018; Ulndreaj et al., 2022). Thus, to evaluate the effect of DCM on T lymphocytes, we measured the frequency of T lymphocytes and three different subsets, CD4⁺ T cells, CD8⁺ T cells and regulatory T cells (Tregs) in the colonic lamina propria and spleen. There were no significant differences in the frequency of T lymphocytes, CD4⁺ T cells, and CD8⁺ T cells between the sham and DCM groups in both female and males. However, DCM females showed a significant increase (*p = 0.045) in the frequency of Tregs compared to sham females (Figure 3C). Non-significant differences were observed in males (Figure 3D). Similarly, in the spleen, T cells, CD4⁺ T cells, CD8⁺ T cells, and Tregs had similar frequencies between the sham and DCM groups in both females and males (Figures 4A,B). In the myeloid family, dendritic cells, neutrophils, and monocytes reached comparable levels in females (Figure 4A), whereas in males the frequency of monocytes there was a 1.2-fold increase (**p = 0.0076) in the DCM group compared with sham males (Figure 4B).

To investigate the effects of DCM on the mouse gut microbiota, we hypothesized that DCM causes gut dysbiosis in a mouse model. To address the latter, we designed an experimental approach in which males and females were randomized into two groups, either sham or DCM (Figure 1A). Fecal samples were longitudinally collected from the same animal at three selected time points: 6-, 9-, and 12-weeks post-DCM induction, herein referred to as early, mid, and late, respectively. The selected time points: early, mid, and late, were chosen based on previous work with the DCM mouse model that showed the presence of symptomatology (e.g., manual dexterity impairment, pain, and locomotor/gait deficits) (Vidal et al., 2017; Laliberte et al., 2021; Ulndreaj et al., 2022). To further characterize the gut microbial community in both females and males during DCM progression, we employed a 16S rRNA sequencing approach, in which DNA was extracted from fecal samples and the bacterial 16S SSU rRNA V4 region was amplified by PCR and further sequenced using Illumina technology, as described in methods section. As a control, we included a sham group for each week of dysbiosis. We sequenced five control samples and six treatment samples per week and per sex, respectively obtaining 66 samples in total (see Supplementary Table S1) which were input into the NextFlow ampliseq nf-core pipeline. First, DADA2 (Callahan et al., 2016) constructed a probabilistic model, based on the inferred sequencing error rates of reads that were used to denoise and collapse these reads as contigs, denominated Amplicon Sequence Variants (ASVs, see Supplementary Figure S1).

Phylogenetic assignament of non-redundant ASVs by DADA2, revealed the presence of different genera, particularly the bacterial members of Lachonospiraceae, Muribaculaceae, Oscillospiraceae, Ruminococcaceae, and Clostridia, among other families, respectively (Figure 5A and Supplementary Table S2). We plotted the abundance of each family across female and male samples, respectively using phyloseq ecologically-oriented ordination methods, based on the NeatMap approach (Rajaram and Oono, 2010). Notably, females did not show a significant difference in alpha diversity between sham and DCM groups, but males did: a clear separation between sham and DCM samples was not evident across female samples when abundances were plotted and ordered across families (Figure 5B, upper heatmap). In contrast, a clear ordination between sham and DCM samples was evident across male samples, when abundance was plotted across all detected families (Figure 5B lower heatmap, respectively). In females, alpha diversity was not significantly different between the sham and DCM groups (p = 0.438), with a tendency of time effect (Figure 5C, upper panel, p = 0.0845). However, a trend is observed at the early time point, where the Shannon diversity is higher in the sham group compared to DCM (p (Day:Treatment = 0.0845)). In males sham and DCM groups presented a clear separation (Figure 5C, lower panel, p = 0.0023) without a time effect (p = 0.4896). Although the time effect was not significant, there is a clear trend where Shannon diversity decreases over time in the sham group, while it increases in DCM. Consistent with the latter, samples ordered according to the Bray–Curtis dissimilarity in a two-dimension Non-metric Multi-dimensional Scaling (NMDS) plot, revealed spatially mixed sham and DCM samples in females, respectively (Figure 5D, left). In contrast, the male sham and DCM samples were spatially ordered with clear separation in the same plots (Figure 5D, right). Figure 5E shows bar plots of the relative abundance of microbial populations belonging to the top 10 genera across sex and disease progression (early, mid, and late DCM, respectively). A clear increase in the abundance of bacterial members belonging to Alloprevotella, Lachonospiraceae, and Helicobacter genera was observed across DCM states in both sexes, but the bacterial population dynamics were different between male and female DCM. In females, early after DCM induction, bacterial members belonging to genera such as Alloprevotella, Helicobacter, and Muribaculum increased in abundance, whereas Alistipes genus diminished in DCM compared to the sham group, except in the mid DCM condition (Figure 5E, left panel). Importantly, the abundance of Lachnospiraceae NK4A136 genus increased during DCM progression at mid and late time points compared to that in the control group (Figure 5E left panel). In contrast, changes in abundance in males were more pronounced than in females: the abundance of bacterial members belonging to Alloprevotella, Bacteroides, Helicobacter, Odoribacter, Lachnospiraceae NK4A136, and Prevotellaceae UCG-001 consistently increased in the DCM group compared to the sham group at all time points (Figure 5E, right panel). Overall, these findings confirm that dysbiosis induced by DCM is sex specific. Although there is a general effect of increased abundance of Alloprevotella and Lachonospiraceae genus in both sexes, more pronounced bacterial abundance changes were observed in male gut microbiomes. In conclusion, our initial analysis in both sexes showed significant changes in the composition of the murine gut microbiota between sham vs. DCM, indicating that DCM induced sexually dimorphic dysbiosis in mice.

We aimed to closely examine the bacterial changes caused by DCM in both male and female gut microbiomes. Our goal was to link each ASV identified during the 16 s rRNA sequencing with a microbial genome from the mouse gut. To achieve this, we used the mouse-gut dereplicated metagenome-associated genomes collection as our reference (iMGMC-hqMAGs-dereplicated_genomes.tar.gz) as the genomic database, a genome collection in FASTA format consisting of 132,958 contigs distributed across 830 high quality microbial genomes (N50 = 40,316, 2,396 Mb) (Lesker et al., 2020). Our primary data source was the DADA2_table.tsv file, which lists all ASV sequences from both sham and DCM conditions in both sexes. For this purpose, we employed a mouse-gut dereplicated metagenome-associated genomes collection. We, then we used the edgeRun package implemented in R to identify statistically significant ASVs that were differentially abundant when compared to sham versus DCM in both male and female samples. Specifically, we found 11 and 41 such ASVs in female and male samples, respectively (see Supplementary Table S3). BLAST results of these ASVs against the iMGMC contig collection revealed seven and 17 unique genomes with ≥95% 16S identity across female and male samples, respectively (Figures 6A,C, respectively, Supplementary Table S4). Using the anvi’o platform for metagenomics along with the GTDB database as a taxonomic reference (October 2022) we calculated N50, percentage of completion, and redundancy, among other parameters across male and female genomes, respectively (Supplementary Table S4). The percentage of completion of male genomes ranged from ~93% with a median N50 of ~38,000 while the percentage of completion of female genomes ranged from ~94% with a median N50 of ~107,000. Importantly, several female genomes belonging to Lachonospiraceae family were abundant in sham but not in DCM and bacteria belonging to Rikenellaceae, Oscillospiraceae, and Muribaculaceae families were exclusively present in female DCM samples, respectively (Figure 6A). GC content, tetranucleotide frequency and taxonomy clustering also supported this observation revealing three Lachonospiraceae genomes present in sham but not in DCM that clustered together and contained lower GC content than bacteria present only in DCM condition (Figure 6B). Of note, one of the two bacteria belonging to the Muribaculaceae family was present in DCM but not in sham (UBA7173 sp002491305), while the bacteria UBA7173 sp001689485 was present in sham but not in DCM (Figure 6B).

The latter observations can be applied to male genomes, with a clear presence of bacterial members belonging to Lachonospiraceae family present in sham, but not in DCM. Conversely, several members of Muribaculaceae family were exclusively present in the DCM group, but not in the sham group (Figure 6C). In males, the phylogenetic separation between bacteria presents in sham but not DCM (and vice versa) was not obvious, but bacterial members of Muribaculaceae family including the Paramuribaculum genus, abundant in DCM, tended to have a higher GC content than bacteria from Lachonospiraceae family that is present in sham but not in DCM (Figure 6D). As observed in females, a male bacterium belonging to Muribaculaceae family, genus UBA7173 (UBA7173 sp001915385), was present in the sham but not in the DCM condition, whereas two bacteria belonging to the same Muribaculaceae family (Paramuribaculum sp001689565 and another bacterium without assigned genus) were present in DCM, but not in sham condition (Figure 6D).

In conclusion, there’s a noticeable shift in GC-rich bacteria from the Muribaculaceae family, present in DCM but not in the sham condition across both sexes, whereas several bacterial members of Lachonospiraceae along with other families tended to be absent when DCM-induced dysbiosis prevailed. Of note, in males, bacteria belonging to Muribaculaceae family were absent in the initial phyloseq analysis, because the ASVs associated with these bacteria did not have an assigned genus (Supplementary Table S4).

To discover conserved gene clusters of genes associated with bacteria present in DCM-induced dysbiosis, we obtained all gene calls from male and female metagenome in gff format in and we parsed gene calls from each MAG. In this manner, we obtained all genes belonging to each MAG in gff format, respectively (see “Pangenome reconstruction of metagenomes” in Materials and Methods section). Then, in both male and female metagenomes, we calculated each pangenome using the Roary pipeline (Page et al., 2015) using gene calls from all MAGs belonging to a given metagenome. From Figure 6A, we considered four and three MAGs in females exclusively present in sham and DCM condition, respectively while from Figure 6B, we considered eight and four MAGs in males exclusively present in sham and DCM condition, respectively (see Supplementary Table S5, sections “female_venn” and “male_venn,” respectively).

In the context of pangenomics, the core genome refers to genes that are present in all strains of a set of analyzed genomes, highlighting the essential functions. On the other hand, the soft-core genome, includes genes found in a high percentage, typically over 95%, of the strains (these might be essential for most but not all strains). Lastly, the cloud genome encompasses genes present in a small subset of the analyzed strains, usually fewer than 15%. These genes shed light on the unique adaptations and potential functionalities of certain strains. In our analysis, almost all gene clusters fell under the umbrella of the cloud genome and none of them were classified as soft-core or core genomes in both male and female pangenomes, despite the close phylogenetic relationship of some MAGs belonging to the families Lachonospiraceae. and Muribaculaceae, present in both pangenomes (see Figures 7A,B, for female and male pangenome matrix, respectively and Supplementary Table S5, sections “common_set_males” and “common_set_females,” respectively). Nevertheless, male MAGs belonging to Muribaculaceae family and Paramuribaculum genus associated with DCM-induced dysbiosis, presented several conserved genes involved in pyruvate pathways (see Supplementary Table S5, section “common_set_males”).

Since we did not find common gene clusters within MAGs associated with DCM-induced dysbiosis in males and females, we asked whether metabolic pathways were common in MAGs present in DCM. Using the anvi’o genomic annotation per genome as starting point, we re-classified each gene using the KEGG ortholog assignment KOfam method (Aramaki et al., 2020), relying on the KEGG BRITE Database (Kanehisa et al., 2017). Then, we detected enriched metabolic modules in MAGs associated with sham or DCM using the anvi’o function “anvi-compute-metabolic-enrichment” as described in Material and Methods, section “Metabolic genes and network reconstruction.” In females, several statistically significant pathways within MAGs associated with DCM were found, including those related to lipid metabolism and the reductive pentose phosphate cycle (Alistipes sp009774895), F-type ATPase, glycogen degradation, and lysine biosynthesis (Lawsonibacter, UBA7173 sp001689485). Notably, the latter three pathways were also present in MAGs associated with sham (see Sankey plot of Figure 7C, and Supplementary Table S6).

In males, more statistically significant pathways were identified. Pathways related to Phosphoribosyl pyrophosphate (PRPP) biosynthesis, ADP-L-glycero-D-manno-heptose biosynthesis, and Nucleotide sugar biosynthesis, among others, were associated with MAGs preferentially present in sham rather than DCM (see Sankey plot of Figure 7D, left and Supplementary Table S7). Regarding pathways present in MAGs associated with DCM-induced dysbiosis, pathways related to lipid A biosynthesis and UDP-L-FucNAm (UDP-N-acetyl-beta-L-fucosamine) biosynthesis were enriched in a MAG belonging to Muribaculaceae family, both pathways leading to the biosynthesis of Lipopolysaccharide (LPS) (Mulrooney et al., 2005; Emiola et al., 2014). In addition, the production of Propanoyl-CoA metabolism is exclusively found in this MAG. A pathway involving Thiamine biosynthesis involved in bacterial butyrate production (Soto-Martin et al., 2020; Park et al., 2022) was enriched in Paramuribaculum sp001689565 and a pathway related to the production of Pyridoxal phosphate biosynthesis enriched in Kineothrix sp000403275 is associated with pathogenesis (Dick et al., 2010; El Qaidi et al., 2013) (see Sankey plot of Figure 7D, right and Supplementary Table S7). This evidence suggests that in both sexes, several pathways are involved in DCM-induced dysbiosis, but other common molecular mechanisms may be present.

To address this, we employed the gapseq program to combine the identified metabolic pathways along metabolic network reconstruction using gap-filling reactions, as described in Materials and Methods (Zimmermann et al., 2021). Male and female reconstructed pathways were parsed according to MAGs present in the sham or DCM groups and are available in Supplementary Table S8. Regarding central metabolism reactions, we reviewed all pyruvate associated pathways across MAGs present in the sham or DCM groups, across both sexes. In female MAGs, we evidenced that the pyruvate fermentation to propanoate and pyruvate decarboxilation to acetyl-CoA reactions were prevalent in DCM-associated MAGs and nearly absent in sham-associated MAGs (Figure 7E). The latter observations can also be applied across male MAGs (Figure 7F). The reactions involved in this process are summarized in Supplementary Table S9. The reaction in which pyruvate ferment to propanoate is complex and involve a fumarase, a fumarate reductase, a propionyl-CoA:succinate CoA transferase, a methylmalonyl-CoA carboxyltransferase and a malate dehydrogenase while the reactions in which pyruvate decarboxile to acetyl-CoA involve the presence of a pyruvate dehydrogenase NADP⁺ (Supplementary Table S9). Conversely, the conversion of pyruvate to butanoate involve three enzymes including a acetyl-CoA-acetyltransferase PaaJ and a Butyrate kinase (Buk). The latter reaction is frequently absent in MAGs present in DCM-induced dysbiosis across both sexes (Figures 7E,F and Supplementary Table S9). Since bacteria can ferment pyruvate down into butyrate via Butyrate kinase activity (Sievers et al., 2019) (diminished in DCM) and shuttle pyruvate into lactate through pyruvate dehydrogenase enzyme (Esquivel-Elizondo et al., 2017) (augmented in DCM), this imply that in the transition from sham to DCM-induced dysbiosis, the pool of butyrate among other metabolites in the microbiome such as propionate and acetate (Louis and Flint, 2017) could largely vary in DCM-induced dysbiosis.

Based on the metabolic reconstruction across sham and DCM-associated MAGs in both sexes, we interrogated whether acetate, propionate and butyrate as fermentation products vary their concentrations between sham and DCM-induced dysbiosis conditions, respectively. This task is complicated in terms of predicting metabolite influx and efflux based solely of bioinformatic gene predictions, due complex metabolic networks and numerous biochemical transformations occurring in a microbiome community. Using gapseq modeling with a gap-filling approach, we simulated at community level the production of the mentioned metabolites using minimal medium composed of glucose and acetate. Across female MAGs present in sham condition, a steady increase of butyrate and acetate concentration was registered along the simulation’s steps with a marked growth increase of UBA7173 sp002491305 (Muribaculaceae family) and Clostridia sensu stricto MAGs (Figure 8A). Conversely, when MAGs associated with DCM-induced dysbiosis were modeled under the same conditions, only acetate was a strongly changed substance with marked growth increase of all MAGs (Figure 8B). Same observations apply in the case of male MAGs modeling. Across male MAGs present in sham condition, both acetate and butyrate metabolites sharply increase along with marked growth of UBA3282 and UMGS1004 sp003334805 MAGs (Figure 8C). Conversely, under DCM-induced dysbiosis, only acetate increased its concentration displaying a marked growth increase of Paramuribaculum sp001689565 (Muribaculaceae family) and Muribaculaceae MAGs (Figure 8D). In these simulations, butyrate metabolite production markedly increases when sham-associated MAGs are present while its production is nearly absent during DCM-induced dysbiosis, irrespective of the sex status. We interrogated the top 10 produced substances predicted by gapseq gap-filling modeling across male and female MAGs associated with sham (left) or DCM condition, respectively (Figure 8E). For better visualization, the predicted measurements in millimoles per gram dry weight per hour (mmol/(gDW*hr)) were normalized using a StandardScaler method in python. In this setting, across both males and females MAGs, the production of acetate is inversely proportional to the production of butyrate, and inversely proportional to the production of CO₂ in less extent. This observation also supports that acetate can be a substrate for butyrate production, since as it has been described that butyrate-producing bacteria found in the human colon can be net users of acetate (Duncan et al., 2004). Overall, our metabolic simulations support the sustained production of butyrate in sham condition but in less extent in DCM condition, with the concomitant use of acetate, not produced by bacteria that synthetize butyrate in the mouse colon.