How early gut bacteria growth relates to brain barrier and thinking development

We examined offspring behavior at four weeks of age. We did not find significant sex differences in any of the behavioral outcomes we measured (Table S2). We first assessed associative learning and memory in the MTerm and MPreterm groups using the contextual and cued fear conditioning test (n_MTerm = 13, n_MPreterm = 22). There was a significant difference in freezing time between MTerm (91.1 ± 28.99 s) and MPreterm (24.71 ± 4.56 s) mice in contextual testing (same environment and no acoustic conditioned stimulus (CS) presentation) (p = 0.0446). Furthermore, MTerm mice had significantly longer freezing time (147.8 ± 6.299 s) than MPreterm mice (126.1 ± 6.906 s) in the second half of the cued fear conditioning test (changed environment and acoustic CS presentation) (p = 0.0421) (Figure 1(b)).

We then employed the Morris water maze to assess spatial learning and memory in the two experimental groups (n_MTerm = 9, n_MPreterm = 16). MTerm and MPreterm mice both showed decreased latencies in finding the escape platform from days one to four (p < 0.0001), demonstrating learning in both groups. Notably, MTerm mice overall took less time to find the escape platform (p = 0.0465). There was no significant interaction (p = 0.539). These findings were determined by two-way ANOVA with predictors of training day and treatment (Figure 1(c)).

In an open field test to assess anxiety-like behaviors, MPreterm offspring demonstrated a reluctance to explore the open space (n_MTerm = 12, n_MPreterm = 22). They had significantly higher overall immobile time (63.85 ± 7.147 s) than that of MTerm (39.37 ± 6.859 s) mice at four weeks of age (p = 0.0276). Further analysis demonstrated that even though total time spent in the border zone of the open field was not significant between MTerm (1033 ± 16.74 s) and MPreterm (1042 ± 15.40 s) (p = 0.643), the immobile time in the border zone was significantly lower in MTerm (37.09 ± 6.89 s) than that of MPreterm (59.25 ± 6.495 s) mice (p = 0.0361). There was no significant difference in total distance traveled between the two groups (p = 0.120) (Figure 1(d)).

However, our DCE-MRI analysis revealed that BBB permeability in twelve of the fourteen brain regions we examined during early development was significantly lower in MTerm mice. The DCE-MRI statistic K^trans is the volume transfer constant (min⁻¹) for the gadolinium-based contrast agent between blood plasma and the tissue extravascular extracellular space (EES), reflecting the combined processes (primarily blood flow and capillary leakage) that determine gadolinium influx rate from plasma into the EES. Higher values of K^trans indicate greater vascular permeability. Based on the curves used to calculate K^trans values (Figure 2(d)), we detected significant functional deficits in BBB permeability in MPreterm mice compared to MTerm mice at two weeks of age (n_MTerm = 8, n_MPreterm = 6). These deficits were observed in the medulla, midbrain, cortex, thalamus, olfactory bulb, white matter, cerebellum, pons, hippocampus, striatum, corpus callosum, and ventricles, though not the hypothalamus or pallidum. The thalamus and white matter, but no other brain regions analyzed, continued to show significantly lower K^trans values in three-week-old MPreterm mice (n_MTerm = 6, n_MPreterm = 9), demonstrating that microbiome-induced developmental BBB disruption occurs during early life (Figure 2(e)).

To uncover broad processes that may participate in the effects of the microbiome on neurodevelopment, we performed gene set enrichment analysis (GSEA) on our scRNA-seq differential gene expression data. GSEA uncovered a complex reprogramming of the BBB following microbiome perturbation, with each cell type showing both unique and coordinated responses. In ECs, we observed significant upregulation of immune defense and interferon-beta signaling pathways in the MPreterm group. Mural cells from MPreterm mice showed significant enrichment in cilium assembly, tRNA methylation, and wound healing processes. Astrocytes from the MPreterm group demonstrated enrichment in cellular structural pathways and metabolic functions compared to those from MTerm mice. Fibroblasts displayed upregulation of axoneme assembly and immune functions such as lymphocyte-mediated immunity in MPreterm mice. Notably, across all cell types examined, synaptic signaling-related pathways showed strong negative enrichment scores, indicating consistent and robust downregulation in the MPreterm group. These synaptic pathways, including trans-synaptic signaling, chemical synaptic transmission, and modulation of chemical synaptic transmission, displayed some of the most significant enrichment patterns in the dataset (Figure 3(c), Table S5).

