Specific pathogen-free (SPF) pregnant mice (n = 3) and germ-free (GF) pregnant mice (n = 3) purchased from GemPharmatech Co., Ltd. were dissected on embryonic day 18 (E18), and the fetal tissues (brain, intestine, and liver) were collected and stored at −80°C for subsequent analyses.

A total of 40 mg mouse fecal pellets were suspended in 200 μl lysis buffer (5 mM EDTA, 0.2% SDS, 0.2M NaCl, and 0.1M Tris-HCl) supplemented with 4 μl of 20 mg/ml proteinase K. The mixtures were disrupted with a grinding rod and then incubated at 56 °C for 6 h. After centrifugation, the supernatant was used for 16S rRNA gene amplification, and the PCR products were visualized on a 2% agarose gel stained with ethidium bromide under UV light. The supernatant from CONV and ABX mice was used for the 16S rRNA gene qPCR analysis. The 16S rRNA gene was detected using two sets of universal bacterial primers: 27F and 1492R; 8F and 1541R. The primers are listed in Supplementary Table S1.

Frozen tissues were homogenized and lysed in RIPA buffer (50 mM Tris-HCl pH 7.5, 1% Nonidet P-40, 0.5% sodium deoxycholate, 0.05% SDS, 1 mM EDTA, and 150 mM NaCl) with freshly added phosphorylase inhibitors and protease inhibitors, and then centrifuged for 20 min at 13,000 × g. The supernatant was aspirated and loaded for the Western immunoblotting analysis. The following antibodies are used: METTL3 (A8370, Abclonal, 1:1,000), METTL14 (HPA038002, Sigma-Aldrich, 1:1,000), FTO (27226-1-AP, Proteintech, 1:1,000), ALKBH5 (16837-1-AP, Proteintech, 1:1,000), and β-actin (66009-1-Ig, Proteintech, 1:5,000).

Fetal mouse tissues were homogenized in 1 ml of TRNzol Universal Reagent (TIANGEN) with glass beads using a LUKYM-I homogenizer, and total RNA was isolated following the manufacturer’s protocol. mRNA was separated from total RNA using a Dynabeads mRNA purification kit (Thermo Fisher Scientific), with two rounds of purification.

LC-MS/MS was performed essentially as described previously (Li et al., 2020). In brief, purified mRNA was digested to nucleosides by nuclease P1 and CIAP, and then it was diluted to 10 ng/μl using nuclease-free water. The samples were filtered and injected into an Agilent Poroshell 120 column coupled online to an AB SCIEX Triple Quad 5500 LC mass spectrometer (Applied Biosystems) in a positive electrospray ionization mode. Concentrations of m⁶A and A were determined based on standard curves of the nucleosides, and the m⁶A/A ratio was calculated.

Total RNA (5 μg) from fetal mouse tissues was reverse-transcribed using a GoScript Reverse Transcription System (Promega), and quantitative real-time PCR was executed using a 2 × RealStar Green Power Mixture (GenStar). The fluorescence intensity of the amplification process was monitored using a LightCycler96 system (Roche). The primers are listed in Supplementary Table S2.

MeRIP experiments were executed as previously reported (Xiao et al., 2019). In brief, approximately 90 µg of total RNA was fragmented into 100- to 300-nucleotide (nt)-long fragments by zinc acetate, followed by the addition of EDTA to terminate the reaction. Then, 5 µg of fragmented RNA was taken as the input control and the remainder was incubated with m⁶A antibodies (4 μg, Abcam, ab151230) in IP buffer (150 mM NaCl, 0.05% NP-40, and 10 mM Tris-HCl) containing RNase inhibitor (Promega), and the mixture was subsequently bound to wash Dynabeads protein G (Invitrogen). After stringent wash, the m⁶A-containing fragments were eluted by competition with 1 mg/ml N6-methyladenosine (Selleck Chemicals). Both the immunoprecipitated RNA fragments and the input RNA were ultimately extracted for library construction using a SMARTer Stranded Total RNA-Seq Kit v2 - Pico Input Mammalian (Takara) following the manufacturer’s protocol. We then performed sequencing using an Illumina Nova platform.

