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 (1–4). 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, 5–9). 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 (10–12). Therefore, many efforts have been made to prevent the onset or alleviate the symptoms of the diseases by intervention to the gut microbiota (13–17).
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 (18–20). 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 (21–23). 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.
Materials and Methods
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.
Methyl-CpG–binding domain protein 2 pull-down assay
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.
Sequence capture assay
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.
Commensal bacteria regulate gene expression in IECs in both MyD88-dependent and -independent manners
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).
Genes encoding antimicrobial peptides, cytokines, and chemokines are differently regulated by commensal bacteria in SIECs and CoECs
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.
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.
Commensal bacteria affect DNA methylation in the 5′ regions of specific genes in CoECs
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.
Aldh1a1 expression is regulated by commensal bacteria in CoECs via DNA methylation of Tmem267 and 3110070M22Rik genes
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).
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).
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.
A part of the 5′ regions of genes suppressed by commensal bacteria is highly methylated
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.
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.
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 (27–30).
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 (31–33). 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:
colonic epithelial cell
DNA Data Bank of Japan
intestinal epithelial cell
methyl-CpG–binding domain protein 2
regenerating islet-derived 3
short chain fatty acid
The authors have no financial conflicts of interest.