As downregulation of synaptic signaling represents a plausible mechanism of cognitive impairment, we decided to further explore the genes involved in this GSEA finding. We focused specifically on differential expression in astrocytes and ECs, as these cell types showed the largest transcriptomic differences between MTerm and MPreterm mice. We examined DEGs with a Gene Ontology (GO) annotation of "synaptic signaling" (GO:0099536) using the Jackson

Laboratory Mouse Genome Informatics database and found 82 such DEGs in astrocytes and 24 in ECs in our transcriptomic data. Here, we quantified the level of differential expression by log₂FC, with negative values indicating downregulation and positive values indicating upregulation in MPreterm animals. In astrocytes, we identified striking downregulation in MPreterm compared to MTerm mice of Nrgn (log₂FC = −2.11, p_adjusted = 1.06 × 10⁻³⁹), Camk2b (log₂FC = −1.80, p_adjusted = 2.13 × 10⁻⁵), and Snca (log₂FC2 = −1.74, p_adjusted = 8.13 × 10⁻¹¹), all of which are genes vital in the control of various synaptic functions.^32–34 Other key genes downregulated in MPreterm astrocytes with "synaptic signaling" annotations include Rab3a (log₂FC = −0.81, p_adjusted = 2.42 × 10⁻¹⁵) and Syt1 (log₂FC = −2.11, p_adjusted = 1.69 × 10⁻¹¹). Both genes are involved in Ca²⁺-mediated vesicle exocytosis, meaning that astrocyte-neuron communication may be affected (Figure 3(d), Table S6).^35–38

We further observed downregulation of several "synaptic signaling"-annotated genes in MPreterm ECs that also play roles in BBB integrity. Rab3a (log₂FC = −0.70, p_adjusted = 1.50 × 10⁻⁴) and Nrgn (log₂FC = −1.41, p_adjusted = 8.98 × 10⁻¹¹) were both downregulated in MPreterm ECs. Rab3a is a regulator of endothelial lumen and tube formation, while Nrgn plays a critical role in controlling BBB permeability through Akt signaling.^39–41Cdh2 (log₂FC = −0.42, p_adjusted = 3.99 × 10⁻⁴) and Trio (log₂FC = −0.38, p_adjusted = 6.56 × 10⁻⁴) were also "synaptic signaling"-annotated genes downregulated in MPreterm ECs. The Cadherin-2 and Trio proteins, likewise, contribute to BBB integrity by forming a complex that promotes the assembly and stabilization of VE-cadherin junctions between ECs.⁴² Both genes are also well-established regulators of synaptic transmission (Figure 3(e–f), Table S6).^43,44 To validate these scRNA-seq findings, we performed an immunofluorescence assay using MTerm and MPreterm brain tissue, co-staining for CD31 (an endothelial cell marker) and Cadherin-2 (n_MTerm = 6, n_MPreterm = 3). As expected, we found no significant difference in CD31 expression but significantly (p < 0.05) higher Cadherin-2 expression at the BBB of MTerm mice compared to that of MPreterm mice (Figure 3(g–h), Figure S5).

Another crucial aspect of BBB function comes from its expression of transporters. Underexpression of transporters can impair the delivery of essential nutrients and removal of waste products, while overexpression can compromise the selectivity of the BBB. Many of these crucial transporters belong to the solute carrier (SLC) superfamily, which encompasses a diverse group of membrane proteins responsible for facilitating the movement of amino acids, neurotransmitters, metabolites, and many other solutes across cellular membranes.⁴⁵ We observed transcriptional dysregulation of SLCs in the MPreterm endothelia. Slc25a25 (log₂FC = −1.17, p_adjusted = 0.00901) and Slc38a2 (log₂FC = −0.493, p_adjusted = 6.36 × 10⁻⁸) were both downregulated in MPreterm ECs. Slc25a25 is a Ca²⁺-sensitive mitochondrial carrier that maintains ATP homeostasis by mediating the exchange of ATP- Mg²⁺ for inorganic phosphate across the inner mitochondrial membrane, and Slc38a2 facilitates the net uptake of small neutral amino acids into the cell in an Na⁺-dependent manner.^46,47 Conversely, Slc29a4 (log₂FC =0.651, p_adjusted = 0.00915) and Slc16a2 (log₂FC =0.459, p_adjusted = 1.23 × 10⁻⁷) were both upregulated in MTerm ECs. Slc29a4 is a transporter of dopamine, serotonin, and other monoamine neurotransmitters, while Slc16a2 mediates the movement of thyroid hormones (triiodothyronine and thyroxine) across cell membranes and is known to function at the BBB (Table S4).^48–51