Prior to mapping, all raw data were filtered to remove adapters, and low-quality reads using Trimmomatic (Bolger et al., 2014). Reads of all samples that mapped to rRNA FASTA sequences from UCSC gene annotation (mm10) using bowtie2 (Langmead and Salzberg, 2012) were discarded, and the remaining reads were aligned to the mouse reference genome (GRCm38) using HISAT2 (Pertea et al., 2016). Then mapped files were filtered to keep unique and high mapping quality reads for further analysis using Picard and SAMtools (Li et al., 2009).

m⁶A peaks were identified using MeTPeak. A custom transcriptome annotation file, assembled by StringTie (Pertea et al., 2016) using all sample reads, was created to include intronic and intergenic m⁶A peaks. All other parameters were set to the default settings. The annotatePeaks.pl script from the Homer software suite (Heinz et al., 2010) was used for m⁶A peak annotation.

m⁶A peaks identified in all samples were merged, and featureCounts (Liao et al., 2014) was used to count the fragments that were mapped to the merged peaks. The normalized fragment counts of each peak in MeRIP-seq (MFPKM) were calculated using (methylated fragment counts mapped to the peak × 10⁹)/(length of the peak × total counts of the mapped fragment), and the normalized fragment counts of each peak in input-seq (IFPKM) were calculated using (input fragment counts mapped to the peak × 10⁹)/(length of the peak × total counts of mapped fragments). The methylation level was then calculated for each peak by dividing the MFPKM by the IFPKM. The Pearson correlation coefficient of log2-scaled m⁶A levels across all samples was calculated using corrplot to represent the similarity of each sample.

m⁶A peaks were used for motif search using the findMotifsGenome.pl script from the Homer software suite, using “-rna” and “-len 5” parameters. The R package Guitar (Cui et al., 2016) was used to analyze and plot the distribution of m⁶A on mRNA.

The regions in which the GF group mean m⁶A level was 1.5 fold higher than the SPF group mean m⁶A level were defined as GF group up regions. Also, the regions in which the GF group mean m⁶A level was 1.5 fold lower than the SPF group mean m⁶A level were defined as GF group down regions.

Differentially methylated regions were assigned to mouse genes using the annotatePeaks.pl script from the Homer software suite. The gene list was used for pathways and GO term enrichment using the clusterProfiler (Wu T. et al., 2021).

The input RNA and the immunoprecipitated RNA fragments from mouse fetal tissues were reverse-transcribed using a GoScript Reverse Transcription System (Promega), and then they were analyzed using real-time qPCR. The ratio of immunoprecipitated RNA to the input of each peak was calculated and normalized to GAPDH. The primers are listed in Supplementary Table S3.

Mouse embryonic stem cell line E14TG2a (mES cells) was cultured with the N2B27 base medium supplemented with 1 mM glutamine (Invitrogen), 1% nonessential amino acids (Invitrogen), 0.15 mM 1-thioglycerol (Sigma), 100 U/ml of penicillin–streptomycin (Invitrogen), 25 μg/ml of BSA (Sigma), 1 μM MEK inhibitor PD0325901 (Selleck Chemicals), 3 μM GSK3β inhibitor CHIR99021 (Selleck Chemicals), 2% KOSR (Thermo Fisher), and 1000 U/ml of ESGRO leukemia inhibitory factor LIF (Millipore) on plates coated with 0.2% gelatin.

The Mettl3 ^–/– mES cell line was generated using CRISPR-Cas9 as described previously (Shalem et al., 2014) and the sgRNA sequences are shown in Supplementary Table S4. In brief, sgRNAs were designed on http://crispr-era.stanford.edu/and cloned into the pXPR_001 plasmid. Then, pXPR_001 plasmid was transfected into mES cells using Lipofectamine 3000 (Invitrogen, L3000015). After 12 h, 3 μg/ml of puromycin was added and resistant cells were plated for single colony isolation. Colonies with the desired mutation were identified by Sanger sequencing.

