A huge number of commensal bacteria inhabit the intestine, which is equipped with the largest immune system in the body. Recently, the regulation of various physiological functions of the host by these bacteria has attracted attention. In this study, the effects of commensal bacteria on gene expression in colonic epithelial cells (CoECs) were investigated with focus on regulation of DNA methylation. RNA sequencing analyses of CoECs from conventional, germ-free, and MyD88−/− mice indicated that, out of the genes affected by commensal bacteria, those downregulated in a MyD88-independent manner were most frequently observed. Furthermore, when the 5′ regions of genes downregulated by commensal bacteria in CoECs were captured using a customized array and immunoprecipitated with the anti-methyl cytosine Ab, a certain population of these genes was found to be highly methylated. Comprehensive analysis of DNA methylation in the 5′ regions of genes in CoECs from conventional and germ-free mice upon pull-down assay with methyl-CpG–binding domain protein 2 directly demonstrated that DNA methylation in these regions was influenced by commensal bacteria. Actually, commensal bacteria were shown to control expression of Aldh1a1, which encodes a retinoic acid–producing enzyme and plays an important role in the maintenance of intestinal homeostasis via DNA methylation in the overlapping 5′ region of Tmem267 and 3110070M22Rik genes in CoECs. Collectively, it can be concluded that regulation of DNA methylation in the 5′ regions of a specific population of genes in CoECs acts as a mechanism by which commensal bacteria have physiological effects on the host.

The significance of gut microbiota in the maintenance of health has recently attracted considerable attention, as increasing evidence demonstrates that the gut microbiota regulates various physiological functions of the host (14). In accordance with this, it has been shown that dysbiosis of the intestinal ecosystem is correlated with a wide variety of diseases, including inflammatory bowel disease, allergy, cancer, autism, and metabolic syndrome (3, 59). Although it is difficult to determine whether such differences in gut microbiota are the cause or the result of these diseases, studies have demonstrated that change in microbiota or specific bacteria is actually involved in the onset, pathogenesis, and prevention of the diseases (1012). Therefore, many efforts have been made to prevent the onset or alleviate the symptoms of the diseases by intervention to the gut microbiota (1317).

As the intestine is equipped with the largest immune system in the body, it is considered that commensal bacteria affect the development and function of the intestinal immune system, particularly in the period immediately after birth when the immune system is in the process of maturation (1820). Moreover, specific mechanisms might be required to establish and maintain symbiosis with commensal bacteria in the intestine considering that commensal bacteria are immunologically recognized as foreign Ags but are not excluded completely. Although induction of immune reactions involves inflammation, inflammatory reactions are strictly controlled at low physiological levels in the intestine despite such a large amount of bacterial Ags. Actually, excessive inflammation is often observed in diseases associated with dysbiosis. Interestingly, intestinal bacteria themselves contribute to this regulation by affecting host cells; however, the precise molecular mechanisms remain to be elucidated.

Based on this context, elucidation of the physiological effects of commensal bacteria on host cells becomes key evidence to clarify the role of commensal bacteria in intestinal homeostasis. For this purpose, the effects of commensal bacteria on gene expression in intestinal epithelial cells (IECs) covering the intestinal mucosa and the underlying regulatory mechanisms were investigated in this study. IECs form a front line of defense by separating the intestinal tract from the internal milieu and are usually exposed to commensal bacteria inhabiting the intestinal tract. Therefore, commensal bacteria have a great impact on gene expression in IECs. In addition, it has been shown that commensal bacteria confer epigenetic effects to specific genes in IECs (2123). DNA methylation, as well as posttranslational histone modifications and noncoding RNA, are important mechanisms mediating epigenetic regulation of gene expression. We previously reported that DNA methylation of the gene encoding TLR4, an innate immune receptor that recognizes LPS of Gram-negative bacteria, in colonic epithelial cells (CoECs) is induced by commensal bacteria (23). As stimulation with LPS through TLR4 induces strong inflammatory reactions, stringent control of TLR4 gene expression is required to prevent excessive inflammation and thereby maintain homeostasis in the intestinal ecosystem; this means that commensal bacteria themselves contribute to establish the symbiotic relationship with the host in the intestine. Although it has been shown that commensal bacteria have potential effects on gene expression and modulate DNA methylation of specific genes, it remains to be clarified for majority of genes how the change in DNA methylation by commensal bacteria is involved in transcriptional regulation. In this study, the effects of commensal bacteria on gene expression and DNA methylation in IECs were analyzed comprehensively to determine their relationship.

Wild-type (WT) and MyD88−/− BALB/c mice were purchased from CLEA Japan (Tokyo, Japan) and Oriental BioService (Kyoto, Japan), respectively. Mice were bred under conventional (CV) or germ-free (GF) conditions as described previously (24). Female mice were used at 9–13 wk of age. All experiments were approved by the Nihon University Animal Care and Use Committee and conducted in accordance with their guidelines.

Small IECs (SIECs) and CoECs were prepared from the whole small intestine and the whole colon of mice, respectively. Epithelial cell preparation was performed as described previously (25). After removing Peyer patches, the tissues were cut into 2–3-mm pieces and washed in HBSS supplemented with 1 mM DTT and 0.5 mM EDTA accompanied by shaking. The tissues were then treated with dispase (BD Biosciences, Franklin Lakes, NJ) to collect single-cell suspensions. Lymphocytes were depleted by MACS using Dynabeads M-450 Streptavidin (Invitrogen, Thermo Fisher Scientific) and biotin-conjugated anti-CD45 Ab (eBioscience, San Diego, CA).

The mouse IEC line CMT-93 developed from rectal carcinoma was purchased from DS Pharma Biomedical (Osaka, Japan) and cultured in DMEM supplemented with 10% (v/v) FBS (Biowest, Nuaillé, France), 100 U/ml penicillin, 100 μg/ml streptomycin, and 5 × 10−5 M 2-ME at 37°C in a humidified incubator with 5% CO2.