The microbiome similarly drove substantial transcriptional reprogramming in astrocyte populations across the brain. Subclustering analysis revealed five functionally distinct astrocyte transcriptional states: insult-mediating astrocytes (IMAs), synaptic connectivity astrocytes (SCAs), synaptic signaling astrocytes (SSAs), adhesion-enriched astrocytes (AEAs), and AMPA

receptor (AMPAR)-expressing astrocytes (ARAs) (Figure 4(d)). The five astrocyte populations revealed distinct functional identities through their gene expression profiles. IMAs were characterized by the expression of developmental regulator Ascl1 and extracellular matrix protein Thbs4, which are both well-known to be upregulated in response to CNS injury.^57,58 SCAs expressed Erbb4, a receptor tyrosine kinase, and scaffolding protein Sntg1.^59,60 SSAs were distinguished by markers Camk2a, which encodes the α subunit of calmodulin-dependent kinase II, and Cnih2, which regulates AMPAR assembly, highlighting their role in synaptic function.^61–63 Finally, AEAs displayed adhesion molecules Mamdc2 and Cntnap4, while ARAs had marker genes Gria4, which encodes the AMPAR subunit GluR4, and Cacng5, which encodes an auxiliary AMPAR subunit and modifies its activity (Figure 4(e)).^64–68 Logistic regression analysis revealed that IMAs are significantly enriched in MPreterm mice compared to controls (OR = 1.88, p_adjusted = 7.98 × 10⁻¹³), while SSAs are significantly more abundant in MTerm mice (OR = 0.35, p_adjusted < 10⁻¹⁶) (Figure 4(f)).

Gene set enrichment analysis of endothelial and astrocyte subclusters revealed distinctive functional adaptations in response to microbiome immaturity, providing further resolution to analyze perturbed pathways. In MPreterm mice, BR-CECs exhibited enrichment of immune-related pathways, including interferon-beta signaling and MHC class I-mediated antigen presentation. MPreterm VM-CECs, on the other hand, were enriched for processes related to mitosis and chromosomal dynamics like spindle checkpoint signaling and chromosome separation but downregulation of processes related to de novo protein folding and development. Both AVECs and BR-CECs exhibited downregulation of processes related to synaptic transmission in MPreterm mice, while VM-CECs did not. All MPreterm astrocyte subclusters, though, exhibited notable downregulation of a diversity of synaptic signaling- related processes including regulation of vesicle fusion, regulation of neurotransmitter transport, and synaptic vesicle exocytosis. MPreterm SSAs were also enriched for terms related to transcription and mRNA processing (Figure 4(g–h), Table S7).

This study was conducted with the approval of the Institutional Review Board of the University of Chicago (Protocol ID: IRB210270), which follows the Declaration of Helsinki and the Belmont Report. The parents of the infants enrolled in the study gave informed consent. The fecal sample donors were a full-term, healthy infant at GA of 39 weeks born via vaginal delivery (MTerm) and a preterm infant at GA of 24 weeks born via cesarean section (MPreterm). Patient post-meconium fecal samples provided at three weeks postnatal were obtained from the Center for the Science of Early Trajectories biorepository at the University of Chicago Comer Children's Hospital, which operates a level IV NICU in Chicago, IL, USA. The infants enrolled had no congenital anomalies.

This study was conducted with the approval of the Institutional Animal Care and Use Committee of the University of Chicago. All animal experiments complied with the Guide for the Care and Use of Laboratory Animals (National Research Council (USA), 2011). All mice in this study were allowed ad libitum access to NIH-31 germ-free chow and water. Germ-free C57BL/6J mice (Jackson Laboratory, Bar Harbor, ME, USA) were maintained at the Gnotobiotic Research Animal Facility at the University of Chicago, where animals are routinely tested for microbes and parasites to ensure they remain germ-free. A mixture of female and male mice was used for all analyses. Sex differences in mouse behavior are strain-dependent and were not statistically significant in this study (Table S2).^130–134 Exact sample sizes and breakdowns by treatment, age, and sex for each experiment are available in Table S10. Information on the use of mice from different litters is available in Table S11.

