A major impact of maternal microbiota in the fetal intestine was indicated by differential gene expression analysis (Table 1 , Additional file 1 : Table S1, Fig. 2 , and Additional file 2 : Fig. S1–S2). GF and SPF expression profiles clustered in principal component analysis (PCA; Fig. 2 A). The groups were completely separated in a supervised orthogonal partial least squares discrimination analysis (OPLS-DA, validated by R 2 cum and Q 2 cum scores and permutation test; Fig. 2 B).
Table 1 Summary of differential expression analysis of GF vs SPF mouse fetal tissues Tissue Genes passing pre-filtering Downregulated in GF ( p.adj < 0.05) Upregulated in GF ( p.adj < 0.05) Fetal intestine 13,493 1249 (9.3%) 941 (7.0%) Fetal brain 14,645 67 (0.46%) 50 (0.34%) Placenta 13,433 409 (3.0%) 309 (2.3%) Fig. 2 Differential gene expression analysis of GF versus SPF fetal intestine. A PCA of 1000 most variable genes. The ellipsoid shows Hotelling’s T2 (95%). B OPLS-DA of 1000 most variable genes. C ORA of genes which were significantly upregulated or downregulated in GF versus SPF ( p.adj 1 and negative or positive S-plot loadings. The top 20 enriched ontology terms are shown; for more details, see Additional file 2 : Fig. S1–2. D Hallmark gene sets enriched in GSEA. A maximum of 10 of the top gene sets with FDR q -value < 25% are shown with negative (blue) and positive (red) net enrichment scores (NES) for GF versus SPF. E Volcano plot of differential gene expression. Genes with negative log 2 fold change values were downregulated in GF fetuses. Dashed lines indicate log 2 fold change ± 0.5 and adjusted p value 0.05
Summary of differential expression analysis of GF vs SPF mouse fetal tissues
Differential gene expression analysis of GF versus SPF fetal intestine. A PCA of 1000 most variable genes. The ellipsoid shows Hotelling’s T2 (95%). B OPLS-DA of 1000 most variable genes. C ORA of genes which were significantly upregulated or downregulated in GF versus SPF ( p.adj 1 and negative or positive S-plot loadings. The top 20 enriched ontology terms are shown; for more details, see Additional file 2 : Fig. S1–2. D Hallmark gene sets enriched in GSEA. A maximum of 10 of the top gene sets with FDR q -value < 25% are shown with negative (blue) and positive (red) net enrichment scores (NES) for GF versus SPF. E Volcano plot of differential gene expression. Genes with negative log 2 fold change values were downregulated in GF fetuses. Dashed lines indicate log 2 fold change ± 0.5 and adjusted p value 0.05
Over-representation analysis (ORA) of differentially expressed (DE) genes ( p.adj < 0.05) by Gene Ontology (GO) terms suggested an impact of maternal microbiota on the fetal intestinal immune system, epithelial cell physiology, and translation (Fig. 2 C; for more details, see Additional file 1 : Fig. S1–S2 and Additional file 1 : Table S4). The GO categories response to virus , brush border , lytic vacuole , lipid localization , microvillus , and response to interferon beta were most strongly over-represented among genes downregulated in the GF intestine (Fig. 2 C). Amide metabolic process , ribonucleoprotein complex biogenesis and binding , and rRNA binding were over-represented among genes upregulated in GF. Over-representation analysis (ORA) of genes with high variable importance in projection (VIP) in OPLS-DA showed similar enrichment.
Gene set enrichment analysis (GSEA) using the hallmark gene sets indicated downregulation of the interferon alpha and interferon gamma response gene sets in the GF fetal intestine (NES < 0, FDR < 5%; Fig. 2 D, Additional file 1 : Table S5a). The leading-edge genes included Irf7 , Irf9 , and Stat transcription factors; B2m (a component of the MHC complex); and interferon-stimulated genes (ISGs) such as Rsad2 , Ifi44 , and Oasl (Additional file 2 : Fig. S5). Negative enrichment scores of gene sets related to interferon signaling and virus response were also seen in the GF fetal intestine in the C2 curated pathways and C7 immunosignature gene sets (Additional file 1 : Table S5b). Proliferation-associated gene sets (such as Myc targets , E2F targets , G2M checkpoint ) and unfolded protein response were upregulated in GF versus SPF fetuses (NES > 0; Fig. 2 D). The hallmark gene sets do not cover all of the GO categories enriched in ORA.
A higher number of genes was downregulated than upregulated in the GF fetal intestine (Table 1 , Fig. 2 E). These were enriched for immunity-related gene sets in ORA (Fig. 2 and Additional file 2 : Fig. S2). Strikingly, 27 out of the 30 genes most strongly downregulated in the GF fetal intestine were involved in immunity, intestinal host-microbe interactions, and xenobiotic metabolism ( p.adj < 0.05 and ranked by fold change). The genes significantly downregulated in GF fetuses included mucins, interferon and cytokine signaling genes such as Stat1-3 ; virus response genes such as Trim and Oasl families and Rsad2 ; interleukins and their receptors (such as Il18 , Il18R , Il34 , and Il10rb ); the major acute phase response gene Saa1 ; several complement genes; 7 MHC-I genes (but none of the MHC-II genes); antimicrobial lectins ( Reg3b , Reg4 , Lgals3 , Lgals4 , Lgals8 , Lgals9 ); the antimicrobial peptide Ang ; and the marker of commensal microbiota associated regulatory T cells and ILC3 cells Rorc . Tight junction component ( Cldn15 , Cldn19 , Tjp3 ) and enteroendocrine cell marker genes Insl5 , Pyy , Gip , and Nts were also downregulated. However, the expression of interferon genes was not detectable in the fetal intestine. In contrast, 80% of the genes annotated for translation, ribosomes, and tRNA metabolism were upregulated in the GF fetal intestine, including several subunits of the translation initiation factor eIF2 and the elongator acetyltransferase complex. The genes most strongly upregulated in the GF fetal intestine included metallothioneins Mt1 and Mt2 , the B cell chemoattractant Cxcl13 , developmental regulators Foxl1 and Hoxc5 , collagen biosynthesis genes P4ha2 and Col11a1 , and various genes involved in metabolism.
Regarding the predicted transcription factor binding sites, the DE genes in the intestine (and in the brain and placenta) were most significantly enriched for E2Fs, ZF5, and FOXN4 (Additional file 1 : Table S4). Predicted targets of interferon regulatory factors (IRFs) and early growth response proteins (EGRs) were also significantly enriched. DE genes in the fetal intestine (but not in the brain and placenta) were significantly enriched for multiple predicted binding sites of AhR and AhR nuclear translocator (Arnt). DE genes in the intestine and placenta (but not in the brain) were also significantly enriched for binding sites for the vitamin D receptor (VDR) and/or FXR.
In the fetal brain, differences in gene expression profiles were clearly less prominent. GF and SPF fetuses were not discriminated by PCA (Fig. 3 A), and no valid OPLS-DA model separating the groups was obtained (best model: 1 + 7 components; R 2 Xcum 0.80, Q 2 cum 0.05, Q 2 permutation intercept −0.27). ORA indicated significant enrichment in GO categories related to neural functions (such as myelin sheath , glial cell projection , modulation of chemical synapse , and axon ) as well as antiviral immunity among genes downregulated in the GF fetal brain (Fig. 3 B and Additional file 2 : Figs. S1, S3). Some neural development categories were upregulated in GF fetuses. Fig. 3 Differential gene expression analysis of GF versus SPF fetal brain. A PCA of 1000 most variable genes. The ellipsoid shows Hotelling’s T2 (95%). B ORA of genes which were significantly upregulated or downregulated in GF versus SPF ( p.adj < 0.05; DE GF up and DE GF down). The top 20 enriched ontology terms are shown; for more details, see Additional file 2 : Fig. S1, S3. C Hallmark gene sets enriched in GSEA. A maximum of 10 of the top gene sets with FDR q -value < 25% are shown with negative (blue) and positive (red) net enrichment scores (NES) for GF versus SPF. D Volcano plot of differential gene expression. Genes with negative log 2 fold change values were downregulated in GF fetuses. Dashed lines indicate log 2 fold change ± 0.5 and adjusted p value 0.05
Differential gene expression analysis of GF versus SPF fetal brain. A PCA of 1000 most variable genes. The ellipsoid shows Hotelling’s T2 (95%). B ORA of genes which were significantly upregulated or downregulated in GF versus SPF ( p.adj < 0.05; DE GF up and DE GF down). The top 20 enriched ontology terms are shown; for more details, see Additional file 2 : Fig. S1, S3. C Hallmark gene sets enriched in GSEA. A maximum of 10 of the top gene sets with FDR q -value < 25% are shown with negative (blue) and positive (red) net enrichment scores (NES) for GF versus SPF. D Volcano plot of differential gene expression. Genes with negative log 2 fold change values were downregulated in GF fetuses. Dashed lines indicate log 2 fold change ± 0.5 and adjusted p value 0.05
Among the hallmark gene sets with negative enrichment scores in GSEA for the fetal GF brain, the most significant were interferon alpha and interferon gamma response sets (Fig. 3 C). The mitotic spindle gene set was most positively enriched in GF versus SPF mice (Fig. 3 C).
Also in the brain, the genes downregulated in GF fetuses were enriched for immunity-related genes (Additional file 1 : Table S4). The most strongly differentially expressed genes included interferon and virus response genes such as Rsad2 , Ifi44 , Ifi27 , Rtp4 , Trim30a , and Lgals3bp ; the glymphatic aquaporin Aqp4 ; and the lymphocyte antigen Ly6a (Fig. 3 D). The expression of interferon genes was not detectable. Multiple genes involved in neuronal development and synaptic signaling were significantly differentially expressed. Nrgn , two synaptotagmins, Shisa6 , Calm1 , the monoamine oxidase gene Maob , and the glial-specific Gfap and Olig2 were significantly downregulated. Genes upregulated in the GF brain included the neural stem cell regulator Phf21b ; the neural transcription factors and activators Nfib , Sox11 , and Eomes ; several cadherins involved in neural system development; and the lncRNAs Snhg6 and Gm47283 .
In the placenta, GF and SPF gene expression profiles were not separated in PCA but could be discriminated by a validated OPLS-DA model (Fig. 4 A, B). ORA indicated differences in gene sets related to development, tissue homeostasis, and immunity. GO categories positive regulation of cell death , regulation of apoptotic signaling pathway , cell adhesion molecule binding , and tube morphogenesis were most significantly enriched among genes downregulated in the GF placenta; collagen-containing extracellular matrix and multicellular organismal-level homeostasis were enriched among upregulated genes (Fig. 4 C; Additional file 2 : Figs. S1, S4). Fig. 4 Differential gene expression analysis of GF versus SPF placenta. A PCA of 1000 most variable genes. The ellipsoid shows Hotelling’s T2 (95%). B OPLS-DA of 1000 most variable genes. C ORA of genes which were significantly upregulated or downregulated in GF versus SPF ( p.adj 1 and negative or positive S-plot loadings. The top 20 enriched ontology terms are shown; for more details, see Additional file 2 : Figs. S1, S4. D Hallmark gene sets enriched in GSEA. The gene sets with FDR q -value < 25% are shown. E Volcano plot of differential gene expression. Genes with negative log 2 fold change values were downregulated in the GF placenta. Dashed lines indicate log 2 fold change ± 0.5 and adjusted p value 0.05
Differential gene expression analysis of GF versus SPF placenta. A PCA of 1000 most variable genes. The ellipsoid shows Hotelling’s T2 (95%). B OPLS-DA of 1000 most variable genes. C ORA of genes which were significantly upregulated or downregulated in GF versus SPF ( p.adj 1 and negative or positive S-plot loadings. The top 20 enriched ontology terms are shown; for more details, see Additional file 2 : Figs. S1, S4. D Hallmark gene sets enriched in GSEA. The gene sets with FDR q -value < 25% are shown. E Volcano plot of differential gene expression. Genes with negative log 2 fold change values were downregulated in the GF placenta. Dashed lines indicate log 2 fold change ± 0.5 and adjusted p value 0.05
In GSEA, only two hallmark gene sets were negatively enriched in the GF placenta ( IL6-JAK-STAT-signaling and TGF-beta-signaling ; Fig. 4 D). No significantly enriched hallmark gene sets in GF over SPF mice were detected in the placenta.
The genes most strongly downregulated in the GF placenta (Fig. 4 E) included immunoglobulin and complement genes, the recently characterized bacterial response gene AW112010 , the interleukin receptor Il1r2 implicated in endometriosis, the diamine oxidase Aoc1 implicated in pregnancy regulation, the prolactin precursor Prl3d1 , the exosomal endonuclease Endod1 , the neuroendocrine peptide Gal , the mitogen Mdk , the mitochondrial pyruvate dehydrogenase kinase Pdk4 , and the lncRNAs Gm7932 and 4933417E11Rik . The expression of the interferon genes Ifnk and Ifne was detectable in the placenta, but it was not statistically significantly stronger in the SPF animals. A majority (23 out of 30) of the most strongly differentially expressed genes were upregulated in the GF placenta. The most upregulated genes included the placental gene regulator Gcm1 , the trypsin inhibitor Pi15 , the killer cell receptor Klra4 , typical epithelial and endothelial genes, hemoglobins, and solute carriers.
Only seven genes were observed to be significantly differentially expressed in all three tissues: the insulin-degrading enzyme Ide ; the virus response genes Rtp4 , Rsad2 , and Isg15 ; the transcription regulator Btaf1 ; the Rho GTPase Rhou ; and the tubulin Tuba4a . A total of 26 genes were significantly differentially expressed both in the fetal intestine and brain, and 129 genes were shared between the intestine and placenta. All these gene sets shared between the tissues were enriched for GO terms for immunity, virus response, and symbiont interaction (not shown). However, more specific GOs (< 500 genes) were largely different.
The DE analysis comparing GF versus SPF fetuses was controlled for sex. When GF and SPF fetuses were compared in each sex separately, male fetuses appeared more sensitive to modulation by maternal microbial status. In the intestine, 1223 genes were significantly differentially expressed in male GF versus SPF fetuses, while only 307 genes were significantly DE in female GF versus SPF fetuses. In the brain, 8 genes were DE in males versus 2 genes in females, and in the placenta, 145 versus 36 genes. The total numbers of significantly DE genes were smaller in these analyses due to the smaller comparison groups. Genes which were significantly DE in GF versus SPF in the male intestine or placenta but not in females were enriched for development, biosynthesis, and histone methylation gene sets (not shown).
We also compared the gene expression in all male fetuses versus female fetuses, controlling for the differences in GF versus SPF animals (Additional file 1 : Table S3). In the fetal brain, all the significantly DE genes were sex chromosomal. In the intestine, several autosomal interferon-inducible genes were significantly upregulated in the male fetuses in comparison with females. In the placenta, 57 autosomal genes were significantly DE in male versus female fetuses; these were enriched for development-related genes but not immunity by ORA.
To explore the potential effects of microbial metabolites on the fetal intestine, brain, and placenta, we analyzed the associations between gene expression and metabolites. We utilized our previously published metabolomics data [ 14 , 23 ] and focused on metabolites which were undetectable in GF fetuses or significantly less abundant than in SPF fetuses. In total, 2200 molecular features were included in the analysis, being significantly more abundant in SPF fetuses in at least one tissue. Ninety-nine of these were only detected in SPF mice. To detect various types of associations between individual metabolites and genes, and between molecular feature groups and gene families or co-regulated pathways, we evaluated (1) correlations of abundances of individual metabolites with expression levels of individual genes, (2) correlations of clusters of metabolites with clusters of genes (by hierarchical clustering), and (3) correlations within biclusters composed of metabolites and genes.
We focused on metabolites and genes which showed significant differences between the experimental groups and strong correlations also within the SPF group (Spearman ρ > 0.9), thus representing metabolite-gene associations which indicate dose–response and are likely due to actual abundances of the compounds, rather than due to the overall differences of GF/SPF physiology. As the small size of the fetal organs required gene expression and metabolite profiling to be done from different fetuses, the association analysis was performed using the averages of two fetuses from each dam. All twin pairs were highly similar in terms of gene expression profiles (Spearman ρ range = 0.969–0.994) and metabolite profiles ( ρ = 0.926–0.991) and significantly more similar than fetuses from different dams (transcriptomes: ρ = 0.897–0.995, significantly lower at p < 0.001; metabolomes: ρ = 0.655–0.981, p < 0.001). Thus, examining the associations calculated across littermates is justified. The approach emphasizes associations to maternally derived metabolites, as the exposure to those is expected to be similar in the fetuses from the same litter.
All molecular features were scored based on significant hits in ORA of the strongly associated genes or gene clusters (Additional file 1 : Table S4). Highest-scoring molecular features were then examined in more detail by ORA of genes which were strongly directly associated with them (Fig. 5 ; for more detail, see Additional file 2 : Figs. S6–S8). An overview of the metabolite-gene associations is shown in Additional file 1 : Table S7. In the intestine and placenta, the expression levels of immunity genes (by GO annotations) were mostly positively correlated with metabolomics signal intensities. In contrast, genes associated with translation correlated mostly negatively with metabolites in the fetal intestine. In the brain, such differences were not observed, but genes with GO annotations for immunity or neurophysiology had more metabolite associations. Fig. 5 Over-representation analysis of genes strongly associated with metabolites in the fetal intestine, brain, and placenta. Highest-scoring metabolites missing from GF fetuses (in bold text) and highest-scoring annotated metabolites more abundant in SPF fetuses (in regular text) are shown. R, RP column; H, HILIC column (positive / negative)
Over-representation analysis of genes strongly associated with metabolites in the fetal intestine, brain, and placenta. Highest-scoring metabolites missing from GF fetuses (in bold text) and highest-scoring annotated metabolites more abundant in SPF fetuses (in regular text) are shown. R, RP column; H, HILIC column (positive / negative)
Examples of associations between individual molecular features and genes are shown as scatterplots in Fig. 6 . Fig. 6 Examples of associations between metabolites and gene expression in the fetal intestine, brain, and placenta. Only one gene per metabolite is shown. Metabolite signal intensity (ion abundance) on the X -axis and gene expression (normalized counts) on the Y -axis. Each data point is the mean of the two fetuses from one dam. Red = GF; turquoise = SPF. Tentatively annotated metabolites are marked with “?”. R, RP column; H, HILIC column (positive/negative)
Examples of associations between metabolites and gene expression in the fetal intestine, brain, and placenta. Only one gene per metabolite is shown. Metabolite signal intensity (ion abundance) on the X -axis and gene expression (normalized counts) on the Y -axis. Each data point is the mean of the two fetuses from one dam. Red = GF; turquoise = SPF. Tentatively annotated metabolites are marked with “?”. R, RP column; H, HILIC column (positive/negative)
A total of 148 metabolites which were not detectable in the GF fetal intestine or brain showed associations with gene expression profiles (Additional file 1 : Table S6). Eighty-four of these were also missing from the GF placenta. Twenty-three of the molecular features missing from GF fetuses were at least putatively characterized. Six of these were aryl sulfates (4-hydroxybenzenesulfonic acid, indoxyl sulfate, pyrocatechol sulfate, and phenyl sulfates). The aryl sulfates were among the metabolites with the strongest associations to the gene expression patterns in all tissues investigated. In the ORA of the fetal intestine, they were associated with immunity (response to biotic stimulus and activation of innate immune response), lipid metabolism, regulation of cell growth, aromatic compound biosynthesis, and tRNA metabolism. In the brain, aryl sulfates were associated with viral infection pathways, dopaminergic synapse, Ras signaling, and response to organocyclic compounds. In the placenta, they did not show significant associations.
3-Indolepropionic acid (IPA) and two unknown compounds (all observed only in SPF fetuses) were associated with the adaptive immune system and Ras signal transduction in the brain. In the fetal intestine, IPA was primarily associated with RNA metabolism, and in the placenta, with the regulation of cell growth.
Kynurenine, tryptophan, and 3-methylhistidine were associated with immunity (largely positively) and RNA metabolism (mostly negatively) in the intestine. 1-Methylhistidine was associated with virus response in the placenta. Several other amino acids and their derivatives were also significantly more abundant in the SPF fetuses.
Aryl sulfates and tryptophan derivatives are typical AhR ligands. As fetal intestinal DE genes were enriched for predicted AhR/Arnt targets, we specifically looked at associations of such annotated metabolites with transcription factor binding sites. Genes strongly associated with tryptophan, kynurenine, IPA, a phenyl sulfate, and several metabolites with tentative annotations as tryptamine and indoleacetic acids were significantly enriched for predicted AhR/Arnt binding sites. This was not observed for indoxyl sulfate, hydroxyindoleacetic acid, 4-hydroxybenzenesulfonic acid, and pyrocatechol sulfates. Instead, indoxyl sulfate, 4-hydroxybenzenesulfonic acid, pyrocatechol sulfates, and metabolites tentatively annotated as indoleacetic acids were associated with predicted VDR targets.
The betaine trimethylamine N-oxide (TMAO; only detected in the SPF fetuses) was associated with brush border, absorption, and lipid metabolism in the SPF fetal intestine. In the brain, the only strong correlation (positive) was with the Ide gene. In the placenta, there were no significant ORA hits for genes strongly associated with TMAO. 5-AVAB was associated with response to biotic stimulus and lipid metabolism in the intestine.
Butyrylcarnitine and several unannotated molecular features were associated with immunity (response to biotic stimulus, virus response, neutrophil degranulation, interleukin-1 production; including the Rorc gene) and lipid metabolism in the intestine. Also, in the placenta, butyrylcarnitine was associated with neutrophil degranulation and with oxidoreductase and peroxidase activities.
A bile acid (with retention time matching the secondary bile acid chenodeoxycholic acid), N-linoleyltyrosine or its isomer, and peiminine or N-oleoylphenylalanine (putative classifications) and two unannotated features were strongly associated with multiple pathways in the intestine. These included immunity (response to biotic stimulus, innate immunity, pattern recognition receptor signaling, T lymphocyte differentiation, several virus response pathways), translation, brush border, and intestinal absorption. These compounds showed no significant ORA hits in the brain or placenta.
The dipeptides Glu-Trp and Glu-Tyr were associated with extracellular response stimulus, virus responses, and other immunity pathways in all tissues investigated.
The largest number of direct associations with genes (306 genes with Spearman ρ > 0.9 in SPF mouse intestine data) were observed for an unidentified metabolite with the probable formula C3H4O5 (possibly 2-hydroxypropanedioic acid alias tartronic acid), which was not detected in the GF intestine but was present in brain and placenta in both experimental groups. In ORA, for the fetal intestine, the associated genes were significantly enriched for immunity pathways (response to virus, cytokine signaling, and adaptive immune response).
Several unidentified compounds were associated with virus response, translation, and ribonucleoprotein pathways in the intestine. One of these (RP +
[email protected]) matched the MS1 mass of queuine. This compound was undetectable in the GF intestine and placenta and significantly less abundant in the GF brain. It was associated with genes primarily enriched for translation and its initiation, ribosomes (including the aminoacyl-tRNA binding protein Rpl8), protein- l -isoaspartate ( d -aspartate) O-methyltransferase activity, nonsense-mediated decay, and E2F-1 targets in the intestine. In the brain, there were multiple unannotated molecular features which were more abundant in SPF fetuses and associated with neural development, synaptic function, and interferon/virus responses. In the placenta, several unidentified compounds were associated with immune responses (inflammatory response, cytokine production, and phagocytosis).
Some genes showed unexpected nonlinear dependencies on metabolites. For example, Slc25a20 was inversely correlated with 4-hydroxybenzenesulfonic acid in the SPF fetal intestine, but the gene was downregulated in the GF intestine, where this metabolite was undetectable (Fig. 6 ; last plot for intestine). This pattern was most commonly observed for aryl sulfates, trimethylated compounds, and other metabolites not detected in GF fetuses. It was more common in the brain than in the intestine or placenta. In the fetal intestine, genes with such nonlinear dependencies were enriched for negative regulation of innate immune response, lipid metabolism, and negative regulation of Notch signaling (not shown). Genes with linear correlations were enriched for positive regulation of defense response, complement activation, and translation. This was not observed in the brain or placenta.