Differential Methylation Profile in Fragile X Syndrome-Prone Offspring Mice after in Utero Exposure to Lactobacillus Reuteri

Environmental factors such as diet, gut microbiota, and infections have proven to have a significant role in epigenetic modifications. It is known that epigenetic modifications may cause behavioral and neuronal changes observed in neurodevelopmental disabilities, including fragile X syndrome (FXS) and autism (ASD). Probiotics are live microorganisms that provide health benefits when consumed, and in some cases are shown to decrease the chance of developing neurological disorders. Here, we examined the epigenetic outcomes in offspring mice after feeding of a probiotic organism, Lactobacillus reuteri (L. reuteri), to pregnant mother animals. In this study, we tested a cohort of Western diet-fed descendant mice exhibiting a high frequency of behavioral features and lower FMRP protein expression similar to what is observed in FXS in humans (described in a companion manuscript in this same GENES special topic issue). By investigating 17,735 CpG sites spanning the whole mouse genome, we characterized the epigenetic profile in two cohorts of mice descended from mothers treated and non-treated with L. reuteri to determine the effect of prenatal probiotic exposure on the prevention of FXS-like symptoms. We found several genes involved in different neurological pathways being differentially methylated (p ≤ 0.05) between the cohorts. Among the key functions, synaptogenesis, neurogenesis, synaptic modulation, synaptic transmission, reelin signaling pathway, promotion of specification and maturation of neurons, and long-term potentiation were observed. The results of this study are relevant as they could lead to a better understanding of the pathways involved in these disorders, to novel therapeutics approaches, and to the identification of potential biomarkers for early detection of these conditions.


Introduction
Several neurodevelopmental disabilities such as ASD cannot be well-characterized by the one-disease one-gene school. Their complex nature is attributed to the interplay between genetic, epigenetic alterations, neuroinflammation, and the innate immune system, that are associated with factors such as alteration in the gut microbiota and other environmental factors [1]. In the past decade or so, the interest of the research community started to increase towards epigenetic modifications such as DNA methylation, histone modifications, RNA interference, and their influence on gene expression [2]. Several studies have shown how different epigenetic modifications may cause behavioral and neuronal changes observed in neurodevelopmental disabilities [3][4][5][6]. Furthermore, it was shown

DNA Methylation Array Processing and Mapping MethylationEPIC Primer to the Mouse Genome
Five hundred ng of genomic DNA was bisulfite-treated using the EZ DNA methylation kit (Zymo research, Irvine, CA, USA). Although Infinium Mouse Methylation BeadChips have been developed recently to measure the methylation profile of mice, there was no mature mouse methylation array when this study started, and thus, the study was carried out using the Illumina Human Infinium Methylation BeadChips. Bisulfite-converted DNA was hybridized with the Illumina Human Infinium MethylationEPIC BeadChip Kit (Illumina, San Diego, CA, USA) and then scanned with the Illumina iScan System (Illumina, CA, USA) using the manufacturer's standard protocol. Samples from the two groups were randomly loaded on the arrays.
The Infinium MethylationEPIC BeadChip Kit included >850,000 individual CpG sites genome-wide at single-nucleotide resolution. As the Infinium MethylationEPIC array was not specifically designed for the mouse genome, many primers for each probe in the array do not target or specifically align within the mouse genome, so their location cannot be identified in the mouse genome with certainty.
Hence, 22,569 out of the 850,000 probes that were identified in previous studies [20,21] to have a unique best alignment score when mapping the probe primers to the mouse genome (GRCm38) were used for the following analysis.

Methylation Data Preprocessing
DNA methylation data were pre-processed and quantile-normalized with R package 'minfi' (Version 4.2.1, Bioconductor, Online Platform) [22]. Probes were filtered according to the following criteria: probes with a detection p-value greater than 0.01 in more than 5% of total samples (≥2).

Identification of Differentially Methylated CpGs (DMs)
The Limma package was used to compare the methylation difference of CpGs that passed the quality control described above between the two groups. M value was used in statistical analysis instead of β value. M value is the ratio of the methylated to the unmethylated intensity and has a statistical advantage over the β value [23].