The microbiota transplantation protocol used for this study was followed as previously described.¹³⁵ 100 mg aliquots of the fecal samples from each donor were thawed. Under anaerobic conditions, fecal samples were resuspended in PBS, homogenates were clarified through 100 µm nylon filters and stored in PBS containing glycerol (concentration of 15% v/v) at −80°C in 600 µL aliquots. Germ-free mice were colonized by gavaging eight- to nine-week-old pregnant females (estimated embryonic day 15–17) with 250 µL aliquots of the prepared MTerm or MPreterm fecal supernatant. Pups in both the MTerm and MPreterm groups were delivered spontaneously and full-term and kept with their mothers until weaning at three weeks of age, allowing them to naturally acquire the transplanted patient-derived microbiota.

Behavioral tests were conducted with both MTerm and MPreterm mice at four weeks of age. Trials were recorded and analyzed with ANY-Maze software (Stoelting Co., Wood Dale, IL, USA).

The same procedure previously described for the contextual and cued fear conditioning test was employed here.⁷² Mice underwent fear conditioning involving pairings of a CS of a 50 dB white noise for 30 s and an unconditioned stimulus (US) of a 0.3 mA foot shock during the final 2 s of CS. On the conditioning day, mice freely explored the chamber for 2 min before experiencing three CS-US pairings at 2 min, 4 min, and 6 min, then remained undisturbed for an additional 90 s. After 24 h, mice were tested for contextual fear by being returned to the original chamber for 5 min without stimuli. 2 h later, cued fear was tested by placing mice in a novel environment for 6 min, with the CS presented during the last 3 min. Fear memory was quantified by measuring freezing duration during the contextual test and the final 3 min of the cued test.

Mice were placed in a circular tank (120 cm in diameter) filled with room-temperature water (22 °C) that had been made opaque by adding nontoxic white tempera paint. High-contrast black and white images were positioned on curtains surrounding the tank to provide visual reference points. During the training phase, mice learned to escape onto a visible 10 cm diameter platform that extended 1 cm above the water surface. This phase consisted of five trials, with the platform location changed between each trial. In the testing phase, mice were challenged to locate a hidden platform submerged 1 cm below the water surface in the southeast quadrant. Testing occurred over four consecutive days, with sessions conducted 24 h apart. Performance was measured by recording the time it took mice to find the escape platform. If a mouse failed to reach the platform within 60 s, a maximum score of 60 s was recorded.

Mice were placed individually in the center of an open field box. Their following measures of spontaneous motor activity were recorded: time spent immobile in the center and border zones, total time in the center and border zones, and distance traveled.

All imaging experiments were performed at the University of Chicago Integrated Small Animal Imaging Research Resource. MTerm and MPreterm mice of two and three weeks of age were used. Mice underwent anesthesia before the experiments and were maintained on 1.5–2.5% isoflurane. Temperature, as well as heart and respiration rates, were monitored using a fiber optic detection system (SA Instruments, Stony Brook, NY, USA). A 9.4 T small animal scanner equipped with 11.6 cm inner diameter actively shielded gradient coils (maximum constant gradient strength across all axes: 230 mT/m) was used. Each mouse was positioned supine on an animal holder and placed into a 30 mm diameter quadrature volume coil (Rapid MR International, Columbus, OH, USA).

Whole-brain coverage was achieved using multi-slice spin echo T2-w imaging in the coronal direction with a Rapid Acquisition with Relaxation Enhancement (RARE) pulse sequence (repetition time (TR) = 4000 ms, echo time (TE) effective = 20 ms, field-of-view (FOV) = 25.6 × 19.2 mm², matrix size = 256 × 192, slice thickness = 0.5 mm, RARE factor = 4, number of excitations (NEX) = 2). Brain regional volume analysis utilized a 3D mouse brain atlas from the VivoQuant (Perceptive, Nottingham, UK) software. This 3D brain atlas is based on the Paxinos–Franklin atlas registered to high-resolution MRI with 100-μm near-isotropic data.^136,137 This yielded 14 brain regions of interest (ROIs): medulla, cerebellum, midbrain, pons, cortex, thalamus, hypothalamus, hippocampus, striatum, pallidum, olfactory bulbs, corpus callosum, white matter, and ventricles.

A RARE-VTR (variable TR) pulse sequence was used (TR = 281, 350, 500, 1000, 1500, 2000, 3000, 5000, 10000 ms, TE = 12.3 ms, RARE factor = 4, FOV = 25.6 × 19.2 mm², matrix size = 128 × 96, thickness = 1.5 mm, number of slices = 9, NEX = 1).