mES cells cultured in 12-well plates at 70–80% confluency were treated with actinomycin D (5 μg/ml final concentration, MCE, HY-17559) for 0, 2, 4, and 8 h before being collected for the extraction of total RNA. RNA was then reverse-transcribed using GoScript Reverse Transcriptase (Promega), and analyzed using real-time qPCR. Expression levels of RNA were calculated and normalized to GAPDH first, and then to the 0 h time point. The mRNA stability of genes was estimated by the half-life of mRNA and calculated using GraphPad Prism 5.0. The primers are listed in Supplementary Table S3.

All of the mice were group-housed in a temperature-controlled (22 ± 1 °C) room with a 12:12-h light:dark cycle, and they had free access to food and water. Mettl3 ^flox/+ mice were generated by Cyagen by inserting loxP sites with the same direction on both sides of exons 2 and 3 of the Mettl3 gene. Male Mettl3 ^flox/+ mice were crossed with female Mettl3 ^flox/+ mice to obtain Mettl3 ^flox/flox mice. Next, Mettl3 ^flox/flox mice were first crossed with DDX4-Cre mice to generate Mettl3 ^flox/+; DDX4-Cre mice, and the latter were then crossed with wild-type mice to generate Mettl3 ^−/+ heterozygous mice. The genotype of each mouse was determined using the genomic DNA extracted from tail tissue.

To mimic GF status, conventional mice (CONV) were treated with antibiotics (ABX), based on methods previously described (Vuong et al., 2020). In brief, 10- to 12 -weeks-old female mice were provided with a mixture of four antibiotics (vancomycin 0.5 g/L, neomycin 1 g/L, ampicillin 1 g/L, and amphotericin-B 0.1 g/L) in their water for 1 week. Female mice were then paired with male mice and gestational day 0.5 was determined by observation of a copulatory plug. Pregnant mice (n = 3) were maintained on ABX in their drinking water until embryonic day 18 (E18), and then dissected to obtain fetal tissues (brain and intestine).

We expressed our measurement data as mean ± SEM. T tests were used for comparisons between two groups. Significant differences were represented by asterisks as follows: *p < 0.05, **p < 0.01, and ***p < 0.001, and ns, not significant.

We initially collected fecal pellets from germ-free (GF, n = 3) and specific pathogen-free (SPF, n = 3) pregnant mice, and the absence of intestinal microbiota in the GF mice was confirmed by 16S rRNA gene amplification (Supplementary Figure S1A). We, then, examined the levels of m⁶A regulators in the mouse fetal brain, intestine, and liver, including writers (METTL3 and METTL14), erasers (FTO and ALKBH5), and readers (YTH-domain family proteins). Using RT-qPCR, we determined that mRNA levels of m⁶A writers and erasers are highly expressed in the fetal brain and intestine from GF pregnant mice (hereafter designated GFB and GFI, respectively) compared to the corresponding tissues from SPF pregnant mice (hereafter designated SPFB and SPFI, respectively). However, the differences in m⁶A reader expression levels are much less marked (Figures 1A,B). Nevertheless, the expression of these proteins is similar in the fetal livers of these two types of mice (Supplementary Figure S1B). A similar tendency in the alteration of protein expression is also uncovered using the Western blotting analysis (Figures 1C,D). Taken together, these results indicated that loss of the maternal microbiome alters the expression of m⁶A writers and erasers in the fetal brain and intestine.

