ABSTRACT
Enterotoxigenic Bacteroides fragilis (ETBF) is a Gram-negative, obligate anaerobe member of the gut microbial community in up to 40% of healthy individuals. This bacterium is found more frequently in people with colorectal cancer (CRC) and causes tumor formation in the distal colon of multiple intestinal neoplasia (Apcmin/+) mice; tumor formation is dependent on ETBF-secreted Bacteroides fragilis toxin (BFT). Because of the extensive data connecting alterations in the epigenome with tumor formation, initial experiments attempting to connect BFT-induced tumor formation with methylation in colon epithelial cells (CECs) have been performed, but the effect of BFT on other epigenetic processes, such as chromatin structure, remains unexplored. Here, the changes in gene expression (transcriptome sequencing [RNA-seq]) and chromatin accessibility (assay for transposase-accessible chromatin using sequencing) induced by treatment of HT29/C1 cells with BFT for 24 and 48 h were examined. Our data show that several genes are differentially expressed after BFT treatment and that these changes relate to the interaction between bacteria and CECs. Further, sites of increased chromatin accessibility are associated with the location of enhancers in CECs and the binding sites of transcription factors in the AP-1/ATF family; they are also enriched for common differentially methylated regions (DMRs) in CRC. These data provide insight into the mechanisms by which BFT induces tumor formation and lay the groundwork for future in vivo studies to explore the impact of BFT on nuclear structure and function.
INTRODUCTION
Bacteroides fragilis is a Gram-negative, obligate anaerobe that is found consistently, but in low numbers, in the gut microbial community (1). Enterotoxigenic B. fragilis (ETBF) is a particular subtype of Bacteroides fragilis, characterized by its production of Bacteroides fragilis toxin (BFT) and its association with diarrhea, inflammatory bowel disease, and colon cancer (2, 3). Since the identification of ETBF, many studies have examined the carriage rate of this potentially pathogenic bacterium. Some studies have probed the mucosa for the bacterium and found an asymptomatic carriage rate as high as 67% in those without colorectal cancer (CRC) (3), but most studies examine the stool and report rates of 10% to 12% (2). So, while we do not yet know exactly how many people are colonized with ETBF, a significant fraction of the population likely experiences prolonged exposure to the bacterium and, thus, BFT.
The bft gene codes for the BFT preproprotein that is processed by ETBF, yielding a secreted active 20-kDa zinc-dependent metalloprotease (2). There are three isotypes of the BFT protein, each of which is produced by a different bft gene (2), with the most potent of the BFT isotypes being BFT2. Expression of BFT2 is required for ETBF-induced tumorigenesis in multiple intestinal neoplasia (Apcmin/+) mice (4). When interacting with colon epithelial cells (CECs), BFT2 binds to an unidentified cell surface receptor, which leads to the rapid but indirect cleavage of E-cadherin, which, in turn, enhances barrier permeability in vivo (5). Because a pool of β-catenin is associated with the intracellular domain of E-cadherin, cleavage of E-cadherin triggers the release of β-catenin and its subsequent nuclear localization with the upregulation of c-Myc expression and enhanced cellular proliferation (6).
Though we know some of the immediate downstream effects of BFT on CECs, we still do not understand how this potent exotoxin causes lasting changes in CECs that lead to tumor formation and growth. A possible explanation is an epigenetic progenitor model to cancer, in which epigenetic dysregulation enables tumor cell proliferation and survival (7). Epigenetics, the study of heritable phenotypic variation that is not caused by base pair changes in the genome, has been linked to various cancers. In colorectal cancer specifically, hypermethylation of CpG islands and global hypomethylation have both been implicated in tumor development and progression (8), and global dysregulation of the methylome has been observed in clinical colon cancer samples (9).
Initial experiments exploring the impact of BFT on the epigenome have been promising. By tracking the localization of specific epigenetic regulators, one study showed that inoculation of C57BL/6J mice with ETBF causes an upregulation of gene silencing complexes at CpG islands (10). Subsequent studies have shown that inoculation of Apcmin/+ mice with ETBF also causes recruitment of DNA methyltransferase 1 (DNMT1), a recruitment potentially mediated by DNA mismatch repair proteins (11). While these studies have established a role for gene regulation and DNA methylation in the effect of in vivo ETBF infection on CECs, studies do not report on how ETBF affects other epigenetic processes, such as chromatin structure. This is of particular importance because chromatin state, as measured by changes in various histone marks, has been linked to colon cancer development (8), and differences in chromatin structure have been shown to impact mutation rates along the genome (12, 13).
Herein, we further tested the hypothesis that the epigenome is part of the mechanism by which BFT acts to increase colon tumor formation by using a reductionist model system to investigate the effect of BFT2 on CEC chromatin structure. Specifically, we examined whether BFT2 alters CEC gene expression and whether this correlates with altered chromatin accessibility using a model human colon epithelial cell line, HT29/C1, that is known to be uniformly sensitive to BFT2 and for which extensive prior data exist (2). We extended our analyses to determine if BFT2-induced functional gene expression changes and/or chromatin changes correlate with transcription factor binding motifs and known cancer-causing mutations in human CRC.
RESULTS
BFT2 treatment induces changes in gene expression that are more pronounced at 24 h than at 48 h.While the effect of BFT2 treatment on select genes at relatively early time points has been studied (14, 15), how this toxin affects genome-wide expression and whether these effects persist are still unknown. So, to examine how BFT2 treatment affects global gene expression in CECs, we performed transcriptome sequencing (RNA-seq) at time points of 24 and 48 h after a single BFT2 treatment in biological triplicate (Fig. 1A).
BFT2 induces changes in gene expression that relate to bacterial interaction with CECs. (A) HT29/C1 cells were treated with BFT2 or left untreated for 24 or 48 h. After each time point, ATAC-seq and RNA isolation were performed. (B) Heatmap showing the change in gene expression at 24 h after BFT2 treatment for those genes with q values of <0.1. RNA-seq was performed on untreated cells and those treated with BFT2 to determine if BFT2 alters gene expression in HT29/C1 cells. Colors represent log2 transcripts per million (log2TPM)-normalized values of gene expression in BFT2-treated samples and untreated samples. Extended lists of differentially expressed genes after BFT2 treatment are provided in Data Sets S1 and S2 in the supplemental material. (C) Bar graphs showing the expression of Ceacam1 and Muc2 after treatment of murine distal colon organoids with BFT2 (100 ng/ml). The expression of each gene after cell treatment was determined using TaqMan qPCR and normalized to Gapdh gene expression (also determined using TaqMan qPCR). Then, the fold change in gene expression was calculated by comparing gene expression before and after BFT2 treatment. Error bars represent the standard error from at least two replicate experiments.
Eight genes were found to be significantly differentially expressed at 24 h after BFT2 treatment (Fig. 1B; see also Data Set S1 in the supplemental material), and none were significantly differentially expressed at 48 h after BFT2 treatment. These data follow trends that have been previously reported, in which BFT2 has an effect on gene expression at very early time points that is less pronounced by 24 h and seems to return to baseline by 48 h after toxin treatment (6, 14, 16). Though the trends that we saw were reinforced by previous studies, the small number of genes that showed statistically significant gene expression changes (with multiple-testing correction) after BFT treatment was unexpected. We believe that the limited statistical significance seen was due to the relatively small changes in gene expression observed and the small sample size. This hypothesis was bolstered by the findings of our quantitative PCR (qPCR) analyses for a subset of genes, the results for several of which lacked statistical significance after correction for multiple testing. For all genes, the direction of the changes in gene expression matched the observations in the RNA-seq data (Fig. S1).
Hence, we decided to examine the genes that were no longer significantly differentially expressed after correction for multiple testing but possessed a P value of less than 0.01. With these criteria, after treatment with BFT2 for 24 h, we found 70 genes that were differentially expressed (Data Set S1). Of these genes, 41 showed a decrease in gene expression, while 29 showed an increase in gene expression. After BFT2 treatment for 48 h, we found 16 differentially expressed genes; 3 showed a decrease in gene expression, and 13 showed an increase in gene expression (Data Set S2). This result follows the previously identified trend that BFT2 has a more potent effect on CECs at 24 h after toxin treatment which diminishes by 48 h. Furthermore, when we expanded our list of potentially relevant genes affected by BFT2, many of the genes identified were related to bacterial interactions with CECs. To follow up on this, we tested the effect of BFT2 on the expression of two genes, Muc2 and Ceacam1, that may mediate CEC mucus production and bacterial cell adhesion when ETBF and pks-positive (pks+) Escherichia coli cocolonize mice in mouse models of CRC (17). We also wanted to determine if BFT2 had the same effect on nontransformed CECs as it had on a transformed CEC cell line, so we utilized organoids isolated from the mouse distal colon (colonoids). BFT2 treatment of colonoids isolated from Apcmin/+ mice for 24 h resulted in increased expression of Ceacam1 (the only mouse Ceacam expressed in the colon) and decreased expression of Muc2 (Fig. 1C). At 48 h after toxin treatment, Muc2 expression was still decreased and Ceacam1 expression was slightly decreased as well (Fig. 1C). To determine if these effects persisted, the colonoids were treated with BFT2 every 48 h for a total of three treatments. The cells were then grown for an additional 7 days before gene expression was assessed. At this 7-day time point, gene expression was unchanged for both Ceacam1 and Muc2 (Fig. 1C). These results suggest that purified BFT2 can alter the expression of Ceacam1 and Muc2 at relatively short time points in vitro in nontransformed CECs. These changes in expression persist only if the toxin is continually renewed, as is thought to happen in the in vivo state, where ETBF colonization is continuous (18).
BFT2 induces dynamic changes in chromatin accessibility throughout the genome.To assay chromatin state, we performed an assay for transposase-accessible chromatin using sequencing (ATAC-seq) on BFT2-treated and untreated cells at 24 and 48 h after BFT2 treatment. Briefly, ATAC-seq leverages the accessibility of chromatin to preferentially attack open chromatin using a transposase (Nextera; Illumina), resulting in a tagmentation reaction. Subsequent preferential amplification and downstream sequencing result in a distribution of sequencing reads which reveals areas of open versus closed chromatin by examination of the depth of coverage (i.e., the number of reads which align to each area of the genome) (19). Each sample contained at least 50 million reads after alignment, filtering, removal of duplicate reads, and removal of mitochondrial reads (Data Set S3).
To identify regions specifically affected by BFT2, we first combined all replicates (n = 3) into one fastq file and utilized the MACS2 algorithm (20) to identify peaks of chromatin accessibility (areas of open chromatin) across all samples. We then used the FeatureCounts package (21) to assess each sample’s depth of coverage in the list of peak regions output by MACS2 and applied the DESeq2 package (22) to identify significant differentially accessible regions between the treated and untreated samples. We classified these peaks into 2 types: BFT2-opened peaks (regions of the genome that were significantly more accessible in our BFT2-treated samples) and BFT2-closed peaks (regions of the genome that were significantly more accessible in our untreated samples). This analysis generated 157 BFT2-opened peaks and 142 BFT2-closed peaks at 24 h (Fig. 2A) and 2 BFT2-opened peaks and 1 BFT2-closed peak at 48 h after cell treatment (Fig. 2B). Although we expected the changes to be modest, the number of peaks generated was less than expected, particularly at 48 h after treatment. We think that this is due to the combination of the subtle changes in chromatin induced by BFT2 and the variability that we saw between replicates (Fig. S2), as well as the variability seen between our 24-h and 48-h baseline samples (Fig. S3).
Chromatin accessibility and gene expression correlate at baseline but not after BFT2 treatment. (A and B) MA plots showing chromatin accessibility changes after BFT2 treatment at 24 and 48 h. Each dot represents a peak region of chromatin accessibility output by MACS2. Peak regions in which chromatin accessibility differs significantly between BFT2-treated and untreated samples are colored in red, while peak regions that do not differ are colored in gray. (C) Violin plot showing the correlation between gene expression and chromatin accessibility at baseline. Genes were sorted by the presence of a peak in their promoter in the 24-h baseline sample. The P value was calculated by comparing the expression of genes with a promoter peak (open) to the expression of those without a promoter peak (closed) in the 24-h baseline sample. (D) Violin plot showing the correlation between gene expression and chromatin accessibility 24 h after BFT2 treatment. Genes were sorted by the presence of a promoter-associated BFT2-opened peak or BFT2-closed peak, identified via DESeq2. P values were calculated by comparing the log2(FC) in gene expression for genes with a promoter BFT2-opened peak to the gene expression for those with a promoter BFT2-closed peak. For all graphs, P values were calculated using the Wilcoxon rank-sum test.
To explore how BFT2 affects the peak distribution throughout the genome, we performed a detailed analysis on the BFT2-opened and BFT2-closed chromatin regions. When we analyzed the distribution of these peaks, we saw that the peak distribution was enriched in promoter regions compared to the distribution of promoters throughout the genome for our 24-h BFT2-opened and 24-h BFT2-closed peaks (Table 1). Our merged peaks data set showed a similar enrichment in promoters, suggesting that the distribution seen in our 24-h BFT2-opened and 24-h BFT2-closed peaks was not a specific attribute of BFT2-induced changes in chromatin accessibility; this enrichment was likely due to the inherent increase in chromatin accessibility at promoter regions resulting from active gene transcription (Table 1).
Distribution of ATAC-seq peaks throughout the genome
We then wanted to specifically interrogate the chromatin state at enhancers; enhancers have shifts in chromatin state associated with their function on proximal genes (23). Using ENCODE data on putative enhancer locations in HCT116 cells (a colorectal cancer cell line) (ENCODE accession number ENCSR931GTJ), we measured the overlap with our regions of altered chromatin accessibility. Fifty-seven enhancers overlapped the BFT2-opened 24-h peaks and only 9 overlapped the BFT2-closed 24-h peaks. The difference in overlap is striking, with 35% of our 159 total BFT2-opened peaks containing an enhancer and only 6% of the 142 total BFT2-closed peaks containing an enhancer (Fisher’s exact test, P = 0.00001).
Overall, our ATAC-seq analysis revealed that BFT2-induced changes in chromatin accessibility are transient and seem to wane by 48 h after toxin treatment in our in vitro assay. The few peaks identified at 48 h after BFT2 treatment were located in intronic or intergenic regions and did not overlap the peaks identified at 24 h after BFT2 treatment. Our data also show that BFT2-opened chromatin is more likely to overlap enhancers, that BFT2 treatment appears to be more likely to open chromatin than to close it, and that chromatin accessibility is affected across the genome by BFT2 (Table 1). We expect that the effect in promoter, intergenic, and coding regions is important, and thus, we next attempted to assess the potential functional significance of each.
Chromatin accessibility correlates with gene expression at baseline but not after BFT2 treatment.We next tested the hypothesis that BFT2-induced changes in chromatin accessibility correlate with BFT2-induced changes in gene expression. First, we classified genes by promoter chromatin accessibility status in our 24-h untreated (baseline) sample. For this, we used MACS2 to call peaks on each individual sample. Only promoter peaks were used because of the intuitive connection between open promoter chromatin and changes in gene expression. Those genes that contained an ATAC-seq peak overlapping their promoter region were categorized as “open,” and those without a peak were categorized as “closed.” We then extracted normalized gene expression values (transcripts per million [TPM]) for each of these genes to determine if the average gene expression values were different for the two groups. We found that, on average, gene expression was statistically significantly higher among genes that were open in our 24-h baseline sample than among those that were closed (P < 2E−16) (Fig. 2C).
Next, we examined the correlation between changes in chromatin accessibility and changes in gene expression after BFT2 treatment at the same time points. To do this, we focused on BFT2-opened and BFT2-closed peaks, identified using DESeq2, which overlapped gene promoter regions. We then calculated the logarithmic fold change [log2(FC)] expression values for each of these genes at 24 or 48 h after BFT2 treatment. Because of the small number of 48-h BFT2-opened and 48-h BFT2-closed peaks, we were unable to assess the correlation between chromatin accessibility and gene expression at 48 h after BFT2 treatment. At 24 h after BFT2 treatment, the median change in gene expression for both BFT2-opened genes and BFT2-closed genes was near zero (Fig. 2D). Moreover, there were no genes with a statistically significant change in promoter chromatin accessibility after BFT2 treatment and a statistically significant change in gene expression after BFT2 treatment. So, changes in gene expression induced by BFT2 do not seem to correlate with changes in chromatin accessibility induced by BFT2.
Treatment with BFT2 causes increased chromatin accessibility at transcription factor binding sites.Gene expression is regulated via multiple mechanisms, including the binding of transcription factors to enhancers and gene promoters. Therefore, we also wanted to test the hypothesis that BFT2-induced changes in chromatin accessibility impact transcription factor binding sites and, thus, gene expression. This can be done by looking for specific transcription factor motifs at sites of altered chromatin accessibility or by calculating the overlap of sites of altered chromatin accessibility with sites of transcription factor binding.
Given the gene expression changes for transcription factors and how chromatin accessibility can play a role in altered DNA-protein binding, we queried our BFT2-opened and BFT2-closed peaks, identified via DESeq2, for specific transcription factor motifs using the haystack_bio package (24). This analysis revealed several transcription factor motifs that were enriched in our samples at 24 h after BFT2 treatment (Table 2; Data Set S4). Notably, these included JUND, JDP2, FOSL1, JUNB, FOS, and FOS::JUN, many of which are downstream of mitogen-activated protein kinase (MAPK) pathways previously shown to be modulated by BFT2 (16, 25, 26).
Transcription factor binding motifs enriched after cell treatment with BFT2 for 24 h
In subsequent analyses, we specifically examined FOSL1, as its expression was upregulated at 24 h after BFT2 treatment. Though the increased gene expression was not statistically significant in our RNA-seq analysis, we confirmed its upregulation at 24 h after BFT2 treatment using qPCR (Fig. S1). For FOSL1 to augment gene expression, it needs to bind to an accessible DNA-binding site, so we next examined the FOSL1-regulated genes and classified them according to accessibility at 24 h after BFT2 treatment, as determined by ATAC-seq. Only 1 gene with a statistically significant change in gene expression after BFT2 treatment is regulated by FOSL1. This gene had an accessible promoter after BFT2 treatment and showed an increase in gene expression after BFT2 treatment. When we expanded our list of differentially expressed genes after BFT2 treatment to include those with a P value of <0.01, we noted that while several FOSL1-regulated genes were upregulated after BFT2 treatment, a nearly equal proportion was also downregulated (Data Sets S1 and S2). We found that 6 of 7 FOSL1-regulated genes with accessible DNA binding sites showed an increase in expression after BFT2 treatment and that only 1 of 6 genes without an accessible DNA binding site showed an increase in expression after BFT2 treatment (Data Sets S1 and S2). Consequently, we concluded that both expression and accessibility data are needed to correctly surmise that transcription factor upregulation, enrichment of transcription factor binding motifs, and the presence of open chromatin in the promoter region all contribute to BFT2-induced gene regulation.
Chromatin accessibility is associated with differential DNA methylation and DNA mutation.The previous analyses helped us to better understand how BFT2 alters chromatin accessibility in CECs and to connect chromatin changes with anticipated gene expression changes. However, they did not help us to understand how changes in chromatin accessibility may contribute to common DNA modifications found in CRC. To explore this, we looked for a correlation between BFT2-induced changes in chromatin accessibility and single nucleotide variants (SNVs) and sites of differential methylation (DMRs) in CRC. SNVs and DMRs were extracted from the COSMIC database (see Materials and Methods). We first calculated the proportion of peaks in our BFT2-opened and BFT2-closed 24-h peak sets that overlapped a common SNV or DMR. We then used Fisher’s exact probability test to determine if the proportion of peaks overlapping an SNV or a DMR differed after treatment of HT29/C1 cells with BFT2.
Correlating our chromatin accessibility changes with SNVs in CRC allowed us to explore the effects of BFT2-induced chromatin accessibility changes in coding regions of the genome where specific chromatin states have been linked to altered mutation rates in cancer (12). Examining DMRs allowed us to explore whether BFT2-induced chromatin accessibility is associated with hypermethylation of CpG islands in promoter regions or hypomethylation of regulatory intergenic regions.
Our SNV analysis revealed a slightly larger proportion of BFT2-opened peaks than BFT2-closed peaks overlapping common SNVs in CRC at 24 h after BFT2 treatment (9/157 versus 3/142), but the difference did not reach statistical significance (P = 0.14). Thus, our data do not support an increased frequency of SNVs commonly found in CRC in regions with increased chromatin accessibility after BFT2 treatment. Our DMR analysis also found that a larger proportion of BFT2-opened peaks than BFT2-closed peaks (5/157 vs 0/142) overlap common DMRs in CRC at 24 h after BFT2 treatment (P = 0.06). When we further parsed the data to analyze regions of hypermethylation and hypomethylation individually, we found that the association between DMRs and chromatin accessibility was due entirely to regions of hypomethylation, as no peaks overlapped hypermethylated DMRs. Together, analysis suggests that CRC DMRs are more likely to occur in chromatin opened by BFT2 than in chromatin closed by BFT2.
DISCUSSION
In previous in vitro studies, treatment of CECs with BFT2 has been shown to affect a number of cellular pathways. Expression changes for genes in these pathways tend to peak from 3 to 6 h and usually return to baseline by 24 h, due to suspected BFT2 degradation intracellularly, yet other data support the suggestion that significant CEC biological changes occur late (6, 14, 25). In this study, we sought to determine if global changes in gene expression and chromatin accessibility that would provide insight into BFT biology and its impact on nuclear structure and function could be identified. Our RNA-seq data corroborate previous observations about the acute onset of BFT2 action on CECs in vitro and in vivo (6, 14, 16). Overall, our chromatin accessibility data mirror the pattern seen in our RNA-seq data: BFT2 induces changes at 24 h, and the effects seem to wane by 48 h. Superficially, one could interpret this as a lack of a durable effect, but there are several pieces of data that suggest that this could have important downstream consequences. In vitro, BFT2 induces cell proliferation, which begins at about 48 h after treatment and lasts until 96 h after treatment (6). Additionally, ETBF-induced tumor formation in Apcmin/+ mice continues to progress even after ETBF elimination from the gut microbiome with cefoxitin treatment. Though antibiotic treatment reduces the number of tumors present, only 2 weeks of ETBF colonization is sufficient to observe gross tumor formation, and only 5 days of ETBF colonization is sufficient to observe microadenoma formation (27). In vivo, it is likely that other factors work together with BFT2 to extend its effects. Factors such as interleukin-17 (IL-17), which is necessary for tumor formation in an ETBF-induced mouse model of CRC (28), may work synergistically with BFT2 to cause changes in the cell that eventually lead to tumor formation. In vivo, CECs are also consistently in contact with ETBF and, thus, are likely exposed to BFT for prolonged periods of time. So, while we observed the transient effects of BFT2 on the expression of specific genes, our data showing the accessibility of specific transcription factor binding sites and the coordination with regions differentially methylated in CRC likely have important consequences in vivo.
We identify several genes that plausibly explain how BFT2 contributes to tumorigenesis by allowing other bacteria to invade the mucus layer and bind to CECs. Specifically, BFT2 upregulated CEACAM6 (as well as Ceacam1 in mice) and downregulated MUC2. CEACAM6 acts as a receptor for adherent-invasive E. coli, an E. coli phenotype implicated in ileal inflammatory bowel disease (29) as well as the development/progression of CRC in individuals with familial adenomatous polyposis (FAP) (17) and sporadic CRC (30). MUC2, the major MUC protein present in human colonic mucus, is notably downregulated after BFT2 treatment; this result is supported by recent studies which show that colonization with ETBF results in decreased mucus thickness in HT29-MTX-E12 cells. pks+ E. coli, on the other hand, is unable to induce mucus degradation in this model but shows enhanced mucosal colonization when coinoculated with ETBF in C57BL/6J mice (17). Both the decreased MUC2 expression and the increased CEACAM6 expression likely contribute to the increased mucosal colonization and increased DNA damage seen when pks+ E. coli cocolonizes with ETBF. Collectively, these results may begin to explain how ETBF works with other bacteria to promote tumor formation and growth.
Areas of increased chromatin accessibility after BFT2 treatment are enriched in binding motifs for many transcription factors belonging to the AP-1/ATF family of transcription factors, including FOSL1, JUN, and JDP. These transcription factors function downstream of the Jun N-terminal protein kinase (JNK) pathway, known to be activated by BFT2 (16, 25, 26). After BFT2 treatment, we expected to see increased expression of genes which are regulated by these particular transcription factors, but chromatin opening alone appeared to be insufficient to drive changes in gene expression of these JNK pathway genes. Hence, we looked for transcription factors that had both an enriched binding motif and increased gene expression after BFT2 treatment. FOSL1 was the primary candidate revealed by our analysis. We found that though the presence of a FOSL1 binding site alone was unable to predict differential gene expression after BFT2 treatment, when coupled to promoter chromatin accessibility, an excellent concordance with increased gene expression was detected. This finding exemplifies how integration of chromatin accessibility and gene expression data can be applied to understand how BFT2 and potentially other bacterial toxins act on CECs.
Because CRC usually results from a series of genetic mutations and can become more invasive via a series of methylation changes, any potential interaction between BFT2, chromatin accessibility, and common SNVs/DMRs in CRC is important. Herein, we found a BFT2-induced increase in chromatin accessibility at sites of common DMRs in CRC. Previously, O’Hagan and colleagues have shown that murine colonization with ETBF induces the binding of DNMT1 (DNA methyltransferase 1) at promoters of low-expressing genes, suggesting increased promoter methylation (10). In contrast, no studies to date have explored the connection between BFT and hypomethylation. We hypothesize that increased chromatin accessibility is, at the least, one step critical to chromatin methylation pattern changes in response to BFT2. Future bisulfite sequencing experiments should be conducted to definitively correlate sites of methylation change after BFT2 treatment with sites of altered chromatin accessibility.
Our study has several limitations that provide avenues for future experiments. Most notably, our experiments were performed on HT29/C1 cells, a colon carcinoma cell line. These cells contain several mutations that may lead them to respond to BFT2 differently than normal CECs. As a result, these experiments should be confirmed in a primary cell culture system or, ideally, in vivo either in a mouse model or from human samples with and without ETBF colonization. Our study, by design, was also performed with transient toxin exposure over only an early time window. In vitro studies with BFT best mimic the earliest in vivo events, whereby ETBF induces CEC morphology changes and colitis by ∼24 h after murine colonization (31). We hypothesize that genes with modest BFT2-induced increases in gene expression may lead to long-term effects on signaling pathways in the cell and may persist in vivo where B. fragilis is continually present in the colon. The persistent CEC exposure to BFT afforded by chronic ETBF colonization in vivo is difficult to reliably model in vitro, due in part to CEC uptake and degradation of BFT (32). This is further demonstrated by our colonoid experiments, where, even with the renewal of short-term BFT2 toxin exposure, we did not observe persistent effects of BFT2 on CECs. Thus, additional in vivo studies are warranted since ETBF murine colonization is persistent and associated with ongoing IL-17-dominant inflammation and CEC hyperplasia after 1 year in C57BL/6 mice (18, 31). Ultimately, we want to know which BFT-induced chromatin and gene expression changes persist and synergize with the carcinogenic impact of increased mucosal IL-17 induced by ETBF.
Our data provide additional insight into the effect of BFT2, a bacterial exotoxin linked to CRC pathogenesis, on CECs. We identified new associations between BFT2-induced chromatin accessibility, gene expression changes, and enhancer locations. We also correlated these changes with previously published data on BFT2-induced signal transduction and DNA modifications that may contribute to tumor formation or growth. Moving forward, we need to better understand how BFT2 affects the genome and epigenome of in vivo CECs to determine how ETBF increases the risk for colon tumorigenesis.
MATERIALS AND METHODS
Culture of HT29/C1 cells.HT29/C1 cells were cultured at 37°C in 10% CO2. Cells were plated at 20% confluence 4 days before addition of BFT2 toxin. At day 0, HT29/C1 cells were washed five times with phosphate-buffered saline (PBS) and placed in minimal medium (Dulbecco’s modified Eagle medium [DMEM]) without fetal bovine serum (FBS) or antibiotics, and toxin was added one time at a concentration of 100 ng/ml (5 nM). After toxin addition, cells were allowed to grow at 37°C in 10% CO2 for 24 or 48 h.
ATAC-seq transposition reaction.For all experiments, cell counts were obtained by trypsinizing cells in 0.5% trypsin for 10 min at 37°C and then counting the cells using a hemocytometer. Separate wells with cells were then scraped, and ∼50,000 cells were added to an Eppendorf tube. ATAC-seq was performed using the protocol outlined by Buenrostro et al. (33). Briefly, cells were washed with 50 μl cold 1× PBS buffer and then centrifuged at 500 × g for 5 min at 4°C. Cells were then resuspended in 50 μl cold lysis buffer [10 mM Tris-HCl (pH 7.4), 10 mM NaCl, 3 mM MgCl2, 0.1% octylphenoxy poly(ethyleneoxy)ethanol (IGEPAL CA-630)], gently lysed to preserve cell nuclei, and centrifuged again at 500 × g for 10 min at 4°C. Cells were then washed 3 times with wash buffer [lysis buffer without octylphenoxy poly(ethyleneoxy)ethanol]. After each wash, cells were centrifuged at 500 × g for 5 min at 4°C. Cell nuclei were then resuspended in the transposition reaction mix and incubated for 30 min at 37°C. DNA was purified using a Qiagen MinElute PCR purification kit and eluted in 10 μl elution buffer. DNA was stored at −20°C until fragments were amplified via PCR.
PCR amplification following transposition.All the DNA purified following transposition was PCR amplified. To do so, the following were combined in a 0.2-ml PCR tube: 10 μl transposed DNA, 10 μl nuclease-free H2O, 2.5 μl 25 μM custom Nextera PCR primer 1, 2.5 μl 25 μM custom Nextera barcoded PCR primer 2, and 25 μl NEBNext high-fidelity 2× PCR master mix. The components were amplified as follows: 1 cycle of 72°C for 5 min and 98°C for 30 s and 5 cycles of 98°C for 10 s, 63°C for 30 s, and 72°C for 1 min. After initial amplification, the number of additional cycles to be run was determined using qPCR. For this, the following were combined in a 0.2-ml PCR tube: 5 μl of previously PCR-amplified DNA, 4.5 μl of nuclease-free H2O, 0.25 μl of 25 μM custom Nextera PCR primer 1, 0.25 μl of 25 μM custom Nextera PCR primer 2, and 5 μl Kapa SYBR Fast qPCR master mix (2×) (total reaction volume, 15 μl). The components were amplified as follows: 1 cycle of 98°C for 30 s and 20 cycles of 98°C for 10 s, 63°C for 30 s, and 72°C for 1 min. To calculate the number of additional cycles of PCR needed, linear normalized reporter value (Rn) versus cycle number was plotted and the cycle number that corresponded with 1/3 of the maximum fluorescent intensity was determined. After the number of additional PCR cycles was determined, the remaining 45 μl of PCR product was run as follows: 1 cycle of 98°C for 30 s and n cycles (as determined via qPCR) of 98°C for 10 s, 63°C for 30 s, and 72°C for 1 min. Finally, the amplified library was purified using 0.9× AMPure beads at room temperature and eluted in 20 μl RNase-free H2O.
ATAC-seq data analysis.The 4 nM pooled ATAC-seq library was sequenced on an Illumina HiSeq 2500 sequencer using 50-bp paired-end sequencing. For each sample, three biological replicates were sequenced. After sequencing, the data were analyzed using a variation of the pipeline developed by the A. Kundaje lab, as outlined on ENCODE (34). Briefly, reads were trimmed, aligned, and filtered using the Bowtie2 program (35). Then, .bam alignment files for all samples were merged to create one consensus file. Using that consensus file, peaks were called using MACS2 (20). The FeatureCounts package (21) was used to sum the number of reads in each sample over the given genomic regions (peaks) specified by MACS2. DEseq2 was then used to determine which of these peaks had different levels of chromatin accessibility by comparing treated and untreated samples for each time point, with the P value threshold being 0.01.
We used the ChIPpeakAnno package to separate peaks into specific regions of the genome (promoter, coding, intergenic) (36). We also used the ChIPpeakAnno package to associate peaks located in promoter regions with their nearest downstream genes. To identify transcription factor binding motifs, we applied the haystack_bio package (24). For this analysis, CG correction was turned off. For the peaks with differential chromatin accessibility measured via DESeq2, we compared the BFT2-opened peak file and the BFT2-closed peak file for both the 24-h and 48-h analyses. A ratio of greater than 1 represents transcription factor binding motifs that are present more frequently in the BFT2-opened peak file, while a ratio of less than 1 represents transcription factor binding motifs that are present more frequently in the BFT2-closed peak file.
RNA-seq assay.Cells were washed once with PBS and collected using a cell scraper, and RNA was extracted from the cells using a Qiagen RNeasy minikit. After RNA extraction, the RNA pellet was flash frozen using 100% ethanol and dry ice. Samples were then stored at −80°C. For library preparation, samples were removed from storage, and mRNA was enriched using a NEBNext poly(A) mRNA magnetic isolation module. Afterwards, a nondirectional RNA-seq library was constructed using the NEBNext Ultra RNA library preparation kit from Illumina. The 2 nM pooled RNA-seq library was sequenced using the Illumina HiSeq sequence. For one sample, 50-bp paired-end sequencing was performed. For the other two samples, 50-bp single-end sequencing was performed.
RNA-seq data analysis.After sequencing, the kallisto program was used to perform pseudoalignment of the raw RNA-seq data (37). Then, the sleuth program was used to quantify gene expression and perform differential expression analyses (38). The Wald test was used to calculate differential expression.
SNV and DMR data extraction.Data on common single nucleotide variants (SNVs) and differentially methylated regions (DMRs) in colorectal cancer were extracted from the COSMIC database. For SNVs, we extracted COSMIC mutation data from genomic screens. This generated a table of coding point mutations identified in genome-wide screens. We then extracted those SNVs with a recorded primary site of “large_intestine.” For DMRs we extracted the COSMIC complete differential methylation data, which are TCGA level 3 methylation data from the International Cancer Genome Consortium (ICGC) portal. We then extracted those DMRs with a recorded primary site of “large_intestine.” We then split that data set into those DMRs with methylation status “H” and “L” to denote hypermethylated and hypomethylated regions, respectively. Both SNV and DMR data were extracted for genome version Hg38.
Colonoid isolation/growth.Animal experiments were approved by the Johns Hopkins University Institutional Animal Care and Use Committee. Colonoids were isolated from the distal colon of Apcmin/+ mice. For isolation, mice were dissected and the distal colon was extracted. The colon was then minced into 1-mm pieces and washed 5 times with 7 ml cold complete chelating solution (CCS); the formula has been described previously (39). The colon pieces, along with 3 ml CCS, were placed in 1 well of a 6-well plate along with 200 μl EDTA and shaken on a rotary shaker at 4°C for 1 h. The 3-ml solution was then combined with 5 ml CCS and 2 ml FBS, and the supernatant containing the crypts was removed and saved. This solution was then spun at 300 × g for 10 min at 4°C. The supernatant was removed; the pellet was resuspended in 10 ml of complete medium without growth factors (CMGF−) and spun again at 300 × g for 10 min at 4°C. The supernatant was again removed, and the crypt pellet was resuspended in Matrigel (25 μl per well of a 24-well plate) and plated. Colonoids were grown at 37°C in 5% CO2 in complete medium with growth factors (CMGF+). The formulas for CMGF− and CMGF+ solutions have been described previously (40).
Colonoid BFT2 treatment.For BFT2 treatment, colonoids were washed once with PBS. BFT2 (100 ng/ml) was added while colonoids were in CMGF− solution. After 6 h, the medium with toxin was removed and the CMGF+ medium was replenished. Colonoids were allowed to grow for 24 h or 48 h after toxin addition. For 2-week toxin treatments, colonoids were treated as described above every 48 h for a total of 3 treatments. After the third treatment, the colonoids were allowed to grow for 7 more days in CMGF+ medium before RNA extraction.
qPCR.For HT29/C1 cells, cells were washed once with PBS and collected using a cell scraper, and RNA was extracted from cells using a Qiagen RNeasy minikit. For colonoids, Matrigel was dissolved in 1 ml organoid harvesting solution. RNA was then extracted using a Qiagen RNeasy minikit. For each sample, first-strand cDNA synthesis was performed using 500 ng of the extracted RNA, 50 ng/μl of random primers, 10 mM deoxynucleoside triphosphates, 25 mM MgCl2, 0.1 M dithiothreitol, 10× reverse transcription buffer, and 200 U/μl SuperScript III reverse transcriptase. Expression of each gene after cell treatment was determined using TaqMan qPCR and normalized to GAPDH (glyceraldehyde-3-phosphate dehydrogenase) gene expression (also determined using TaqMan qPCR). Then, the fold change in gene expression was calculated by comparing gene expression before and after BFT2 treatment.
Statistics.To calculate the correlation between chromatin accessibility (ATAC-seq) and gene expression (RNA-seq) at baseline and after BFT2 treatment, the Wilcoxon rank-sum test was used with an alpha value of 0.05. To calculate the enrichment of transcription factor binding motifs in BFT2-treated and untreated samples, Fisher’s exact test with an alpha value of 0.01 was used. To determine if there was an association between the location of BFT2-treated peaks and common SNVs and DMRs in CRC, Fisher’s exact test was used, in which BFT2-opened peaks were compared to BFT2-closed peaks.
Accession number(s).The complete experimental data set was deposited in the Gene Expression Omnibus (GEO) database under accession number GSE113220.
ACKNOWLEDGMENTS
This work was supported by the Bloomberg Philanthropies (to C.L.S.), National Institutes of Health R01CA179440 (to C.L.S.), the Johns Hopkins Department of Medicine (to C.L.S.), and the Howard Hughes Medical Institute (to J.A.).
We thank K. D. Hansen for helpful discussions.
C.L.S. reports grant support from Bristol Myer Squibb; all other authors declare no competing interests.
FOOTNOTES
- Received 1 October 2018.
- Returned for modification 1 November 2018.
- Accepted 13 March 2019.
- Accepted manuscript posted online 18 March 2019.
Supplemental material for this article may be found at https://doi.org/10.1128/IAI.00447-18.
- Copyright © 2019 American Society for Microbiology.