DCE-MRI was used to model the movement of gadolinium contrast in the brain as a method to measure BBB permeability. T1-w DCE-MRI data (TR/TE = 52.05/1.56 ms, FOV = 25.6 × 19.2 mm², matrix size = 128 × 96, flip angle = 30°, slice thickness = 1.5 mm, number of slices = 9) were acquired with a temporal resolution of 5.0 s. The DCE-MRI data were continually acquired before (for 1 min), during, and after a bolus intravenous injection of ~0.2 mmol/kg (body weight) of Omniscan (generic name: gadodiamide) (GE HealthCare, Chicago, IL, USA) for a total duration of 8 min (96 frames).

DCE-MRI data were analyzed in MATLAB. Specifically, the 14 brain ROIs obtained from T2-w images were applied to the DCE-MRI to determine the permeability of BBB at these regions. The contrast agent concentration curve C(t) as a function of time t was calculated using a previously published method based on the gradient echo signal equation (non-linear model).¹³⁸ The standard Tofts model was used to fit C(t) to obtain maps of physiological parameters of the volume transfer constant K^trans between blood plasma and the extravascular extracellular space. The arterial input function was traced from the internal carotid artery.

Whole-brain samples were collected from MTerm and MPreterm mice at two weeks of age. Before tissue dissociation, the meninges were removed. To isolate brain cells, the Neural Tissue Dissociation Kit – Postnatal Neurons (Miltenyi Biotec, Bergisch Gladbach, Germany) was used. The manufacturer's protocol was followed with minimal modifications. First, a 40 µm strainer was used instead of a 70 µm one. This allowed for more precise filtering out of debris and, thus, better cell isolation. Moreover, the enzymatic tissue digestion incubation times were decreased from 15 to 8 min during the initial dissociation and from 10 to 6 min during the secondary dissociation to reduce potential oxidative stress, preserve RNA quality, and increase cell viability.

scRNA-seq libraries were prepared using the Chromium Next GEM Single Cell 3' v3.1 reagent kit (10x Genomics, Pleasanton, CA, USA). Cells were barcoded using a 10x Genomics Chromium Controller. Eleven polymerase chain reaction cycles were applied to generate cDNA, with thirteen cycles for final library generation. The final scRNA-seq libraries were sequenced with the Illumina NovaSeq X (Illumina, San Diego, CA, USA).

After receiving the raw FASTQ sequencing files, Cell Ranger's (10x Genomics) count function was used to align sequencing reads to the GRCm38 Mus musculus reference genome. Count matrices were imported to R using the Seurat (v5.1.0) package's Read10X function, and Seurat objects were created with sample-specific metadata.¹³⁹ Cells were filtered based on quality control (QC) metrics (nCount_RNA > 800, nFeature_RNA > 200, mitochondrial percentage < 10%), data normalized with NormalizeData, and variable features identified with FindVariableFeatures. The data were scaled with ScaleData and dimensionality reduction was performed with RunPCA and RunUMAP. To remove doublets, DoubletFinder (v2.0.4) was used, with an expected doublet rate of 0.8% per 1,000 cells.¹⁴⁰

To integrate data from multiple samples and minimize batch effects, batch correction was performed using the Harmony (v1.2.1) package with the "sample" variable as the integration factor.¹⁴¹ This approach preserves the underlying biological variation in gene expression within each cell type across treatment conditions and assumes that treatment condition does not drastically alter broad cell type identity (EC, astrocyte, mural cell, or fibroblast), while still allowing for large differences in subcluster composition due to the downstream analysis. Next, UMAP visualization and clustering on the Harmony-corrected embeddings (RunUMAP, FindNeighbors with dims = 1:50, FindClusters with resolution = 0.2) was performed. Clusters were annotated for cell types using canonical marker genes from PanglaoDB and CellMarker 2.0.^142,143 BBB-relevant cell types were identified by expression of the following genes: Cldn5, Flt1, and Cd93 for ECs; Gfap, Aqp4, and Sox9 for astrocytes; Rgs5, Pdgfrb, and Kcnj8 for mural cells; Col1a1, Col1a2, and Lum for fibroblasts. Marker genes for all other cell types can be found in Figure S4. For clearer visualization for the UMAP plot in Figure 3(a), the data were subset to include only BBB-relevant clusters (ECs, astrocytes, mural cells, and fibroblasts) and re-embedded using PCA (50 dimensions). The colorblind-friendly ggpubfigs (v0.0.1) palette was used for visualization of BBB-relevant cell types and some of the other cell types in Figure 3 and Figure S4.¹⁴⁴