Total RNA was prepared from SIECs and CoECs of WT-CV, MyD88−/−-CV, and WT-GF mice by pooling from six to eight mice per group using the High Pure RNA Isolation Kit (Roche, Basel, Switzerland). RNA sequencing (RNA-seq) was performed by INFOBIO (Tokyo, Japan). Briefly, mRNA was extracted, reverse transcribed, treated with restriction enzyme NlaIII, and ligated with adaptor sequences. The library was prepared by processing the tag with EcoP15I and analyzed using Genome Analyzer IIx (Illumina, San Diego, CA). DNA sequences originated from the library were extracted from the resulting reads, and CATG sequence was added to the 5′ ends. The acquired tags (26 bp from the 5′ end) were then mapped to Mus musculus mRNA reference sequences (National Center for Biotechnology Information RefSeq; ftp://ftp.ncbi.nlm.nih.gov/refseq/M_musculus/mRNA_Prot/) using Maq 0.7.1 (Wellcome Trust Sanger Institute). To compare mRNA expression among CV, MyD88, and GF mice, tag numbers obtained for each gene were corrected to the total numbers of tags obtained (22,577,586, 21,832,965, and 24,226,181 for CV, GF, and MyD88−/−, respectively). If these corrected values differed by more than 3-fold between CV and GF mice, except for genes with unquantifiable or low expression to which only one digit tags were mapped in both CV and GF mice, the expression of the corresponding gene was considered to be affected by commensal bacteria. MyD88 dependency for regulation of these genes by commensal bacteria was determined if the values from results of MyD88−/− mice were close to those of CV or GF mice. Next-generation sequencing (NGS) data have been deposited in the DNA Data Bank of Japan (DDBJ) Sequence Read Archive (https://www.ddbj.nig.ac.jp/dra) under accession number DRA008906.

For the quantitative analysis of mRNA expression, total RNA was prepared from cells using the High Pure RNA Isolation Kit (Roche). Total RNA was then reverse-transcribed using SuperScript IV Reverse Transcriptase (Invitrogen, Thermo Fisher Scientific) and oligo(dT)20 primers (Invitrogen, Thermo Fisher Scientific). SYBR Green I Master reagent (Roche) was used to quantify cDNA by real-time PCR using LightCycler 480 (Roche) with the cycling condition of 95°C for 10 s, 58°C for 10 s, and 72°C for 12 s. The relative expression of each gene was calculated from the calibration curve created using a dilution series of standard samples followed by normalization to that of Gapdh. Specific amplification of each gene was confirmed by melting curve analysis. Information on primers is provided in Supplemental Table I.

Genomic DNA was prepared from CoECs of CV or GF mice by pooling from two independent experiments (six to seven mice per experiment) and fragmented by nebulization (N2 0.23 MPa, 6 min) in shearing buffer (Tris–EDTA [TE; pH 8], 10% glycerol). Sheared DNA was end repaired, followed by dA-tailing and adaptor ligation using the NEBNext DNA Sample Prep Master Mix Set (New England Biolabs, Ipswich, MA) according to the manufacturer’s instructions. DNA fragments obtained were electrophoresed on agarose gel, cut out from the gel for purification using a Gel Extraction Kit (QIAGEN, Hilden, Germany). Then, methylated DNA fragments were precipitated with methyl-CpG–binding domain protein 2 (MBD2)–coupled beads using the MethylMiner Methylated DNA Enrichment Kit (Applied Biosystems, Thermo Fisher Scientific), as described in the protocol supplied by the manufacturer. DNA fragments collected were amplified by PCR using adapter sequences as primers according to the user’s guide of the NimbleGen SeqCap EZ Library SR (Roche). The library thus prepared was analyzed using HiSeq 2000 (Illumina) to acquire 101-bp single-end reads at INFOBIO. Identified reads were mapped to the mouse genome (mm10) using the Burrows–Wheeler Aligner (version 0.6.2-r126) and then annotated to the regions near the transcription start site of each gene (from nt −1000 to +501 on both the same and the opposite strands; nucleotide numbers are counted from the transcription start site as +1). Extracted reads were compared quantitatively between CV and GF mice. NGS data have been deposited in the DDBJ Sequence Read Archive (https://www.ddbj.nig.ac.jp/dra) under accession number DRA008905.

Genomic DNA was prepared from cells using a PureLink Genomic DNA Mini Kit (Invitrogen, Thermo Fisher Scientific), and methylation status of CpG motifs was analyzed as follows. Genomic DNA (1 μg) was denatured with 6 N NaOH and modified by the bisulfite conversion reaction using a BisulFast DNA modification kit (TOYOBO, Osaka, Japan). The 388-bp 5′ region (nt −102/+286; nucleotide numbers are counted from the transcription start site as +1.) of the Tmem267 gene was amplified by PCR from the modified DNA. Sequences of the synthetic oligonucleotides used as PCR primers are as follows: 5′-GGTGGAGAGTTGAGGTTTTTTGTG-3′ (forward) and 5′-CCCCTAACCCAAACTACCTTCATC-3′ (reverse). PCR products were cloned into the pCR2.1 vector, and nucleotide sequences of 16–17 independent clones were analyzed.

To construct plasmids for RNA interference (RNAi) against Tmem267 and 3110070M22Rik RNAi, synthetic oligo DNA listed below was respectively annealed by denaturing at 95°C for 4 min and subsequently cooling: Tmem267 forward, 5′-TGCTGTACATCAGGAATGAGCACAGGGTTTTGGCCACTGACTGACCCTGTGCTTTCCTGATGTA-3′; Tmem267 reverse, 5′-CCTGTACATCAGGAAAGCACAGGGTCAGTCAGTGGCCAAAACCCTGTGCTCATTCCTGATGTAC-3′; 3110070M22Rik RNAi forward, 5′-TGCTGTCATGATCCAGCCTTGAACTTGTTTTGGCCACTGACTGACAAGTTCAACTGGATCATGA-3′; and 3110070M22Rik reverse, 5′-CCTGTCATGATCCAGTTGAACTTGTCAGTCAGTGGCCAAAACAAGTTCAAGGCTGGATCATGAC-3′

Double-stranded oligo DNA thus obtained was introduced into the pCDNA6.2-GW/EmGFP-miR vector, respectively, using the BLOCK-iT PolII miR RNAi Expression Vector Kit (Invitrogen, Thermo Fisher Scientific). Resulting clones or the negative control plasmid provided with the kit were introduced into CMT-93 cells using the X-tremeGene HP Reagent (Roche). After selecting transfected cells with 4 μg/ml blasticidin for 10–14 d, total RNA was prepared for quantitative RT-PCR (qRT-PCR) analysis.

Cells were washed with ice-cold PBS and incubated on ice for 30 min in the lysis buffer (20 mM Tris [pH 7.5], 150 mM NaCl, 1 mM EDTA, 60 mM n-octyl-β-d-glucoside, and 1% Nonidet P-40) supplemented with protease inhibitor mixture (Nacalai Tesque, Kyoto, Japan). After centrifugation at 4°C, 20,000 × g, for 10 min, the supernatants were collected. The protein content of the lysates was measured using the BCA Protein Assay Kit (Pierce, Thermo Fisher Scientific), and equal amounts of protein were analyzed by immunoblotting using anti-ALDH1A1 (ab52492; Abcam, Cambridge, U.K.) and anti–β-Actin (ab49900; Abcam) Abs.

Small intestine and colon were surgically excised from mice. Small intestines were divided into five pieces of equal length, and the first (proximal), the third (medial), and the fifth (distal) parts were used for further experiments. Luminal content was gently washed out with Ca2+- and Mg2+-free Hanks’ balanced salt solution (Sigma-Aldrich, St. Louis, MO) containing 5% FBS and was replaced with SCEM-L1 embedding medium (Section-lab, Hiroshima, Japan). Tissues were then embedded in SCEM-L1 and cut into 16-μm sections using the HM550 cryostat (Thermo Fisher Scientific). Tissue sections were fixed with ethanol and stained with toluidine blue. Laser microdissection (LMD) was performed using LMD7000 (Leica Microsystems) to capture segments including IECs at upper (apical) and bottom (basolateral) portions of villi (as shown in Supplemental Fig. 1). Total RNA was prepared from collected tissue segments of approximately equal total area using RNeasy Micro Kit (QIAGEN) for quantification of Aldh1a1 expression by qRT-PCR.

A customized array with DNA sequences covering the nt −1000/+501 regions of genes showed <1/3 expression in CoECs of CV mice than that seen in those of GF mice by RNA-seq was designed and manufactured by Roche (relevant genes are listed in Supplemental Table II). Genomic DNA was prepared from CoECs of CV mice by pooling from two independent experiments (four mice per experiment). Shearing, end repairing, dA tailing, adaptor ligation, and purification were carried out as described above. DNA fragments were hybridized onto the customized array. Hybridization, washing, and collection of hybridized DNA fragments were performed using the SeqCap EZ Reagent (Roche) according to the manufacturer’s instructions. Collected DNA fragments were denatured at 95°C for 5 min and then mixed with the anti–5-methylcytosine Ab (clone 33D3; Abcam) at 4°C overnight with rotation in immunoprecipitation (IP) buffer (10 mM Tris-HCl [pH 8], 150 mM NaCl, 1 mM EDTA, and 0.05% Triton X-100) in the presence of random 20-mer DNA. Ab-bound DNA was then mixed with salmon sperm DNA–treated protein G Dynabeads (Invitrogen, Thermo Fisher Scientific) at 4°C for 2 h with rotation. After five washes with the IP buffer and a subsequent wash with TE, immunoprecipitated DNA was eluted in TE (pH 8) supplemented with 0.25% SDS and 500 μg/ml proteinase K at 55°C for 2 h and was purified using the MinElute Reaction Cleanup Kit (QIAGEN). The eluents and input control before IP were amplified by PCR using adapter sequences as primers and were purified with the MinElute Reaction Cleanup Kit in the presence of random 20-mer DNA as described in the NimbleGen SeqCap EZ Library SR User’s Guide. Libraries obtained were subjected to NGS to acquire 50-bp, single-end reads using HiSeq 2000 (Illumina) at INFOBIO. Identified reads were mapped to the mouse genome (mm9) using Maq (version 0.7.1) and then annotated to nt −1000/+501 regions of genes whose expression was <1/3 in CoECs of CV mice than in those of GF mice. NGS data have been deposited in the DDBJ Sequence Read Archive (https://www.ddbj.nig.ac.jp/dra) under accession number DRA008907.

Differences between two or more groups were analyzed by two-tailed Student t test or one-way ANOVA followed by Dunnett test or Tukey test, respectively. A p value <0.05 was considered statistically significant.

As IECs receive the stimulation from commensal bacteria at the front line and more than 99% of the intestinal bacteria inhabit the colon, the effect of commensal bacteria on gene expression in CoECs was first examined comprehensively. For this purpose, mRNA expression in CoECs was compared among WT mice bred under CV and GF conditions and MyD88−/− mice bred under CV condition (described as CV, GF, and MyD88−/−, respectively) by RNA-seq. The numbers of genes showing more than three-fold difference in expression between CV and GF mice are summarized in Table I. MyD88 dependency for up- and downregulation by commensal bacteria is also specified in Table I based on whether gene expression in MyD88−/− mice is similar to that in CV or GF mice. We found that the number of genes downregulated by commensals was larger than those upregulated by commensals (644 and 388, respectively). Interestingly, of the upregulated genes, MyD88-dependent regulation was predominant (261/388), whereas, of the downregulated genes, MyD88-independent regulation was predominant (517/644).

Table I.
Up- and downregulated genes in CoECs by commensals
Upregulation (CV > GF)Downregulation (CV < GF)
MyD88 independent 119 517 
Intermediate 11 
MyD88 dependent 261 116 
Total 388 644 
Upregulation (CV > GF)Downregulation (CV < GF)
MyD88 independent 119 517 
Intermediate 11 
MyD88 dependent 261 116 
Total 388 644 

Gene expression was systematically analyzed by RNA-seq using RNA prepared from CoECs of CV, GF, and MyD88−/− mice. Numbers of genes expressed more than 3.0-fold (upregulation by commensal bacteria) or <1/3 (downregulation by commensal bacteria) in CV mice than in GF mice are shown. MyD88 dependency was considered based on whether the gene expression in MyD88−/− mice is closer to that in CV mice (MyD88 independent) or GF mice (MyD88 dependent) or just middle of these (intermediate).

Tables II and III respectively show representative genes that were upregulated and downregulated by commensals as seen in the RNA-seq results from CoECs. For example, expression of genes encoding antimicrobial peptides such as regenerating islet-derived 3 (Reg3) and α-defensin were markedly affected by commensal bacteria. Notably, expression of Reg3β and Reg3γ was induced by commensals in a MyD88-dependent manner. In contrast, several α-defensins, including 5, 20, 1, and 24, were downregulated by commensals; the former two were suppressed in a MyD88-independent manner, and in contrast, the latter two were suppressed in a MyD88-dependent manner.

Table II.
Genes induced by commensals in CoECs
GeneCVGFMyD88−/−
MyD88-independent 
 Delta-like non-canonical Notch ligand 1 (Dlk1126 83 
 ATP/GTP binding protein-like 2 (Agbl2107 78 
 Angiogenin, RNase A family, member 4 (Ang42773 269 43,569 
 Gasdermin C4 (Gsdmc43240 408 4102 
 Phospholipase A2, group IIA (Pla2g2a285 49 544 
 Phospholipase A2, group III (Pla2g3148 34 116 
 TGF, β 1 (Tgfβ168 14 56 
 NO synthase 2, inducible (Nos246 12 89 
 Solute carrier family 30, member 10 (Slc30a101080 287 940 
MyD88-dependent 
 Regenerating islet-derived 3 β (Reg3β25,426 2111 1908 
 Regenerating islet-derived 3 γ (Reg3γ71 
 Chemokine (C-X-C motif) ligand 1 (Cxcl1908 252 239 
 TNF (Tnf707 145 164 
 IL-17C (Il17c39 12 
 Fucosyltransferase 2 (Fut21199 340 569 
 Cytochrome P450, family 3, subfamily a, polypeptide 44 (Cyp3a44290 11 64 
 Cytochrome P450, family 2, subfamily c, polypeptide 55 (Cyp2c5533,206 4752 10,046 
 Aquaporin 4 (Aqp45745 862 1128 
 Carbonyl reductase 2 (Cbr2533 87 232 
 Chloride intracellular channel 6 (Clic698 13 30 
GeneCVGFMyD88−/−
MyD88-independent 
 Delta-like non-canonical Notch ligand 1 (Dlk1126 83 
 ATP/GTP binding protein-like 2 (Agbl2107 78 
 Angiogenin, RNase A family, member 4 (Ang42773 269 43,569 
 Gasdermin C4 (Gsdmc43240 408 4102 
 Phospholipase A2, group IIA (Pla2g2a285 49 544 
 Phospholipase A2, group III (Pla2g3148 34 116 
 TGF, β 1 (Tgfβ168 14 56 
 NO synthase 2, inducible (Nos246 12 89 
 Solute carrier family 30, member 10 (Slc30a101080 287 940 
MyD88-dependent 
 Regenerating islet-derived 3 β (Reg3β25,426 2111 1908 
 Regenerating islet-derived 3 γ (Reg3γ71 
 Chemokine (C-X-C motif) ligand 1 (Cxcl1908 252 239 
 TNF (Tnf707 145 164 
 IL-17C (Il17c39 12 
 Fucosyltransferase 2 (Fut21199 340 569 
 Cytochrome P450, family 3, subfamily a, polypeptide 44 (Cyp3a44290 11 64 
 Cytochrome P450, family 2, subfamily c, polypeptide 55 (Cyp2c5533,206 4752 10,046 
 Aquaporin 4 (Aqp45745 862 1128 
 Carbonyl reductase 2 (Cbr2533 87 232 
 Chloride intracellular channel 6 (Clic698 13 30 

Representatives of genes induced by commensals in CoECs in MyD88-independent (upper) and -dependent (lower) manner. Numbers of reads mapped to each corresponding gene in CV, GF, and MyD88−/− are shown accordingly. Values are corrected by the total reads obtained for each sample.

Table III.
Genes suppressed by commensals in CoECs
GeneCVGFMyD88−/−
MyD88-independent 
 Defensin α5 (Defα5143 
 Defensin α20 (Defα2047 
 C-C chemokine ligand 25 (Ccl25125 934 150 
 C-C chemokine ligand 11 (Ccl1145 20 
 Amylase 2a (Amy2a2, 3, 493 1025 
 Sucrase-isomaltase (Sis) (α-glucosidase) 25 868 36 
 Solute carrier family 2 member 5 (Slc2a5) (fructose transporter) 47 158 15 
 Apolipoprotein A-I (Apoa1146 
 Fatty acid binding protein 4 (Fabp417 129 
 Carboxypeptidase A1 (Cpa192 
 Cytochrome P450 (Cyp4b1724 2322 670 
 Adenylate cyclase 9 (Adcy934 132 65 
 IL-22R (Il22ra110 38 19 
 Calcium/calmodulin-dependent protein kinase II a (Camk2a12 
MyD88-dependent 
 Defensin α1 (Defα11321 873 
 Defensin α24 (Defα241341 873 
 Adenosine A1R (Adora139 80 
 Inositol 1,4,5-trisphosphate 3-kinase A (Itpka13 47 31 
 Folate receptor 1 (Folr112 
 H+-K+-ATPase (Atp12a793 2479 2373 
GeneCVGFMyD88−/−
MyD88-independent 
 Defensin α5 (Defα5143 
 Defensin α20 (Defα2047 
 C-C chemokine ligand 25 (Ccl25125 934 150 
 C-C chemokine ligand 11 (Ccl1145 20 
 Amylase 2a (Amy2a2, 3, 493 1025 
 Sucrase-isomaltase (Sis) (α-glucosidase) 25 868 36 
 Solute carrier family 2 member 5 (Slc2a5) (fructose transporter) 47 158 15 
 Apolipoprotein A-I (Apoa1146 
 Fatty acid binding protein 4 (Fabp417 129 
 Carboxypeptidase A1 (Cpa192 
 Cytochrome P450 (Cyp4b1724 2322 670 
 Adenylate cyclase 9 (Adcy934 132 65 
 IL-22R (Il22ra110 38 19 
 Calcium/calmodulin-dependent protein kinase II a (Camk2a12 
MyD88-dependent 
 Defensin α1 (Defα11321 873 
 Defensin α24 (Defα241341 873 
 Adenosine A1R (Adora139 80 
 Inositol 1,4,5-trisphosphate 3-kinase A (Itpka13 47 31 
 Folate receptor 1 (Folr112 
 H+-K+-ATPase (Atp12a793 2479 2373 

Representatives of genes suppressed by commensals in CoECs in MyD88-independent (upper) and -dependent (lower) manner. Numbers of reads mapped to each corresponding gene in CV, GF, and MyD88−/− are shown accordingly. Values are corrected by the total reads obtained for each sample.

Next, the expression of α-Defensin 1, α-Defensin 5, and Reg3β in SIECs and CoECs of CV, GF, and MyD88−/− mice were quantified by qRT-PCR (Fig. 1A). Results from CoECs were similar to those obtained by RNA-seq. However, expression patterns were different between SIECs and CoECs; α-Defensin 1 and 5 showed lower expression in CoECs of CV mice than in those of GF mice, but they were expressed at higher levels in SIECs of CV mice than in those of GF mice. In addition, the suppression of α-Defensin 1 by commensal bacteria in CoECs was MyD88 dependent, whereas that of α-Defensin 5 was MyD88 independent. As α-defensin is mainly produced in the epithelium of the small intestine, which contains Paneth cells, the expression levels of α-Defensin 1 and 5 were much higher in SIECs than in CoECs, whereas the expression of Reg3β was almost comparable in these cells. Expression of genes encoding some chemokines such as CCL25, CCL11, and CXCL1 and cytokines such as TGF-β1, TNF, and IL-17C were also affected by commensals (Tables II, III). When the expression of Ccl25 and Cxcl1 in SIECs and CoECs was quantified by qRT-PCR (Fig. 1B), the effect of commensal bacteria on the expression of these chemokine genes was not found to be as significant in SIECs as compared with that seen in CoECs.

FIGURE 1.

Gut microbiota differently regulates antimicrobial peptide and chemokine expression in SIECs and CoECs.

Total RNA was prepared from SIECs and CoECs of CV, GF, and MyD88−/− mice to quantify mRNA expression of antimicrobial peptides (A) and chemokines (B) by qRT-PCR. Values are normalized using GAPDH mRNA levels and expressed as percentages relative to their expression in SIECs of CV mice. Results are expressed as mean ± SD of four to five independent experiments. Five or six mice were used for each experiment. *p < 0.05, **p < 0.01, ***p < 0.005, ****p < 0.001, *****p < 0.0005 by Tukey test.

FIGURE 1.

Gut microbiota differently regulates antimicrobial peptide and chemokine expression in SIECs and CoECs.

Total RNA was prepared from SIECs and CoECs of CV, GF, and MyD88−/− mice to quantify mRNA expression of antimicrobial peptides (A) and chemokines (B) by qRT-PCR. Values are normalized using GAPDH mRNA levels and expressed as percentages relative to their expression in SIECs of CV mice. Results are expressed as mean ± SD of four to five independent experiments. Five or six mice were used for each experiment. *p < 0.05, **p < 0.01, ***p < 0.005, ****p < 0.001, *****p < 0.0005 by Tukey test.

Close modal

We previously reported that commensal bacteria induce DNA methylation of the gene encoding TLR4, which acts as a sensor for Gram-negative bacteria, in CoECs (23). Thus, we focused on DNA methylation as a mechanism underlying the regulation of gene expression by commensal bacteria. Genomic DNA fragments obtained from CoECs of CV and GF mice were pulled down using MBD2-coupled beads and analyzed by NGS to compare CpG methylation in the 5′ regions of genes between these mice. The reads mapped to the region from nt −1000 to nt +501 were extracted and compared between CV and GF mice for each gene (Fig. 2). Out of 30,395 transcription start sites, 23 and 453 sites showed more than 2.0-fold and 1.5-fold change in the annotation frequency between CoECs from CV and GF mice, respectively (Table IV). These included both types of genes methylated at higher and lower levels in CoECs of CV mice than in those of GF mice. The results directly demonstrated that commensal bacteria affect DNA methylation in the 5′ regions of specific genes.

FIGURE 2.

Gut microbiota alters DNA methylation frequencies in the 5′ regions of genes.

To compare the DNA methylation patterns, genomic DNA was prepared from CoECs of CV and GF mice by pooling from two independent experiments (six to seven mice per experiment). Genomic DNA fragments were pulled down using MBD2-coupled beads and analyzed by NGS. Identified reads were annotated to 1.5-kb regions near the transcription start site (nt −1000/+501) of each gene. The ratio of annotated read number of GF mice relative to that of CV mice is shown for each gene.

FIGURE 2.

Gut microbiota alters DNA methylation frequencies in the 5′ regions of genes.

To compare the DNA methylation patterns, genomic DNA was prepared from CoECs of CV and GF mice by pooling from two independent experiments (six to seven mice per experiment). Genomic DNA fragments were pulled down using MBD2-coupled beads and analyzed by NGS. Identified reads were annotated to 1.5-kb regions near the transcription start site (nt −1000/+501) of each gene. The ratio of annotated read number of GF mice relative to that of CV mice is shown for each gene.

Close modal
Table IV.
Effect of commensals on DNA methylation in the 5′ regions of genes in CoECs
No. of Genes
>1.5-Fold change 453 
>2.0-Fold change 23 
No. of Genes
>1.5-Fold change 453 
>2.0-Fold change 23 

Genomic DNA fragments of CoECs from CV and GF mice were pulled down using MBD2-coupled beads and analyzed by NGS. Values represent numbers of genes with indicated quantitative differences in the reads mapped to their 5′ region between CV and GF mice.

The 5′ regions of Tmem267 and 3110070M22Rik encoding a transmembrane protein 267 and a noncoding RNA, respectively, showed higher methylation in CV mice than in GF mice (Table V). These genes are closely located on different DNA strands of chromosome 13 in mice as shown Fig. 3A. Methylation frequencies of CpG motifs present in the overlapping 5′ region of these genes in CoECs from CV and GF mice were further analyzed by bisulfite sequencing. Similar to the results of the MBD2 pull-down assay, methylation frequency was significantly higher in CV mice than in GF mice (Fig. 3B, 3C). In accordance with these results, gene expression of both Tmem267 and 3110070M22Rik was found to be lower in CV mice than in GF mice (Fig. 3D).

Table V.
Comparison of DNA methylation in the 5′ region of Tmem267 and 3110070M22Rik genes in CoECs between CV and GF mice
CV Strand
GF Strand
SameOppositeTotalSameOppositeTotal
Tmem267 301 240 541 97 88 185 
3110070M22Rik 240 298 538 86 96 182 
CV Strand
GF Strand
SameOppositeTotalSameOppositeTotal
Tmem267 301 240 541 97 88 185 
3110070M22Rik 240 298 538 86 96 182 

Results for the MBD2 pull-down analysis showing DNA methylation in the 5′ region (nt −1000/+501) of Tmem267 and 3110070M22Rik in CoECs from CV and GF mice are summarized. The numbers of reads mapped respectively to the same and the opposite strand as the gene, and those of the total are shown.

FIGURE 3.

Gut microbiota suppresses expression of Tmem267 and 3110070M22Rik genes through induction of DNA methylation in their overlapping 5′ region.

(A) Positional relationship of Tmem267 and 3110070M22Rik genes on chromosome 13. (B and C) Genomic DNA was prepared from CoECs of CV and GF mice by pooling from two independent experiments (five to six mice per experiment), respectively. Methylation of 27 CpG motifs in the 5′ region of Tmem267 and 3110070M22Rik genes was analyzed by the bisulfite reaction. Analysis was performed for 16–17 independent clones. In (B), filled and open circles present methylated and unmethylated motifs, respectively, and clones are arranged in order of their methylation frequency. Data are summarized in (C) as mean ± SD of methylated CpG motifs for the total 27 motifs in the region. (D) Total RNA was prepared from CoECs of CV and GF mice to analyze expression of Tmem267 and 3110070M22Rik by qRT-PCR. Results are expressed as mean ± SD of seven independent experiments. Each experiment was conducted using five to six mice per group. *p < 0.05, **p < 0. 005 by two-tailed Student t test.

FIGURE 3.

Gut microbiota suppresses expression of Tmem267 and 3110070M22Rik genes through induction of DNA methylation in their overlapping 5′ region.

(A) Positional relationship of Tmem267 and 3110070M22Rik genes on chromosome 13. (B and C) Genomic DNA was prepared from CoECs of CV and GF mice by pooling from two independent experiments (five to six mice per experiment), respectively. Methylation of 27 CpG motifs in the 5′ region of Tmem267 and 3110070M22Rik genes was analyzed by the bisulfite reaction. Analysis was performed for 16–17 independent clones. In (B), filled and open circles present methylated and unmethylated motifs, respectively, and clones are arranged in order of their methylation frequency. Data are summarized in (C) as mean ± SD of methylated CpG motifs for the total 27 motifs in the region. (D) Total RNA was prepared from CoECs of CV and GF mice to analyze expression of Tmem267 and 3110070M22Rik by qRT-PCR. Results are expressed as mean ± SD of seven independent experiments. Each experiment was conducted using five to six mice per group. *p < 0.05, **p < 0. 005 by two-tailed Student t test.

Close modal

As the function of Tmem267 is unknown, the effect of Tmem267 knockdown by RNAi in a CoEC line CMT-93 was analyzed. As shown in Fig. 4A, it was confirmed that expression of Tmem267 and 3110070M22Rik was significantly suppressed by each corresponding RNAi. In addition, unexpectedly, RNAi against 3110070M22Rik also caused decreased expression of Tmem267. The effect of RNAi against Tmem267 or 3110070M22Rik on expression of various genes in IECs was then analyzed. Suppression of Tmem267 markedly decreased the expression of Aldh1a1 encoding RALDH1, an enzyme mediating the conversion of retinal, derived from vitamin A, to retinoic acid, whereas it did not affect tight junction-related proteins (Occuldin and Zo-1), pattern recognition receptors (Tlr2 and Tlr4), and negative regulator of TLR signaling (Tollip) (Fig. 4B, Supplemental Fig. 2). In addition, suppression of 3110070M22Rik similarly caused a decrease in Aldh1a1 expression. RALDH1 protein expression was also significantly decreased by knockdown of Tmem267 and 3110070M22Rik (Fig. 4C).

FIGURE 4.

Tmem267 induces expression of Aldh1a1, which produces retinoic acid, an immunomodulator in the intestine, in CoECs.

The mouse IEC line CMT-93 was transfected with Tmem267 or 3110070MRik small interfering RNA (siRNA) expression vector. After selection of transfected cells with blasticidin for 10–14 d, cells were collected for qRT-PCR and Western blotting analyses. (A) To confirm the effect of RNAi, expression of Tmem267 or 3110070MRik in transfected cells was quantified. Results are expressed as mean ± SD of three independent experiments. (B) The effect of RNAi for Tmem267 and 3110070MRik on Aldh1a1 mRNA expression was analyzed by qRT-PCR. Results are expressed as mean ± SD of three independent experiments. (C) The effect of RNAi for Tmem267 and 3110070MRik on RALDH1 protein expression was analyzed by Western blotting. Representative blots of an experiment performed in duplicate (left) and the mean ± SD of the relative band intensities normalized to β-actin from three independent experiments (right) are shown. (D) Tissue sections of apical (a) and basolateral (b) portions of villus epithelia of proximal (Pro), medial (Med), and distal (Dis) small intestine and colon (Co) were collected from CV (n = 4, open bars) and GF (n = 4, filled bars) mice by LMD and pooled so that the total area was almost equal. Total RNA was extracted from the pooled tissue sections of each intestinal portion to analyze Aldh1a1 mRNA expression by qRT-PCR. *p < 0.05, **p < 0.01 by Dunnett test.

FIGURE 4.

Tmem267 induces expression of Aldh1a1, which produces retinoic acid, an immunomodulator in the intestine, in CoECs.

The mouse IEC line CMT-93 was transfected with Tmem267 or 3110070MRik small interfering RNA (siRNA) expression vector. After selection of transfected cells with blasticidin for 10–14 d, cells were collected for qRT-PCR and Western blotting analyses. (A) To confirm the effect of RNAi, expression of Tmem267 or 3110070MRik in transfected cells was quantified. Results are expressed as mean ± SD of three independent experiments. (B) The effect of RNAi for Tmem267 and 3110070MRik on Aldh1a1 mRNA expression was analyzed by qRT-PCR. Results are expressed as mean ± SD of three independent experiments. (C) The effect of RNAi for Tmem267 and 3110070MRik on RALDH1 protein expression was analyzed by Western blotting. Representative blots of an experiment performed in duplicate (left) and the mean ± SD of the relative band intensities normalized to β-actin from three independent experiments (right) are shown. (D) Tissue sections of apical (a) and basolateral (b) portions of villus epithelia of proximal (Pro), medial (Med), and distal (Dis) small intestine and colon (Co) were collected from CV (n = 4, open bars) and GF (n = 4, filled bars) mice by LMD and pooled so that the total area was almost equal. Total RNA was extracted from the pooled tissue sections of each intestinal portion to analyze Aldh1a1 mRNA expression by qRT-PCR. *p < 0.05, **p < 0.01 by Dunnett test.

Close modal

The expression of Aldh1a1 in the epithelium from various intestinal regions including proximal, medial, and distal small intestine and colon of CV and GF mice was further analyzed by LMD, followed by RNA extraction and qRT-PCR (Fig. 4D). Aldh1a1 expression was lower in the epithelial tissues of CV mice than in those of GF mice in all intestinal portions. These results indicated that commensal bacteria suppress Aldh1a1 expression in the intestinal epithelium in vivo, probably through the suppression of Tmem267 and 3110070M22Rik expression via DNA methylation of 5′ regions of these genes. Furthermore, Aldh1a1 expression tended to be lower in the distal small intestine and colon, which are inhabited by a large number of commensal bacteria, than in the proximal and medial small intestine. However, because this tendency was observed in both CV and GF mice, additional factors other than the commensal bacteria, such as food components in the upper small intestine, could also influence the expression of Aldh1a1. In addition, Aldh1a1 expression was tended to be higher in apical portions than in basolateral portions of the villus epithelium, particularly in the proximal small intestine.

As specific gene expression seemed to be regulated by DNA methylation in the 5′ region in CoECs, we examined the methylation status in the 5′ regions of the 644 genes downregulated to <1/3 by commensals, as shown in Table I. For this analysis, genomic DNA fragments from CoECs of CV mice were hybridized using an array with sequences in the 5′ regions (nt −1000/+501) of 562 out of 644 genes whose transcription start sites could be identified. The captured DNA fragments were then immunoprecipitated with the anti–5-methylcytosine Ab and analyzed by NGS. Figure 5 shows frequencies of the reads mapped to the nt −1000/+501 region of each target gene relative to the total number of reads before (input) and after IP. Reads that were mapped to the nt −1000/+501 region were detected for 510 out of 562 genes after IP, indicating that these genes were possibly methylated. Moreover, reads were mapped with high frequency to some specific genes in IP sample, whereas remarkable difference in their frequency was not observed among the genes in input sample, showing that specific genes among those downregulated by commensals in CoECs were highly methylated. In contrast, as reads that were mapped to 35 genes were obtained from the input sample, but not from the IP sample, these genes were judged to be unmethylated. Reads mapped to the remaining 17 genes were not detected in input sample.

FIGURE 5.

A part of gene population downregulated by commensals is methylated at relatively high levels in their 5′ regions.

Genomic DNA was prepared from CoECs of CV mice by pooling from two independent experiments (four mice per experiment). The genomic DNA fragments corresponding to the sequences in nt −1000/+501 regions of 562 genes that were downregulated by commensals were captured using a customized array, immunoprecipitated with the anti-methylcytosine Ab, and analyzed by NGS. Frequencies of the reads mapped to each target gene relative to the total number of reads before (input, red) and after (IP, blue) IP are shown.

FIGURE 5.

A part of gene population downregulated by commensals is methylated at relatively high levels in their 5′ regions.

Genomic DNA was prepared from CoECs of CV mice by pooling from two independent experiments (four mice per experiment). The genomic DNA fragments corresponding to the sequences in nt −1000/+501 regions of 562 genes that were downregulated by commensals were captured using a customized array, immunoprecipitated with the anti-methylcytosine Ab, and analyzed by NGS. Frequencies of the reads mapped to each target gene relative to the total number of reads before (input, red) and after (IP, blue) IP are shown.

Close modal

Representative genes showing relatively high levels of methylation are indicated in Table VI. These include genes related to the cytoskeleton (Fgd1, Gas7, and Mark4) and apoptosis (Tnfrsf25). In addition, Table VII summarizes the methylation status in the 5′ region of the 20 genes listed in Table III. Out of these 20 genes, six genes showed relatively higher levels of methylation, nine were possibly methylated, and one was not methylated detectably. Among the six highly methylated genes, Ccl11 encoding the C-C chemokine ligand 11, which recruits eosinophils; Folr1 encoding folate receptor α, which is known to be overexpressed in various epithelial-derived cancers including colorectal cancer; and Slc2a5 encoding a fructose transporter were included. The DNA fragments corresponding to the 5′ regions of the remaining four genes were not detected in the input sample and were therefore considered not efficiently collected by the sequence capture array.

Table VI.
Representative genes with relatively high levels of DNA methylation in their 5′ region
GeneInput (%)IP (%)
FYVE, RhoGEF, and PH domain containing 1 (Fgd10.117 1.448 
TNFR superfamily, member 25 (Tnfrsf250.325 3.325 
Growth arrest specific 7 (Gas70.358 2.378 
MAP/microtubule affinity-regulating kinase 4 (Mark40.183 1.021 
Secreted and transmembrane 1A (Sectm1a0.316 1.476 
Src-like adaptor 2 (Sla20.291 1.270 
GeneInput (%)IP (%)
FYVE, RhoGEF, and PH domain containing 1 (Fgd10.117 1.448 
TNFR superfamily, member 25 (Tnfrsf250.325 3.325 
Growth arrest specific 7 (Gas70.358 2.378 
MAP/microtubule affinity-regulating kinase 4 (Mark40.183 1.021 
Secreted and transmembrane 1A (Sectm1a0.316 1.476 
Src-like adaptor 2 (Sla20.291 1.270 

Genomic DNA fragments corresponding to the 5′ regions of genes downregulated by commensals were captured using a customized array and then immunoprecipitated with anti-methylcytosine Ab. Representative genes with relatively high levels of DNA methylation are listed with percentage of reads obtained relative to the total reads for the input and IP samples.

Table VII.
DNA methylation status in the 5′ regions of genes downregulated by commensal bacteria
GeneInput (%)IP (%)
Relatively high methylation 
 Adenylate cyclase 9 (Adcy90.2331 0.6520 
 Carboxypeptidase A1 (Cpa10.2081 0.5687 
 C-C chemokine ligand 11 (Ccl110.4162 0.9115 
 Folate receptor 1 (Folr10.2081 0.4215 
 Solute carrier family 2 member 5 (Slc2a50.1748 0.2822 
 Calcium/calmodulin-dependent protein kinase II a (Camk2a0.3912 0.6104 
Possible methylation 
 IL-22R (Il22ra10.2913 0.1760 
 C-C chemokine ligand 25 (Ccl250.1665 0.0856 
 Inositol 1,4,5-trisphosphate 3-kinase A (Itpka0.1665 0.0353 
 Apolipoprotein A-I (Apoa10.2913 0.0463 
 Defensin a5 (Defa50.0250 0.0035 
 Cytochrome P450 (Cyp4b10.3246 0.0401 
 Fatty acid binding protein 4 (Fabp40.1332 0.0057 
 H+-K+-ATPase (Atp12a0.3163 0.0054 
 Defensin a24 (Defa240.0666 0.0005 
Methylation Not Detected 
 Adenosine A1R (Adora10.2580 — 
GeneInput (%)IP (%)
Relatively high methylation 
 Adenylate cyclase 9 (Adcy90.2331 0.6520 
 Carboxypeptidase A1 (Cpa10.2081 0.5687 
 C-C chemokine ligand 11 (Ccl110.4162 0.9115 
 Folate receptor 1 (Folr10.2081 0.4215 
 Solute carrier family 2 member 5 (Slc2a50.1748 0.2822 
 Calcium/calmodulin-dependent protein kinase II a (Camk2a0.3912 0.6104 
Possible methylation 
 IL-22R (Il22ra10.2913 0.1760 
 C-C chemokine ligand 25 (Ccl250.1665 0.0856 
 Inositol 1,4,5-trisphosphate 3-kinase A (Itpka0.1665 0.0353 
 Apolipoprotein A-I (Apoa10.2913 0.0463 
 Defensin a5 (Defa50.0250 0.0035 
 Cytochrome P450 (Cyp4b10.3246 0.0401 
 Fatty acid binding protein 4 (Fabp40.1332 0.0057 
 H+-K+-ATPase (Atp12a0.3163 0.0054 
 Defensin a24 (Defa240.0666 0.0005 
Methylation Not Detected 
 Adenosine A1R (Adora10.2580 — 

Results obtained from sequence capture NGS analysis for the status of methylation in the 5′ region of the genes downregulated by commensals (listed in Table III) are summarized. Percentage of reads obtained relative to the total reads for the IP sample [IP (%)] and the input sample [Input (%)] is shown. Relatively high methylation indicates that the proportion of reads relative to the total reads for the IP sample was higher than that for the input sample; possible methylation indicates that reads were obtained for the IP sample, but the proportion of the reads relative to the total reads for the IP sample was not higher than that for the input sample; and methylation not detected indicates that reads were obtained for the input sample but were not obtained for the IP sample. —, not detected.

Further, methylation levels in the 5′ region of the genes analyzed in this assay were compared between GF and CV mice by linking the data of MBD2 pull-down assay. The average ratio of methylation levels in GF mice to that in CV mice was 1.01 for the entire gene population analyzed in this assay. Thus, it can be inferred that no difference was found in the methylation levels between these mice. Next, we focused on the genes categorized as unmethylated and methylated by this assay. Out of the 35 unmethylated genes, data of methylation in GF and CV mice were obtained for 25 genes by MBD2 pull-down assay. When the GF versus CV ratio of methylation levels were compared between these 25 genes and the top 25 methylated genes selected in descending order of IP versus input ratios of read count percentages, they were significantly different (p < 0.05). The average ratios of GF to CV were 1.06 for the genes categorized as unmethylated and 0.964 for the genes categorized as methylated. Taken together, these results indicate that for a part of the gene population analyzed, DNA methylation is induced by commensal bacteria in the 5′ regions, showing that induction of DNA methylation is one mechanism by which commensal bacteria mediate the downregulation of genes in CoECs.

Comprehensive analysis of mRNA expression in CoECs by comparison of results from CV and GF mice showed that a greater number of genes was downregulated by commensal bacteria than those upregulated by them. In addition, a large proportion of the downregulated genes was suppressed by commensals in a MyD88-independent manner (Table I). These results suggest that the expression of a large number of genes in CoECs is suppressed by commensal bacteria through their metabolites rather than their cellular constituents. In fact, the expression of α-Defensin 5 in CoECs was suppressed by commensal bacteria almost completely in a MyD88-independent manner (Fig. 1), indicating that bacterial metabolites suppress expression of this gene in CoECs. Interestingly, it was shown that the expression of α-Defensin 5 gene was upregulated by commensals through their constituents in SIECs (Fig. 1). Given that commensal bacteria are found in greater abundance in the colon as compared with the small intestine and the composition of microbiota differs between these intestinal portions, CoECs and SIECs might have different regulatory mechanisms for genes that encode key molecules involved in host defense and immune responses.

In addition to immune-related genes, those encoding enzymes of the digestive system including amylase, glucosidase, and carboxypeptidase, fructose transporter, and lipid metabolism-related molecules such as apolipoprotein and fatty acid binding protein were also found to be downregulated in CoECs by commensals in a MyD88-independent manner (Table III), indicating that commensal bacteria might affect immune and inflammatory responses and also metabolism of the host through the regulation of gene expression in IECs. As digestion and absorption mainly occur in the small intestine rather than the colon, regulation of these genes in the small intestine needs to be analyzed. Other types of regulation by commensal bacteria than MyD88-independent downregulation of genes might also contribute to intestinal homeostasis. For example, MyD88-dependent induction of fucosyltransferase 2 (Fut2) in CoECs was observed (Table II) as previously reported (26). Fut2 fucosylates glycans on the surface of the intestinal epithelium and is known to serve as an attachment receptor and a nutrient source for specific bacteria, suggesting its role in establishing intestinal symbiosis and preventing diseases (2730).

The DNA methylation status of a specific population of genes was shown to be altered by gut microbiota. For instance, commensal bacteria-induced methylation of CpG motifs in the overlapping 5′ region of Tmem267 and 3110070M22Rik suppressed the expression of these genes (Fig. 3). Furthermore, these genes were found to regulate the expression of RALDH1; this is the first report, to our knowledge, showing the function of Tmem267 in IECs. The mechanism by which this transmembrane protein regulates Aldh1a1 expression is yet to be elucidated. Suppression of 3110070M22Rik upon RNAi decreased not only the expression of 3110070M22Rik but also, surprisingly, that of Tmem267 (Fig. 4). Although it is not clear why suppression of 3110070M22Rik leads to the decrease in Tmem267 expression, the effect of RNAi for 3110070M22Rik on Aldh1a1 expression is probably caused by the suppression of Tmem267 expression. Many studies have focused on RALDH2 expressed in dendritic cells and found that retinoic acid produced by this enzyme plays a crucial key role in IgA production, T cell homing, and regulatory T cell induction in the intestine (3133). In contrast, however, it has been also suggested that RALDH1, which is highly expressed in IECs, also affects the intestinal immune system because retinoic acid produced from IECs was shown to be important for dendritic cell conditioning (34, 35). It is considered that gut microbiota epigenetically controls specific host genes, Tmem267 and 3110070M22Rik, and thereby regulate RALDH1 expression in CoECs to control the concentration of retinoic acid, a key immunomodulator in the intestine. Interestingly, dietary fiber and short chain fatty acids (SCFAs) were shown to increase the expression of Aldh1a1 in SIECs (36). In particular, the effect was observed specifically in the proximal SIECs, not in CoECs, although gut microbiota is known to produce SCFAs by metabolizing dietary fiber mainly in the colon. Aldh1a1 expression in IECs was lower in CV mice than in GF mice in our study, indicating that gut microbiota might affect Aldh1a1 expression through both SCFA- and Tmem267-mediated manners differently in SIEC and CoECs. Recently, Grizotte-Lake et al. (37) reported that commensal bacteria suppress retinoic acid synthesis by IECs to control IL-22 activity and prevent dysbiosis. They showed that expression of Rdh7, another enzyme involved in the retinoic acid synthesis, is suppressed by commensals, suggesting that commensal bacteria regulate expression of multiple enzymes involved in vitamin A metabolism in IECs.

Focusing on the population of genes, whose expression was downregulated to <1/3-fold by commensals, a part of these genes was shown to be actually methylated (Fig. 5, Tables VI, VII). The remaining genes might be suppressed by commensals through mechanisms other than DNA methylation, such as the binding of some suppressive transcription factors. These findings suggest that commensal bacteria contribute to intestinal homeostasis via DNA methylation of specific genes in CoECs of the host. Thus, commensal bacteria maintain intestinal symbiosis by themselves through epigenetic control of a specific population of genes in CoECs, which might act as a mechanism supporting the intestinal ecosystem.

Only a few studies within the recent 1–2 y have comprehensively analyzed the status of DNA methylation in IECs in combination with mRNA expression. Howell et al. (38) evaluated differences in DNA methylation and transcription patterns in IECs between inflammatory bowel disease patients and controls by genome-wide DNA methylation and transcriptome analyses. They found that IECs from inflammatory bowel disease children before treatment showed alternations in DNA methylation and transcription and suggested that it might explain variations in disease outcomes and thus might be used as prognostic biomarkers. The correlation between clinical outcome and DNA methylation, or transcription, was analyzed not in combination, but separately in their study. In contrast, methylome and transcriptome analyses of SIECs from CV and GF mice were performed to examine changes during postnatal development (39). The results showed the presence of both microbiota-dependent and -independent processes and further helped to identify 126 genomic loci at which microbiota-dependent alternation of DNA methylation and transcription were both detected. In our study, integrated analysis of DNA methylation and transcription by focusing on DNA methylation profiles around the transcription start sites of genes in CoECs, several genes encoding key molecules for the maintenance of intestinal homeostasis such as cytoskeleton-related molecules, chemokine, receptor, and transporter were identified as those regulated by microbiota through DNA methylation (Tables VI, VII). Interestingly, it was reported that genomic DNA in CoECs of GF mice was hypermethylated as compared with that of CV mice (40). Although significant difference in methylation frequencies was not observed in regions around the transcription start sites (Fig. 2), and a greater number of genes were expressed at lower levels in CoECs of CV mice than in those of GF mice (Table I), how the genome-wide methylation status, including that of the regions apart from the transcription start sites, affects the function of IECs is yet to be elucidated. Our experiments focusing on DNA methylation in the 5′ regions of genes and mRNA expression in CoECs from CV and GF mice show that modulation of DNA methylation in the 5′ regions of host genes by commensal bacteria, as well as their metabolites, act as one mechanism for the regulation of gene expression in IECs. However, we cannot completely exclude the possibility that the difference of irradiated and nonirradiated diet between the GF and CV mice, respectively, affects the status of DNA methylation or gene expression. Further analyses will clarify detailed molecular mechanisms that explain how gut microbiota affects physiological functions of the host by regulating gene expression in IECs via DNA methylation, which leads to the establishment of targets and strategies for preventing diseases by interference with intestinal ecosystems.

We are very grateful to Dr. Ametani and Dr. Yurino (INFOBIO Co., Ltd.) for kind discussion of the NGS data analysis.

This study was supported in part by grants from the Japan Society for the Promotion of Science (KAKENHI 17K07801) and Nagase Science and Technology Foundation (to K.T.).

The sequences presented in this article have been submitted to DNA Data Bank of Japan Sequence Read Archive (https://www.ddbj.nig.ac.jp/dra) under accession numbers DRA008905, DRA008906, and DRA008907.

The online version of this article contains supplemental material.

Abbreviations used in this article:

     
  • CoEC

    colonic epithelial cell

  •  
  • CV

    conventional

  •  
  • DDBJ

    DNA Data Bank of Japan

  •  
  • GF

    germ-free

  •  
  • IEC

    intestinal epithelial cell

  •  
  • IP

    immunoprecipitation

  •  
  • LMD

    laser microdissection

  •  
  • MBD2

    methyl-CpG–binding domain protein 2

  •  
  • NGS

    next-generation sequencing

  •  
  • qRT-PCR

    quantitative RT-PCR

  •  
  • Reg3

    regenerating islet-derived 3

  •  
  • RNAi

    RNA interference

  •  
  • RNA-seq

    RNA sequencing

  •  
  • SCFA

    short chain fatty acid

  •  
  • SIEC

    small IEC

  •  
  • TE

    Tris–EDTA

  •  
  • WT

    wild-type.

1
Parker
A.
,
S.
Fonseca
,
S. R.
Carding
.
2020
.
Gut microbes and metabolites as modulators of blood-brain barrier integrity and brain health.
Gut Microbes
11
:
135
157
.
2
Zaiss
M. M.
,
R. M.
Jones
,
G.
Schett
,
R.
Pacifici
.
2019
.
The gut-bone axis: how bacterial metabolites bridge the distance.
J. Clin. Invest.
129
:
3018
3028
.
3
Stephen-Victor
E.
,
T. A.
Chatila
.
2019
.
Regulation of oral immune tolerance by the microbiome in food allergy.
Curr. Opin. Immunol.
60
:
141
147
.
4
Ducarmon
Q. R.
,
R. D.
Zwittink
,
B. V. H.
Hornung
,
W.
van Schaik
,
V. B.
Young
,
E. J.
Kuijper
.
2019
.
Gut microbiota and colonization resistance against bacterial enteric infection.
Microbiol. Mol. Biol. Rev.
83
: e00007-19.
5
Lo Presti
A.
,
F.
Zorzi
,
F.
Del Chierico
,
A.
Altomare
,
S.
Cocca
,
A.
Avola
,
F.
De Biasio
,
A.
Russo
,
E.
Cella
,
S.
Reddel
, et al
.
2019
.
Fecal and mucosal microbiota profiling in irritable bowel syndrome and inflammatory bowel disease.
Front. Microbiol.
10
:
1655
.
6
Herrera
S.
,
J.
Martínez-Sanz
,
S.
Serrano-Villar
.
2019
.
HIV, cancer, and the microbiota: common pathways influencing different diseases.
Front. Immunol.
10
:
1466
.
7
Sudo
N.
2019
.
Role of gut microbiota in brain function and stress-related pathology.
Biosci. Microbiota Food Health
38
:
75
80
.
8
Tilg
H.
,
N.
Zmora
,
T. E.
Adolph
,
E.
Elinav
.
2020
.
The intestinal microbiota fuelling metabolic inflammation.
Nat. Rev. Immunol.
20
:
40
54
.
9
Vatanen
T.
,
E. A.
Franzosa
,
R.
Schwager
,
S.
Tripathi
,
T. D.
Arthur
,
K.
Vehik
,
Å.
Lernmark
,
W. A.
Hagopian
,
M. J.
Rewers
,
J. X.
She
, et al
.
2018
.
The human gut microbiome in early-onset type 1 diabetes from the TEDDY study.
Nature
562
:
589
594
.
10
Tanoue
T.
,
S.
Morita
,
D. R.
Plichta
,
A. N.
Skelly
,
W.
Suda
,
Y.
Sugiura
,
S.
Narushima
,
H.
Vlamakis
,
I.
Motoo
,
K.
Sugita
, et al
.
2019
.
A defined commensal consortium elicits CD8 T cells and anti-cancer immunity.
Nature
565
:
600
605
.
11
Torres
J.
,
J.
Hu
,
A.
Seki
,
C.
Eisele
,
N.
Nair
,
R.
Huang
,
L.
Tarassishin
,
B.
Jharap
,
J.
Cote-Daigneault
,
Q.
Mao
, et al
.
2020
.
Infants born to mothers with IBD present with altered gut microbiome that transfers abnormalities of the adaptive immune system to germ-free mice.
Gut
69
:
42
51
.
12
Liao
L.
,
K. M.
Schneider
,
E. J. C.
Galvez
,
M.
Frissen
,
H. U.
Marschall
,
H.
Su
,
M.
Hatting
,
A.
Wahlström
,
J.
Haybaeck
,
P.
Puchas
, et al
.
2019
.
Intestinal dysbiosis augments liver disease progression via NLRP3 in a murine model of primary sclerosing cholangitis.
Gut
68
:
1477
1492
.
13
Allegretti
J. R.
,
B. H.
Mullish
,
C.
Kelly
,
M.
Fischer
.
2019
.
The evolution of the use of faecal microbiota transplantation and emerging therapeutic indications.
Lancet
394
:
420
431
.
14
Sanders
M. E.
,
D. J.
Merenstein
,
G.
Reid
,
G. R.
Gibson
,
R. A.
Rastall
.
2019
.
Probiotics and prebiotics in intestinal health and disease: from biology to the clinic. [Published erratum appears in 2019 Nat. Rev. Gastroenterol. Hepatol. 16: 642]
Nat. Rev. Gastroenterol. Hepatol.
16
:
605
616
.
15
Zhu
W.
,
M. G.
Winter
,
M. X.
Byndloss
,
L.
Spiga
,
B. A.
Duerkop
,
E. R.
Hughes
,
L.
Büttner
,
E.
de Lima Romão
,
C. L.
Behrendt
,
C. A.
Lopez
, et al
.
2018
.
Precision editing of the gut microbiota ameliorates colitis.
Nature
553
:
208
211
.
16
Depommier
C.
,
A.
Everard
,
C.
Druart
,
H.
Plovier
,
M.
Van Hul
,
S.
Vieira-Silva
,
G.
Falony
,
J.
Raes
,
D.
Maiter
,
N. M.
Delzenne
, et al
.
2019
.
Supplementation with Akkermansia muciniphila in overweight and obese human volunteers: a proof-of-concept exploratory study.
Nat. Med.
25
:
1096
1103
.
17
Zhu
W.
,
N.
Miyata
,
M. G.
Winter
,
A.
Arenales
,
E. R.
Hughes
,
L.
Spiga
,
J.
Kim
,
L.
Sifuentes-Dominguez
,
P.
Starokadomskyy
,
P.
Gopal
, et al
.
2019
.
Editing of the gut microbiota reduces carcinogenesis in mouse models of colitis-associated colorectal cancer.
J. Exp. Med.
216
:
2378
2393
.
18
Selma-Royo
M.
,
M.
Tarrazó
,
I.
García-Mantrana
,
C.
Gómez-Gallego
,
S.
Salminen
,
M. C.
Collado
.
2019
.
Shaping microbiota during the first 1000 days of life.
Adv. Exp. Med. Biol.
1125
:
3
24
.
19
Pronovost
G. N.
,
E. Y.
Hsiao
.
2019
.
Perinatal interactions between the microbiome, immunity, and neurodevelopment.
Immunity
50
:
18
36
.
20
Cerdó
T.
,
E.
Diéguez
,
C.
Campoy
.
2019
.
Early nutrition and gut microbiome: interrelationship between bacterial metabolism, immune system, brain structure, and neurodevelopment.
Am. J. Physiol. Endocrinol. Metab.
317
:
E617
E630
.
21
Woo
V.
,
E. M.
Eshleman
,
T.
Rice
,
J.
Whitt
,
B. A.
Vallance
,
T.
Alenghat
.
2019
.
Microbiota inhibit epithelial pathogen adherence by epigenetically regulating C-type lectin expression.
Front. Immunol.
10
:
928
.
22
Martin-Gallausiaux
C.
,
P.
Larraufie
,
A.
Jarry
,
F.
Béguet-Crespel
,
L.
Marinelli
,
F.
Ledue
,
F.
Reimann
,
H. M.
Blottière
,
N.
Lapaque
.
2018
.
Butyrate produced by commensal bacteria down-regulates indolamine 2,3-dioxygenase 1 (IDO-1) expression via a dual mechanism in human intestinal epithelial cells.
Front. Immunol.
9
:
2838
.
23
Takahashi
K.
,
Y.
Sugi
,
K.
Nakano
,
M.
Tsuda
,
K.
Kurihara
,
A.
Hosono
,
S.
Kaminogawa
.
2011
.
Epigenetic control of the host gene by commensal bacteria in large intestinal epithelial cells.
J. Biol. Chem.
286
:
35755
35762
.
24
Nakata
K.
,
Y.
Sugi
,
H.
Narabayashi
,
T.
Kobayakawa
,
Y.
Nakanishi
,
M.
Tsuda
,
A.
Hosono
,
S.
Kaminogawa
,
S.
Hanazawa
,
K.
Takahashi
.
2017
.
Commensal microbiota-induced microRNA modulates intestinal epithelial permeability through the small GTPase ARF4.
J. Biol. Chem.
292
:
15426
15433
.
25
Sugi
Y.
,
K.
Takahashi
,
K.
Kurihara
,
K.
Nakata
,
H.
Narabayashi
,
Y.
Hamamoto
,
M.
Suzuki
,
M.
Tsuda
,
S.
Hanazawa
,
A.
Hosono
,
S.
Kaminogawa
.
2016
.
Post-transcriptional regulation of Toll-interacting protein in the intestinal epithelium.
PLoS One
11
: e0164858.
26
Meng
D.
,
D. S.
Newburg
,
C.
Young
,
A.
Baker
,
S. L.
Tonkonogy
,
R. B.
Sartor
,
W. A.
Walker
,
N. N.
Nanthakumar
.
2007
.
Bacterial symbionts induce a FUT2-dependent fucosylated niche on colonic epithelium via ERK and JNK signaling.
Am. J. Physiol. Gastrointest. Liver Physiol.
293
:
G780
G787
.
27
Nanthakumar
N. N.
,
D.
Meng
,
D. S.
Newburg
.
2013
.
Glucocorticoids and microbiota regulate ontogeny of intestinal fucosyltransferase 2 requisite for gut homeostasis.
Glycobiology
23
:
1131
1141
.
28
Goto
Y.
,
T.
Obata
,
J.
Kunisawa
,
S.
Sato
,
I. I.
Ivanov
,
A.
Lamichhane
,
N.
Takeyama
,
M.
Kamioka
,
M.
Sakamoto
,
T.
Matsuki
, et al
.
2014
.
Innate lymphoid cells regulate intestinal epithelial cell glycosylation.
Science
345
: 1254009.
29
Pickard
J. M.
,
C. F.
Maurice
,
M. A.
Kinnebrew
,
M. C.
Abt
,
D.
Schenten
,
T. V.
Golovkina
,
S. R.
Bogatyrev
,
R. F.
Ismagilov
,
E. G.
Pamer
,
P. J.
Turnbaugh
,
A. V.
Chervonsky
.
2014
.
Rapid fucosylation of intestinal epithelium sustains host-commensal symbiosis in sickness.
Nature
514
:
638
641
.
30
Omata
Y.
,
R.
Aoki
,
A.
Aoki-Yoshida
,
K.
Hiemori
,
A.
Toyoda
,
H.
Tateno
,
C.
Suzuki
,
Y.
Takayama
.
2018
.
Reduced fucosylation in the distal intestinal epithelium of mice subjected to chronic social defeat stress.
Sci. Rep.
8
:
13199
.
31
Mora
J. R.
,
M.
Iwata
,
B.
Eksteen
,
S. Y.
Song
,
T.
Junt
,
B.
Senman
,
K. L.
Otipoby
,
A.
Yokota
,
H.
Takeuchi
,
P.
Ricciardi-Castagnoli
, et al
.
2006
.
Generation of gut-homing IgA-secreting B cells by intestinal dendritic cells.
Science
314
:
1157
1160
.
32
Iwata
M.
,
A.
Hirakiyama
,
Y.
Eshima
,
H.
Kagechika
,
C.
Kato
,
S. Y.
Song
.
2004
.
Retinoic acid imprints gut-homing specificity on T cells.
Immunity
21
:
527
538
.
33
Coombes
J. L.
,
K. R.
Siddiqui
,
C. V.
Arancibia-Cárcamo
,
J.
Hall
,
C. M.
Sun
,
Y.
Belkaid
,
F.
Powrie
.
2007
.
A functionally specialized population of mucosal CD103+ DCs induces Foxp3+ regulatory T cells via a TGF-beta and retinoic acid-dependent mechanism.
J. Exp. Med.
204
:
1757
1764
.
34
Edele
F.
,
R.
Molenaar
,
D.
Gütle
,
J. C.
Dudda
,
T.
Jakob
,
B.
Homey
,
R.
Mebius
,
M.
Hornef
,
S. F.
Martin
.
2008
.
Cutting edge: instructive role of peripheral tissue cells in the imprinting of T cell homing receptor patterns.
J. Immunol.
181
:
3745
3749
.
35
Iliev
I. D.
,
E.
Mileti
,
G.
Matteoli
,
M.
Chieppa
,
M.
Rescigno
.
2009
.
Intestinal epithelial cells promote colitis-protective regulatory T-cell differentiation through dendritic cell conditioning.
Mucosal Immunol.
2
:
340
350
.
36
Goverse
G.
,
R.
Molenaar
,
L.
Macia
,
J.
Tan
,
M. N.
Erkelens
,
T.
Konijn
,
M.
Knippenberg
,
E. C.
Cook
,
D.
Hanekamp
,
M.
Veldhoen
, et al
.
2017
.
Diet-derived short chain fatty acids stimulate intestinal epithelial cells to induce mucosal tolerogenic dendritic cells.
J. Immunol.
198
:
2172
2181
.
37
Grizotte-Lake
M.
,
G.
Zhong
,
K.
Duncan
,
J.
Kirkwood
,
N.
Iyer
,
I.
Smolenski
,
N.
Isoherranen
,
S.
Vaishnava
.
2018
.
Commensals suppress intestinal epithelial cell retinoic acid synthesis to regulate interleukin-22 activity and prevent microbial dysbiosis.
Immunity
49
:
1103
1115.e6
.
38
Howell
K. J.
,
J.
Kraiczy
,
K. M.
Nayak
,
M.
Gasparetto
,
A.
Ross
,
C.
Lee
,
T. N.
Mak
,
B. K.
Koo
,
N.
Kumar
,
T.
Lawley
, et al
.
2018
.
DNA methylation and transcription patterns in intestinal epithelial cells from pediatric patients with inflammatory bowel diseases differentiate disease subtypes and associate with outcome.
Gastroenterology
154
:
585
598
.
39
Pan
W. H.
,
F.
Sommer
,
M.
Falk-Paulsen
,
T.
Ulas
,
P.
Best
,
A.
Fazio
,
P.
Kachroo
,
A.
Luzius
,
M.
Jentzsch
,
A.
Rehman
, et al
.
2018
.
Exposure to the gut microbiota drives distinct methylome and transcriptome changes in intestinal epithelial cells during postnatal development.
Genome Med.
10
:
27
.
40
Poupeau
A.
,
C.
Garde
,
K.
Sulek
,
K.
Citirikkaya
,
J. T.
Treebak
,
M.
Arumugam
,
D.
Simar
,
L. E.
Olofsson
,
F.
Bäckhed
,
R.
Barrès
.
2019
.
Genes controlling the activation of natural killer lymphocytes are epigenetically remodeled in intestinal cells from germ-free mice.
FASEB J.
33
:
2719
2731
.

The authors have no financial conflicts of interest.

This article is distributed under the terms of the CC BY 4.0 Unported license.

Supplementary data