Pathway Analysis
Before the pathway analysis, the significance of each gene was determined. In the case of genes with one CpG site, the p-value of the CpG site was used as that of the gene. The p-values of genes with multiple CpG sites in the gene were calculated by combining the p-values of all CpG sites in the gene, using Fisher's method [24]. R package 'fgsea' was used to carry out gene set enrichment analyses on gene sets based on Gene Ontology and metabolic pathways [25]. In Gene Ontology, gene sets are divided into three categories: biological process (BP), molecular function (MF), and cellular component (CC). The metabolic pathway was collected from four different sources: EHMN, HUMANCYC, KEGG, and REACTOME.

Constructing a Functional Association Network for Some of the Proteins Involved in Brain Functions Using String
Different proteins that are encoded by the genes involved in different brain functions based on the conducted pathway analysis were inserted into the String database to see the possible protein-protein interactions between them [26]. A set of proteins were entered in Genes 2022, 13, 1300 4 of 14 the String database, and an identifier mapping on the input was performed and displayed as a network covering all the mapped proteins and interconnections between them, if any.

Western Blot Analysis
Lyophilized brain tissues from Group A of treated mice (n = 8) and Group B of untreated mice (n = 8) were homogenized with a tip-sonicator (Q-Sonica) in homogenization buffer containing 20 mM of HEPES buffer (Thermo Fisher Scientific, Waltham, MA, USA), 0.25 M of sucrose (MiliporeSigma, Burlington, MA, USA), and 1X protease inhibitor cocktail (MiliporeSigma, Burlington, MA, USA). Homogenates were centrifuged at 200,000× g rpm for 45 min at 4 • C and cytoplasmic fractions (supernatant) were collected and further quantified using the Bradford assay (BioRad Laboratories, Inc. Hercules, CA, USA). Ten µg of protein were loaded on 4-12% Bis-Tris gels (BioRad Laboratories, Inc. Hercules, CA, USA) and run at 100 V for 60 min and then 130 V for 60 min. Resolved proteins were then transferred onto nitrocellulose membranes using the Trans-Blot Turbo transfer system (BioRad Laboratories, Inc. Hercules, CA, USA) at 25 V, 1.0 A, for 30 min. Membranes were stained with Ponceau to test for transfer efficiency, blocked with 3% milk for 1 h at room temperature, followed by incubation with 1:500 diluted Agap3 primary antibodies (sc-390362, Santa Cruz Biotechnology, Inc. Dallas, TX, USA), 1:2000 diluted Shank3 primary antibodies (ab 264347, abcam, Cambridge, MA, USA), and 2.7:1000 diluted Dlg2 primary antibodies (ab 261634, abcam, Cambridge, MA, USA), respectively, overnight at 4 • C. Membranes were then washed in 1X-TBS and incubated with HRP-linked secondary antibody diluted 1:10,000 (Catalog # 1706516, BioRad Laboratories, Inc., Hercules, CA, USA) for Agap3, and HRP-linked secondary antibody diluted 1:2000 (cell signaling, Catalog 7074) for Shank3 and Dlg2, for 1 h at room temperature. Bands were then visualized using a Chemiluminescent substrate, Super Signal West Dura (Thermo Fisher Scientific, Waltham, MA, USA). Densitometry analysis of bands for relative protein quantification was performed using the α Innotech Gel Imaging System (Cambridge Scientific, Watertown, MA, USA).

Mouse Models
Group A mice were male progeny that received probiotic L. reuteri 6475 exposures in utero, and they were clinically normal. Group B mice were male progeny whose mothers received placebo water prenatally. Group B mice had FXS-like symptoms including physical features such as dimorphic head, ears, and enlarged testicles, and behavioral features such as head bobbing, hyperactivity, and stereotypy, similarly to the typical features observed in humans with FXS.

Identification of Perfectly Mapping Probes to the Mouse Genome
Although 128,259 probes passed QC, we only included 17,735 of them with primers that had high mapping quality to the mouse genome for the following analysis (Supplementary Figure S1). The β value distribution of all probes that passed QC has a large peak around 0.25 (Supplementary Figure S2a), while the β value distribution of probes that passed quality control and had a unique best alignment score has two peaks, one around 0 and the other around 0.9 (Supplementary Figure S2b), similar to what was reported by Garcia-Prieto et al. [27]. The probes that passed QC without a unique best alignment score were excluded from the analysis, as they were likely targeting multiple locations, making the measurement of the levels of methylation not reliable.

Identification and Distribution of DMs among the Two Groups
Between the 2 groups, 570 probes were found to be differentially methylated (p-value < 0.05, Supplementary Table S1). The distribution of the top 100 CpGs with the largest variations is shown in Figure 1 ing multiple locations, making the measurement of the levels of methylation not reliabl 3.1.3. Identification and Distribution of DMs among the Two Groups Between the 2 groups, 570 probes were found to be differentially methylated ( value < 0.05, Supplementary Table S1). The distribution of the top 100 CpGs with t largest variations is shown in Figure 1. Interestingly, four samples (20-473, 20-675, 2 676, and 20-678) in Group B showed a very different CpG methylation profile than t other samples in Group B and those in Group A.

Differentially Methylated Genes in FXS-Like Mice Are Involved in Brain Function
DMs were mapped to several genes involved in different neurological pathways ≤ 0.05). Among the key functions, synaptogenesis, neurogenesis, synaptic modulatio synaptic transmission, reelin signaling pathway, promotion of specification and matur tion of neurons, and long-term potentiation were observed. Furthermore, addition