To further investigate whether the maternal microbiome participates in modulating the m⁶A epitranscriptome of offspring, we first detected total m⁶A levels of mouse fetal tissues. We did not observe an apparent change in the global mRNA m⁶A levels between SPF and GF mice as revealed by LC-MS/MS (Supplementary Figure S2A). We, thus, characterized m⁶A methylomes of both mouse fetal brain and intestine (SPFB and GFB and SPFI and GFI—using two independent biological replicates for both) by an m⁶A-immuno-coprecipitation sequencing (MeRIP-seq) analysis. The samples of the same tissue type were clustered well (Figure 2A) and the classic GGAC motif was observed in the fetal brain and intestine (Figure 2B). In agreement with previous studies (Dominissini et al., 2012; Roundtree et al., 2017), the distribution of m⁶A signals around mRNA in the two types of fetal tissue samples was mostly presented in the CDS and 3’UTR, and to a lesser extent in the 5’UTR (Figure 2C). We identified the numbers of m⁶A peaks from these fetal tissues (17,526 in SPFB, 16,885 in GFB, 14,436 in SPFI, and 13,781 in GFI), and we ascertained that approximately three-fourths of the m⁶A peaks overlapped in both fetal brain and intestine (Figure 2D). Compared with SPFB, GFB showed some changes in patterns of m⁶A peaks, with a relative elevation in exonic (SPF 26.81% vs. GFB 28%) and intronic regions (SPF 27.86% vs. GFB 29.38%), and a relative diminution in the 3’untranslated region (3’UTR) from 18 to 16.96% (Figure 2E). Compared with SPFI, the GFI also showed some alterations in patterns of m⁶A peaks with a relative augmentation in exonic regions (SPFI 30.68% vs. GFI 32.5%), and a relative reduction in intronic regions from 30.12 to 29.18% and intergenic regions of 7.7–6.41% (Figure 2F).

To investigate the dynamic characteristics of m⁶A methylation, we further analyzed the differential m⁶A peaks in mouse fetal tissues. As shown in Figure 3A, GFB manifested 2072 upregulated m⁶A peaks and 583 downregulated m⁶A peaks (with the criterion of fold-change ≥1.5). In further examination of the genomic distribution in all three mRNA regions of differential m⁶A peaks, we demonstrated that a majority of the differential m⁶A peaks were in CDS and 3’UTR (Supplementary Figure S3A). Mapping these reads of differential m⁶A peaks to the genome, we identified 1147 genes with upregulated m⁶A peaks and 496 genes with downregulated m⁶A peaks (Supplementary Figure S3B). To further study the biological significance of dysregulated m⁶A modifications in the fetal brain, we conducted GO analyses of differentially m⁶A-methylated genes (Figure 3B and Supplementary Figure S3C). We concentrated on the function of m⁶A-hypermethylated genes and showed that these genes were significantly enriched in pathways related to neurodevelopment, such as synapse formation and axonogenesis. The read coverage plot of a representative gene Cabp1 associated with neurodevelopment was depicted in Figure 3C, and the m⁶A levels of genes (Sema4c, Cobl, Cabp1, Insr, Ntng2, Gabrg2, and Plxna3) were increased in GFB as revealed by using the MeRIP-qPCR analysis (Figure 3D). In addition, the transcript levels of these genes were confirmed by using the RT-qPCR analysis (Supplementary Figure S3D). Collectively, these data suggest that the maternal microbiome regulates the m⁶A of neurodevelopment genes in the mouse fetal brain.