Subclustering analyses were then performed on ECs and astrocytes (RunPCA, FindNeighbors, FindClusters with resolution = 0.1 for ECs, 0.05 for astrocytes) to identify distinct transcriptional states within these populations. FindMarkers with the Model-based Analysis of Single-cell Transcriptomics (MAST) test was used to identify subcluster marker genes, comparing each individual cluster against all other cells within the same main cell type.¹⁴⁵ The following parameters were used: latent variables of sex, mitochondrial percentage, and ribosomal percentage, a minimum expression threshold of 0.1 (expressed in at least 10% of cells), p_adjusted < 0.05, and log₂FC ≥0.585 (at least a 50% difference in average expression compared to other clusters). These subclusters were named according to their marker genes (see Results). To assess subcluster enrichment between conditions, logistic regression was performed, controlling for sex.

For differential expression analysis between MTerm and MPreterm conditions for both overall cell types and subclusters, FindMarkers was used with the MAST test, with sex, mitochondrial percentage, and ribosomal percentage set as latent variables. Genes with adjusted p-value < 0.05 and |log₂FC| ≥ 0.3785 (at least a 30% difference in average expression between conditions) were considered to be significantly differentially expressed. A list of all genes with a GO annotation of "synaptic signaling" (GO:0099536) were retrieved from the Jackson Laboratory Mouse Genome Informatics database to identify which DEGs were implicated in synaptic signaling.^146–148

Gene set enrichment analysis was performed on the expression data, comparing between MTerm and MPreterm conditions. As the input for GSEA, the complete FindMarkers output (not only the DEGs) was used, creating ranked gene lists based on the log₂FC values for each cell type and subtype. Gene symbols were converted to ENTREZ IDs using the org.Mm.eg.db (v3.19.1) mouse genome database package and gseGO from the clusterProfiler (v4.12.6) package was run with the following parameters: ont = "BP" (Biological Process), minGSSize = 10, and maxGSSize = 500.^149,150 This allowed for the identification of enriched biological pathways and processes that were upregulated in either MTerm or MPreterm conditions, as indicated by positive or negative normalized enrichment scores, respectively, with significant pathways defined by p_adjusted < 0.05 (Benjamini–Hochberg procedure).