Differentially Methylated Genes in FXS-Like Mice Are Involved in Brain Functions
DMs were mapped to several genes involved in different neurological pathways (p ≤ 0.05). Among the key functions, synaptogenesis, neurogenesis, synaptic modulation, synaptic transmission, reelin signaling pathway, promotion of specification and maturation of neurons, and long-term potentiation were observed. Furthermore, additional genes implicated in schizophrenia, intellectual disabilities, social deficits, repetitive behaviors, deafness, pituitary adenomas, and dysmorphic facies were also identified ( Table 1).

Genes and Pathway Enrichment Analysis
Methylation analysis was obtained from 17,735 sites for the whole mouse genome, which identified DMs located in 6418 genes. The results of the Gene Ontology enrichment analysis datasets shown in Figure 2 indicate that several pathways were significantly enriched between the two groups of mice. As for the biological processes (BP), among the top pathways was that involved in the activation of cysteine-type endopeptidase activity involved in apoptotic processes (GO:0006919). Another pathway was that involved in retrograde transport of membrane-bounded vesicles from endosomes back to the trans-Golgi network, where they are recycled for further rounds of transport (GO:0042147). The cellular component that was found to be of significance was the endocytic vesicle, which is formed by invagination of the plasma membrane around an extracellular substance and fuses with early endosomes to deliver the cargo for further sorting (GO:0030139)). The last subset of Gene Ontology studies was related to molecular function (MF), with the most significant being the ligation activity of acids to amino acids with the concomitant hydrolysis of ATP (GO:0016881).
dase activity involved in apoptotic processes (GO:0006919). Another pathway was that involved in retrograde transport of membrane-bounded vesicles from endosomes back to the trans-Golgi network, where they are recycled for further rounds of transport (GO:0042147). The cellular component that was found to be of significance was the endocytic vesicle, which is formed by invagination of the plasma membrane around an extracellular substance and fuses with early endosomes to deliver the cargo for further sorting (GO:0030139)). The last subset of Gene Ontology studies was related to molecular function (MF), with the most significant being the ligation activity of acids to amino acids with the concomitant hydrolysis of ATP (GO:0016881). In addition, there was an array of statistically significant pathways, including a positive regulation of MAPK activity (p = 0.005), embryonic placenta development (p = 0.009), and the positive regulation of inflammatory response (p = 0.015). A list of GO pathways that were found to be of significance and the corresponding genes involved are listed in Supplementary Table S2. We focused mainly on pathways that are involved in brain development and maintenance of the central nervous system. A total of five pathways were found to be statistically significant, including midbrain development (p = 0.004), positive regulation of neuron apoptotic process (p = 0.01), negative regulation of neuron projection development (p-value = 0.02), cerebral cortex neuron differentiation (p = 0.04), and neuromuscular junction development (p = 0.05) ( Table 2). In addition, there was an array of statistically significant pathways, including a positive regulation of MAPK activity (p = 0.005), embryonic placenta development (p = 0.009), and the positive regulation of inflammatory response (p = 0.015). A list of GO pathways that were found to be of significance and the corresponding genes involved are listed in Supplementary Table S2. We focused mainly on pathways that are involved in brain development and maintenance of the central nervous system. A total of five pathways were found to be statistically significant, including midbrain development (p = 0.004), positive regulation of neuron apoptotic process (p = 0.01), negative regulation of neuron projection development (p-value = 0.02), cerebral cortex neuron differentiation (p = 0.04), and neuromuscular junction development (p = 0.05) ( Table 2).