As shown in Figure 4A, GFI reflected 2068 upregulated m⁶A peaks and 184 downregulated m⁶A peaks (with a fold-change ≥1.5). Further examination of the genomic distribution in all three mRNA regions of the differential m⁶A peaks revealed that most of the differential m⁶A peaks were in CDS and 3′UTR (Supplementary Figure S4A). When we mapped these reads of differential m⁶A peaks to the genome, we identified 1590 genes with upregulated m⁶A peaks and 166 genes with downregulated m⁶A peaks (Supplementary Figure S4B). To further assess the biological significance of dysregulated m⁶A modification in the fetal intestine, we executed GO analysis of differentially m⁶A-methylated genes (Figure 4B and Supplementary Figure S4C). When we concentrated on the functions of m⁶A-hypermethylated genes, we found that they were significantly enriched in the Wnt signaling pathway. The read coverage plot of a representative gene Wnt4 is shown in Figure 4C. The differential m⁶A levels of representative genes (Wnt4, Fzd5, Fzd8, Sulf1, Sox13, Axin2, and Abl2) were confirmed by the MeRIP-qPCR analysis (Figure 4D), and their transcript levels were all attenuated in GFI compared to SPFI as revealed by the RT-qPCR analysis (Figure 4E). This indicates that differential m⁶A modifications in these two types of fetal intestines are correlated with the expression of genes enriched in the Wnt signaling pathways. Next, we knocked out Mettl3 (Mettl3 ^-/-) in the mES cell line using CRISPR/Cas9, and we consistently found that Mettl3 knockout significantly decreased m⁶A levels of representative genes while increasing mRNA expression levels (Figures 4F,G). We further investigated whether the changes in m⁶A methylation would affect mRNA levels of representative genes in mESC. We observed that in the presence of actinomycin D (an inhibitor of mRNA transcription), Mettl3 knockout retards the degradation of representative genes mRNAs (Figure 4H). Collectively, these data suggest that the maternal microbiome regulates fetal intestinal m⁶A-modified genes in the Wnt signaling pathway.

To confirm the aforementioned results, we treated CONV pregnant mice with a mixture of four antibiotics (vancomycin, neomycin, ampicillin, and amphotericin-B) to mimic germ-free status (ABX mice) and validated that intestinal microbiota were almost exhausted by the 16S rRNA gene qPCR analysis (Supplementary Figure S5A). Similar to our previous experimental results, the mRNA expression levels of Mettl3 and Fto in the ABX fetal brain were slightly higher than those in the CONV fetal brain (Figure 5A), while the mRNA expression levels of both m⁶A writers and erasers in the ABX fetal intestine were significantly increased compared to the CONV fetal intestine (Figure 5B). In addition, the expression of these proteins remained unchanged in fetal livers from both CONV and ABX (Supplementary Figure S5B). As for the protein expression levels of m⁶A writers and erasers, we noted a universal tendency for them to increase in the ABX fetal brain and intestine (Figures 5C,D). We then determined the m⁶A levels and the expression of representative genes regulated by the maternal microbiome in the ABX and CONV fetal brain and intestine, and we found that the m⁶A levels of these genes in the ABX brain and intestine were also increased relative to CONV (Figures 5E,F), and their transcript levels were confirmed by the RT-qPCR analysis (Figures 5G,H). Collectively, these data show that antibiotic treatment mostly recapitulates m⁶A alterations in the mouse fetal intestine and brain.

To further confirm that the expression of developmental genes was regulated by m⁶A as programed by the maternal microbiome, we generated Mettl3 ^−/+ heterozygous mice (Supplementary Figures S6A,B). Because the homozygous knockout of Mettl3 was embryonically lethal, we crossed Mettl3 heterozygous knockout male mice (Mettl3 ^−/+) with wild-type (WT) female mice. The latter were provided with water (i.e., the offspring of CONV and Mettl3 ^−/+ mice) or ABX (i.e., the offspring of ABX and ABX + Mettl3 ^−/+mice). As expected, there were no significant differences in m⁶A levels of representative genes between CONV and Mettl3 ^−/+ fetal intestines, however, m⁶A levels of representative genes increased in the ABX fetal intestine but not in the ABX + Mettl3 ^−/+ fetal intestine, compared with the CONV fetal intestine (Figure 6A). Correspondingly, mRNA expression levels of representative genes showed no differences between CONV and Mettl3 ^−/+ fetal intestines, while they were significantly reduced in the ABX fetal intestine but not in the ABX + Mettl3 ^−/+ fetal intestine, compared with the CONV fetal intestine (Figure 6B). Our collective results, therefore, indicate that the maternal microbiome affects the developmental gene expression via m⁶A modifications.