Mouse brains were freshly collected at two weeks postnatal and embedded in optimal cutting temperature compound before freezing. An 8 µm thick section from each mouse brain was fixed with ice-cold methanol for 20 min at −20°C. Samples underwent permeabilization with 0.2% Triton-X in PBS (PBST) for 15 min, followed by incubation in blocking solution (5% goat serum in PBST) for 1 h at room temperature. Brain sections were then incubated overnight with 50 μL of primary antibody solution at 4°C. After four PBST washes of 10 min each, sections were incubated with fluorophore-conjugated secondary antibodies for 1 h at room temperature. Primary antibodies used were rat anti-CD31 (cat#NB600–1475, Novus Biologicals, Centennial, CO, USA) and mouse anti-Cadherin-2 (cat#ab98952, Abcam, Cambridge, UK). Secondary antibodies used were goat anti-rat IgG, Alexa Fluor 488 (cat#A-11006, Thermo Fisher Scientific, Waltham, MA, USA) and goat anti-mouse IgG, Alexa Fluor Plus 647 (cat#A32728, Thermo Fisher). Nuclei were counterstained using ProLong Gold Antifade Mountant with DNA Stain DAPI (cat#P36935, Thermo Fisher). A STELLARIS confocal microscope (Leica Microsystems, Wetzlar, Germany) was used to capture five to seven images from each mouse brain. Image processing and analysis were performed using ImageJ.¹⁵¹ The BBB was identified via expression of CD31. CD31+ polygonal areas were outlined (Figure S5) and used for comparisons of fluorescence in Figure 3.

The data collection approach for the metabolomics experiment was followed as previously described.⁷² Serum metabolites were extracted by adding 100% methanol containing heavy-labeled internal standards (stored at −80 °C) to serum samples at a 1:10 ratio (100 µL serum:1,000 µL solvent). Samples were vortexed for 1 min, incubated overnight at −20°C, then centrifuged at −10°C, 20,000 × g for 15 min. Following centrifugation, 100 µL of supernatant was dried in a Genevac EZ-2 Elite (ATS Scientific Products, Warminster, PA, USA). Samples were resuspended in 200 µL of a 50/50 water/methanol solution, mixed at 4°C, 1,000 rpm for 15 min in an Eppendorf (Hamburg, Germany) ThermoMixer C, and centrifuged at 4°C, 20,000 × g for 15 min, after which 100 µL of clarified supernatant was transferred to MS vials. Analysis was performed with a Thermo Fisher liquid chromatography system coupled to an Orbitrap IQ-X mass spectrometer (Thermo Fisher) in positive mode, injecting 2 µL of sample onto a Cortecs UPLC T3 column (1.2 µm, 2.1 × 100 mm) (Waters Corporation, Milford, MA, USA) at 30°C, using water with 0.1% formic acid (mobile phase A) and 95% acetonitrile with 0.1% formic acid (mobile phase B). The gradient (0% B for 1 min, linear ramp to 100% B over 10 min, hold 2 min, re-equilibrate 2 min) flowed at 0.35 mL/min. Electrospray ionization parameters included 3.5 kV spray voltage, 400°C vaporizer temperature, and 133–2000 m/z detection window. MS2 scans employed dynamic exclusion after three selections of the same precursor ion within 15 s (exclusion duration 10 s), with an isolation window of 0.7 m/z (without offset) and fixed normalized collision energy of 35%.

The metabolomics data analysis pipeline was applied as previously described.⁷² Feature lists were generated from raw data using MZmine (v2.53), beginning with mass detection set to noise cutoff thresholds of 8.0 × 10³ for MS2 and 2.0 × 10⁴ for MS1 scans.¹⁵² Extracted ion chromatograms were built using the ADAP chromatogram builder (minimum group size: 5 scans; group intensity threshold: 2.0 × 10⁴; minimum highest intensity: 4.0 × 10⁴; m/z tolerance: 0.01 Da or 3 ppm). Chromatograms were deconvoluted using the ADAP wavelets algorithm (S/N threshold: 10; minimum feature height: 6 × 10⁴; coefficient/area threshold: 110; peak duration: 0–1 min; RT wavelet range: 0.01–0.10 min). Isotopic peaks were grouped (maximum charge: 3; representative isotope: lowest m/z; tolerance: 0.01 Da or 3 ppm m/z, 0.1 min retention time (RT)), and features across samples were aligned (join aligner: 0.01 Da or 5 ppm m/z, 0.1 min RT). Duplicate features were filtered out, and missing peaks were gap-filled (peak-finder algorithm; 10% intensity tolerance, 0.01 Da or 5 ppm m/z, 0.1 min RT). Feature lists were exported for Global Natural Products Social Molecular Networking (GNPS) analysis.¹⁵³ Feature-Based Molecular Networking¹⁵⁴ via GNPS was used to filter MS/MS data (precursor and fragment tolerances: 0.01 Da; cosine score cutoff: 0.7; ≥ 4 matching peaks; mutual top-10 node connections; maximum family size: 100). Analog searching was set at a precursor mass difference ≤100 Da with identical filtering for library matches. GNPS-derived feature lists were further filtered in Excel based on blanks and QC injections to eliminate low-quality data. Peak areas within samples were normalized against the average area of internal standards (D₄-cholic acid, D₄-taurocholic acid, D₄-glycocholic acid, D₄-deoxycholic acid, D₄-chenodeoxycholic acid,¹³C₉-tyrosine, ¹³C₉-phenylalanine, and ¹³C₁₁-tryptophan). Three QC injections, using pooled samples, were used to assess reproducibility by excluding features with coefficient of variation > 20%. This resulted in a final list of 1,571 features; after removing known contaminants and polyether polymers, 231 putatively identified features remained. MetaboAnalyst was utilized to perform statistical analysis and visualization for both full and putatively identified feature lists, removing only features with single or constant values, without additional filtering or normalization.¹⁵⁵

Fecal samples were collected from mice at four weeks postnatal and preserved at −80°C. Genomic DNA extraction and shotgun metagenomic sequencing were performed at the Duchossois Family Institute (DFI) Microbiome Metagenomics Facility using the Illumina NextSeq platform. Paired-end libraries were created in batches of 60–72 samples, with insert sizes of approximately 350 bp per sample. High-throughput sequencing generated about seven to eight million paired-end reads per sample with 150-bp read lengths.

The raw sequencing data were processed by first trimming adapters off from raw reads and performing quality assessment and control with Trimmomatic (v0.39).¹⁵⁶ Human genomic sequences were identified and removed using KneadData (v0.7.10).¹⁵⁷ Taxonomic classification of the resultant high-quality reads was performed using MetaPhlAn (v4.0.2), and relative abundances were tabulated.¹⁵⁸ Host-removed sequencing reads were also processed using HUMAnN (v3.8), and pathway abundances were normalized to relative abundances using the built-in normalization option.¹⁵⁷

These 16S rRNA-seq protocols were followed as previously described.^69,72,159 Fecal samples from the human infant donors were stored at −80°C and subsequently submitted for genomic DNA extraction and subsequent 16S rRNA-seq via the GS-FLX platform with XLR70 sequencing chemistry (Roche, Basel, Switzerland) at the Institute for Genomics and Systems Biology at Argonne National Laboratory.^160,161 PCR primers targeting the V3–V4 regions of the 16S rRNA gene (nucleotide positions 338–802 in Escherichia coli) were used. These primers were constructed with 454 sequencing-compatible adapter sequences and unique 8-nucleotide barcode identifiers. The resulting sequence data were processed and merged with DADA2 in QIIME2 (v2019.10).^162,163 Taxonomic assignment used the QIIME2 naïve Bayes classifier trained on the SILVA database (release 138).¹⁶⁴

Mouse fecal samples were collected at two weeks postnatal, four weeks postnatal, and from dams. These samples were stored at −80°C. They were then submitted to the DFI Microbiome Metagenomics Facility for genomic DNA extraction before 16S rRNA-seq using the Illumina MiSeq. Sequencing reads were processed in R using DADA2 (v1.18.0) with minor modifications.¹⁶² Forward and reverse reads were truncated to 190 bp to remove low-quality tails, and chimeric sequences were discarded via DADA2's consensus algorithm. Amplicon sequencing variants between 320 bp and 365 bp were retained as high-quality. Taxonomy was assigned with the RDP classifier (v2.13, training set 18/release 11.5) with 80 as the minimum bootstrap confidence score.¹⁶⁵

For both the fear conditioning and open field tests, the Shapiro-Wilk test was used to determine whether behavioral outcomes were normally distributed at α = 0.05 for both MTerm and MPreterm mice. If a given behavioral outcome was determined to be non-normally distributed for at least one of the treatment groups, then the data was log-transformed for hypothesis testing. If not, the data remained unchanged. Then, Welch's t-test was performed. For the Morris water maze, a two-way ANOVA was performed with treatment and training day as predictors. The same statistical tests were performed to test for sex differences as were used to test for differences due to treatment but with Benjamini–Hochberg adjustment of the p-values for the fear conditioning and open field test readouts to account for the same outcome being queried twice (once in the MTerm comparison and once in the MPreterm comparison).

For determination of DEGs in scRNA-seq data, FindMarkers with MAST p-value adjustment was used.¹³⁹ Gene set enrichment analysis was performed with gseGO, with the Benjamini–Hochberg procedure to adjust p-values.¹⁵⁰ For comparison of subclusters, a logistic regression model was fit to the subclustering data to retrieve odds ratios. The corresponding p-values were calculated via the Wald test and Bonferroni-adjusted. Wilcoxon rank sum tests were used for comparison of fluorescence in immunostaining images.

In comparisons of brain regional volume and BBB permeability, the Shapiro–Wilk test was applied to assess whether measurements were normally distributed at α = 0.05 for both treatment groups. If, for a given brain region, the data were non-normally distributed for at least one of the MTerm or MPreterm groups, then the data were log-transformed for hypothesis testing. The data remained unchanged otherwise. Welch's t-test was subsequently performed.

For the identification of microbial metabolites present at altered levels in MTerm and MPreterm mice, unpaired t-tests were used before adjustment by the Benjamini, Krieger, and Yekutieli procedure.¹⁶⁶

In the analysis of shotgun metagenomic data, the Wilcoxon rank sum test with Benjamini–Hochberg adjustment was utilized to compare α-diversity and relative abundances of bacterial taxa and genomic pathways. For further comparison between MTerm and MPreterm gut microbiota, a Bray-Curtis dissimilarity matrix was computed using vegdist from the vegan package.¹⁶⁷ PERMANOVA (adonis2) with 10⁵ permutations was used, with treatment as the independent variable.

For 16S rRNA-seq analysis, the modified Wilcoxon rank sum test for zero-inflated data¹⁶⁸ was used for comparisons of relative abundance of bacterial taxa, and p-values were adjusted via the Benjamini–Hochberg procedure. To assess β-diversity, a Bray-Curtis dissimilarity matrix was constructed, and PERMANOVAs were performed with 10⁵ permutations, with either treatment or age as the independent variable.