String-Based Protein-Protein Interaction Network Analysis
Several proteins that were identified through the pathway enrichment analysis to be involved in different brain functions were found to interact together. Some of those proteins were found to interact with several others, such as LIM homeobox transcription factor 1-β (LMX1B) which interacts with Homeobox protein OTX2 (OTX2), Sonic hedgehog protein (SHH), Homeobox protein MSX-1 (MSX1), and Homeobox protein engrailed-1 (EN1). Another protein that was found to interact with other different proteins involved in different functions in the brain is Fibroblast growth factor receptor 2 (FGFR2), which interacts with MSX1, SHH, Wingless-type MMTV integration site family, member 1 (WNT1), and Fibroblast growth factor receptor 1 (FGFR1) (Figure 3). Several proteins that were identified through the pathway enrichment analysis to be involved in different brain functions were found to interact together. Some of those proteins were found to interact with several others, such as LIM homeobox transcription factor 1-β (LMX1B) which interacts with Homeobox protein OTX2 (OTX2), Sonic hedgehog protein (SHH), Homeobox protein MSX-1 (MSX1), and Homeobox protein engrailed-1 (EN1). Another protein that was found to interact with other different proteins involved in different functions in the brain is Fibroblast growth factor receptor 2 (FGFR2), which interacts with MSX1, SHH, Wingless-type MMTV integration site family, member 1 (WNT1), and Fibroblast growth factor receptor 1 (FGFR1) (Figure 3). . Protein-protein interactions network. Protein interactome network for the proteins that were shown to be differentially methylated in FXS-like mouse models compared to controls using the STRING software. All the colored nodes represent the query proteins and first shell of interactions. Each node represents all the proteins produced by a single, protein-coding gene locus. Filled nodes mean that the 3D structure of that protein is known or predicted. The colored lines represent different types of associations. The navy-blue line indicates known interactions from curated databases. The fuchsia lines indicate known interactions that are experimentally determined. The dark green lines represent predicted interactions of neighboring genes, the red lines represent gene fusions, and the royal blue line represents gene co-occurrence. The light green line represents text mining, the black line represents co-expression, and the purple line represents protein homology.

DLG2, SHANK3, and AGAP3 Protein Expression
To validate our methylation findings, among the differentially methylated genes between the two groups, we measured the expression level of a subset of the correspondent encoded proteins, including AGAP3, SHANK3, and DLG2. Western blot analysis in the group of mice treated with probiotics (Group A, n = 8) and in those without any treatment (Group B, n = 8) revealed a significant difference in expression levels of DLG2 between both groups (p = 0.026), while the difference in AGAP3 and SHANK3 expression levels did not reach significance (Figure 4). nodes mean that the 3D structure of that protein is known or predicted. The colored lines represent different types of associations. The navy-blue line indicates known interactions from curated databases. The fuchsia lines indicate known interactions that are experimentally determined. The dark green lines represent predicted interactions of neighboring genes, the red lines represent gene fusions, and the royal blue line represents gene co-occurrence. The light green line represents text mining, the black line represents co-expression, and the purple line represents protein homology.

DLG2, SHANK3, and AGAP3 Protein Expression
To validate our methylation findings, among the differentially methylated genes between the two groups, we measured the expression level of a subset of the correspondent encoded proteins, including AGAP3, SHANK3, and DLG2. Western blot analysis in the group of mice treated with probiotics (Group A, n = 8) and in those without any treatment (Group B, n = 8) revealed a significant difference in expression levels of DLG2 between both groups (p = 0.026), while the difference in AGAP3 and SHANK3 expression levels did not reach significance (Figure 4).

Discussion
During pregnancy and in early childhood, many factors such as hormones, stress, genetics, and diet can affect the brain development of a child. Studies in recent years have demonstrated that gut microbiota and maternal obesity can also influence a child's neurodevelopment [54]. Changes in the microbiota profile that occur within the maternal gut could affect the metabolic profile of the mother and contribute to the immune system and metabolism of their children [55]. Furthermore, microbiota establishment and neurodevelopment were found to be vulnerable to damage during similar developmental windows, with altered developing gut microbiota therefore affecting neurodevelopment and potentially leading to mental health problems [56]. Studies have also

Discussion
During pregnancy and in early childhood, many factors such as hormones, stress, genetics, and diet can affect the brain development of a child. Studies in recent years have demonstrated that gut microbiota and maternal obesity can also influence a child's neurodevelopment [54]. Changes in the microbiota profile that occur within the maternal gut could affect the metabolic profile of the mother and contribute to the immune system and metabolism of their children [55]. Furthermore, microbiota establishment and neurodevelopment were found to be vulnerable to damage during similar developmental windows, with altered developing gut microbiota therefore affecting neurodevelopment and potentially leading to mental health problems [56]. Studies have also shown that mothers who have experienced maternal immune activation (MIA) have an increased risk of their offspring having ASD [57].
Animal experiments indicate that interleukin-6 cytokines (IL-6) are key to the induction of ASD by MIA [58][59][60] and that increased levels of both IL-6 and IL-17a in maternal serum and fetal brains were halted by oral probiotics, as were the loss of parvalbumin-positive neurons and the declines in γ-aminobutyric acid levels in adult offspring brain [61]. The metagenomic whole shotgun sequencing of maternal high-fat diet (MHFD) fecal microbiota showed a significant reduction in the abundance of the commensal L. reuteri [62,63], while synaptic plasticity within the ventral tegmental area, hypothalamic oxytocin, and social behavior in maternal high-fat diet offspring were restored by microbial reconstitution with L. reuteri [64].
In this study, we compared the epigenetic profile of two groups of mice: those that were exposed to L. reuteri in utero and those that were exposed to placebo. A total of 23 genes known to have a role in neurological function (Table 1) showed a significant differential methylation between the two groups. We recognize that we cannot rule out the possibility that observed methylation differences may stem from the genetic variability introduced by our mouse model [65]. However, given the significant differences found across control and experimental groups testing randomly sorted litter mates, our results are likely to be ascribable to an altered gut microbiota.
Among them were the Ncdn, the Agpat1, the Dctn1, the Fbx011, the Emx2, the Dcc gene, and the Gart genes, all encoding for proteins associated with seizures, abnormal hippocampal brain development ID, ASD, ADHD, anxiety, developmental delays, aggression, hyperactivity, altered social skills, delays in speech and language skills, and childhood epilepsy [46,51,66,67]. A number of identified genes were associated with ASD and listed in the SFARI database (Table 1, labeled with an asterisk). A study showed that a mouse model with the Pax6 mutation presented with aggression in the reciprocal social interaction test and hyperactivity signs in the open-field test and the tail suspension test [68]. The Dlgap2 KO mice had deficits in learning, abnormal social behavior, and intense aggressive behavior, which are behavioral features observed in both ASD and FXS [63,64,69].
Among the identified differentially methylated genes, the Dlg2, Shank3, and Agap3 were significantly hypermethylated in the FXS-like group compared to the control group (Table 1).
Dlg2 deficiency in mice contributes to aberrant synaptic transmission accompanied by reduced sociability and increased repetitive behavior [35]. In neurons, specific Dlg2 isoforms bind to NMDA (N-methyl-aspartate) receptors and K channels, mediating their clustering on the postsynaptic membrane, and they are differentially expressed in plasmacytoid dendritic cells (pDCs), potentially having different functions in the neuronal postsynaptic membrane [36]. Furthermore, two very similar proteins, PSD-93 and PSD-95, which developed from a duplication of a single ancestral gene, play an opposite role in the maturation of synapses from immature to active, with PSD-95 promoting maturation, while PSD-93 slows it down [34]. Thus, fine-tuning and a balance between these proteins appear to be essential during developmental critical periods, suggesting that impairment in their expression in either direction could be relevant or constitutes a potential leading cause of neurodevelopmental disorders.
SHANK3, highly expressed in both the peripheral and the central nervous system, is a scaffolding protein, which plays a key role in synaptic development and function. Impaired gene expression, mutations, methylation, and alternative splicing of SHANK3 have been reported in ASD [41]. Moreover, SHANK3-haploinsufficiency results in altered function, which negatively affects synaptic development and function, and its disruption contributes to ASD etiology and is also observed in individuals affected by Phelan-McDermid syndrome [70][71][72]. In this study, we observed hypermethylation, lower protein expression levels, and thus haploinsufficiency of SHANK3 expression in the FXS-like group compared to the controls, which strongly suggest its potential role in the phenotype displayed by these mice.
Centaurin-gamma3, or AGAP3, is an essential signaling component of the NMDA receptor complex that links NMDA receptor activation to AMPA receptor trafficking [28]. Furthermore, a study conducted by Gross et al. [25] showed that PIKE reduction rescued PI3K-dependent and -independent neuronal defects in FXS. Furthermore, it was determined that the genetic reduction of CenG1A, Drosophila's ortholog of PIKE, rescued excessive PI3K signaling and mushroom body defects [29]. In this study, even though CenG1A was hypermethylated, its protein levels were higher in the FXS-like group compared to the control group. This might be another factor contributing to the development of FXS-like symptoms in that group, especially since it was shown in the above study that the increased expression of PIKE, which is encoded by CenG1A, mediates deficits in synaptic plasticity and behavior in Fragile X syndrome [29].
On the other hand, some genes were not associated with ASD nor FXS but had significant roles in brain functions. For example, Nbea is implicated in developmental diseases with early generalized epilepsy phenotypes [49]. Overexpression of Gart improves the IQ of patients with Down syndrome [50]. Hypermethylation of Agap1 is associated with neurodevelopmental disorders through interference with the cellular structure and internal functioning of neurons and through interactions with other genes, particularly the Dtnbp1 gene [52]. Considering the neurological disorders that are connected to hypermethylation of Agap1, one potential impact of L. reuteri that can be inferred from the presented mouse experiment is that it reduced the methylation of Agap1 in mice that received the probiotic in utero, making the mice more neurotypical. In addition, hypermethylation of Agpat1, as well as its loss-of-function mutations, have been found to correlate with metabolic and reproductive disorders, seizures, and abnormal hippocampal brain development [66]. Fbxo11 was shown to be involved in several neurological disorders. Damage to this gene correlates with intellectual disabilities, autism spectrum disorder, ADHD, and anxiety [67].
Thus, these findings may help to explain why these different proteins interact together and might collectively contribute to the observed FXS-like symptoms in the cohort of mice included in this study.
It is believed that administration of L. reuteri might be transferred from mothers to pups, which could halt alteration in the expression of the genes discussed, hence preventing the development of FXS-like symptoms in the offspring. In our studies (see the companion manuscript in this same issue), C-section rederivation with cross-fostering was used to eliminate the possibility of mother-to-infant transfer of microbiota.
Our findings are of importance as the observed Fragile X-like syndrome (FXLS), in mice, is based on the maternal microbiome, rather than the presence of the Fmr1 CGG expansion, and because a link between maternal gut dysbiosis and the development of FXS has never been reported.
Thus, the gut microbiome might contribute to the development of behavioral and physical features associated with FXS, although the mechanisms by which altered intestinal bacteria may affect brain development and function are not currently known.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes13081300/s1, Figure S1: Overlap of CpGs that passed quality control (PassedQC) and with unique best alignment score (UniqueAlign). Only the CpGs with unique best alignment scores and that passed quality control are used in the analysis to find different methylated CpGs. Figure S2: β value distribution of (A) all CpGs that passed quality control and (B) CpGs that passed quality control and with a unique best alignment score. Table S1: List of differentially methylated CpG sites and their corresponding genes. Table S2: List of significant Gene Ontology pathways of three main categories: biological processes (BP), cellular components (CC), and molecular functions (MF).  Data Availability Statement: Data and results generated from this project will be fully available upon request. Biological samples from subjects included in this study will be available under MTA agreement accordingly to the University of California, Davis policy.