Altered Interactions between the Gut Microbiome and Colonic Mucosa Precede Polyposis in APCMin/+ Mice

Mutation of the adenomatous polyposis coli (APC gene), an early event in the adenoma-carcinoma sequence, is present in 70-80% of sporadic human colorectal adenomas and carcinomas. To test the hypothesis that mutation of the APC gene alters microbial interactions with host intestinal mucosa prior to the development of polyposis, culture-independent methods (targeted qPCR assays and Illumina sequencing of the 16S rRNA gene V1V2 hypervariable region) were used to compare the intestinal microbial composition of 30 six-week old C57BL/6 APCMin/+ and 30 congenic wild type (WT) mice. The results demonstrate that similar to 12-14 week old APCMin/+ mice with intestinal neoplasia, 6 week old APCMin/+ mice with no detectable neoplasia, exhibit an increased relative abundance of Bacteroidetes spp in the colon. Parallel mouse RNA sequence analysis, conducted on a subset of proximal colonic RNA samples (6 APCMin/+, 6 WT) revealed 130 differentially expressed genes (DEGs, fold change ≥ 2, FDR <0.05). Hierarchical clustering of the DEGs was carried out by using 1-r dissimilarity measurement, where r stands for the Pearson correlation, and Ward minimum variance linkage, in order to reduce the number of input variables. When the cluster centroids (medians) were included along with APC genotype as input variables in a negative binomial (NB) regression model, four of seven mouse gene clusters, in addition to APC genotype, were significantly associated with the increased relative abundance of Bacteroidetes spp. Three of the four clusters include several downregulated genes encoding immunoglobulin variable regions and non-protein coding RNAs. These results support the concept that mutation of the APC gene alters colonic-microbial interactions prior to polyposis. It remains to be determined whether interventions directed at ameliorating dysbiosis in APCMin/+mice, such as through probiotics, prebiotics or antibiotics, could reduce tumor formation.


Introduction
Alterations in the gut microbiome (dysbiosis) have been reported in human colonic neoplasia [1][2][3][4][5][6]. However it remains unclear as to whether dysbiosis represents a response to tumorigenesis or whether it precedes tumor formation. One of the most prominent genetic mutations associated with the pathogenesis of sporadic and hereditary colorectal cancers (CRC) lies in the tumor suppressing adenomatous polyposis coli (APC) gene [7][8][9][10][11][12][13]. A germ-line mutation of the APC gene causes familial adenomatous polyposis (FAP), which results in the development of multiple colorectal adenomas at an early age that unequivocally lead to CRC if no surgical interventions are taken. APC mutations also represent an early event in the adenoma-carcinoma sequence and are present in about 70-80% of sporadic human colorectal adenomas and carcinomas.
The multiple intestinal neoplasia (Min) mouse model of FAP carries a truncation mutation at codon 850 of the Apc gene [14]. Studies comparing the number of intestinal polyps in germfree and conventionally raised C57Bl/6 APC Min/+ mice suggest that the gut microbiome may promote development of intestinal neoplasia [15,16]. One study reported decreased incidence of polyps in only the mid small intestinal segment, however a subsequent study reported a significant reduction of intestinal adenomas in both the small and large intestine of germ-free mice compared with conventionally raised mice. Antibiotic treatment of C57BL/6 APC Min/+-MSH2 -/mice, which carry both the APC mutation and an HNPCC DNA mismatch repair mutation, reduced the number of polyps in both the small and large intestine [17].
We hypothesize that mutation of the APC gene results in alterations in host-microbiota interactions prior to tumor formation. To test this hypothesis, gut microbial composition was compared between 6 week-old C57Bl/6 APC Min/+ , prior to the development of detectable neoplasia [18], and congenic WT mice.

Animal Type and Housing
All of the mice were acclimated for two weeks in order to reduce stress from traveling. Carbon dioxide was used during euthanasia of the mice. This study was approved by the Institutional Animal Care and Use Committee (#202449) and Division of Laboratory Animal Resources at Stony Brook University. Three shipments of 10 four-week-old female C57BL/6J APC Min/+ and 10 four-week-old female C57BL/6J WT mice were received from The Jackson Laboratory (Bar Harbor, ME) between June 2012 and May 2013. APC Min/+ mice and WT mice were housed separately in groups of three to four in specific pathogen free (SPF) cages for two weeks prior to euthanization. All the experiments strictly followed guidelines from the Institutional Animal Care and Use Committee and Division of Laboratory Animal Resources at Stony Brook University.

Tissue and Luminal Content Sample Collection
All of the mice were euthanized at 6 weeks of age using carbon dioxide. Immediately after sacrifice, the gastrointestinal tract was divided along its cephalocaudal axis as previously described [19]. The segments analyzed included the ileum, cecum, proximal colon, and distal colon. Each small intestinal segment was washed in sterile phosphate buffered saline to remove the luminal content. A 1.0-1.5-cm section was obtained from the proximal ends of duodenum, jejunum, distal ends of ileum, proximal colon, and distal colon, and placed into RNAlater solution (Life Technologies, Grand Island, NY, USA) for RNA/DNA studies. The cecum was placed in its entirety in RNAlater. Three pellets of distal colonic luminal content (formed stool) were collected from the distal colon and stored in RNAlater. All the samples were kept at room temperature overnight and then stored in -80°C. In the first cohort, the remainder of each intestinal segment was processed into "swiss rolls" and fixed in 10% buffered formalin for histological analyses. In the second and third cohort, the remainder of each intestinal segment was stained with 0.25% methylene blue and inspected under a dissecting microscope for adenomas and aberrant crypt foci [20]. Stools were collected from nine 12-14 week APC Min\+ female mice with intestinal neoplasia [18] and six WT female mice, and placed into RNAlater.

DNA and RNA Extraction of Intestinal Tissue and Luminal Content Samples
Total RNA and DNA, (host and associated bacterial mixed community) were extracted from the duodenum, jejunum, ileum, cecal pouch, proximal colon and distal colon tissues using TRI Reagent (Sigma, St. Louis, MO) according to the manufacturer's recommendations. For distal colonic luminal content samples, DNA was extracted using the UltraClean Fecal Kit (Mo BIO Laboratories, Inc., Carlsbad, CA).
Quantitative PCR (qPCR) for targeted bacterial clades QPCR assays were performed using established primers for Bacteroides-Prevotella-Porphyromonas [21], Lachnospiraceae [21], and total bacteria [22] on all tissue and distal colonic luminal content samples as previously described [23,24]. The relative abundance of taxa within the Bacteroidetes phylum measured by ΔCt = Ct (threshold cycle) total bacteria -Ct Bacteroides-Prevotella-Porphyromonas . The relative abundance of taxa within the Lachnospiraceae clade (i.e., Clostridia Group Xi'an clade) was measured by ΔCt = Ct (threshold cycle) total bacteria -Ct Lachnospiraceae as previously described [23,24]. All assays were carried out in triplicate. Plasmid quantification standards were prepared from representative clones of the target organisms. The Mann-Whitney test using GraphPad Prism 5 (La Jolla, California) was performed to compare APC Min/+ and WT ΔCt values. The Bonferroni correction was made to correct for multiple comparison, thus significance required p< 0.025.

16S rRNA Amplicon Library Construction and Illumina V1V2 Sequencing Analysis
Bacterial profiles were determined by broad-range amplification and sequence analysis of 16S rRNA genes following our previously described methods [25,26]. In brief, amplicons were generated using primers that target approximately 300 bp. of the V1V2 variable region of the 16S rRNA gene (primers 8F and 338R, modified by addition of Illumina adapter and dual index sequences). PCR products were normalized using a SequalPrep kit (Invitrogen, Carlsbad, CA), pooled, lyophilized, purified and concentrated using a DNA Clean and Concentrator Kit (Zymo, Irvine, CA). Pooled amplicons was quantified using Qubit Fluorometer 2.0 (Invitrogen, Carlsbad, CA). The pool was diluted to 4nM and denatured with 0.2 N NaOH at room temperature. The denatured DNA was diluted to 15pM and spiked with 25% of the Illumina PhiX control DNA prior to loading the sequencer. Illumina paired-end sequencing was performed on the Miseq platform with version v2.3.0.8 of the Miseq Control Software and version v2.3.32 of MiSeq Reporter, using a 600 cycle version 3 reagent kit.
Paired-end sequences were sorted by sample via barcodes in the paired reads with a python script [25]. Sorted paired end sequence data were deposited in the NCBI Short Read Archive under BioProject Accession Number: PRJNA270112 (www.ncbi.nlm.nih.gov/bioproject/ PRJNA270112)PRJNA270112. The sorted paired reads were assembled using phrap [27,28] and paired reads that did not assemble were discarded. Assembled sequence ends were trimmed over a moving window of 5 nucleotides until average quality met or exceeded 20.
Trimmed contigs with more than 1 ambiguity or shorter than 200 nt were discarded. Potential chimeras identified with Uchime (usearch6.0.203_i86linux32) [29] using the Schloss [30] Silva reference sequences were removed from subsequent analyses. Assembled sequences were aligned and classified with SINA (1.2.11) using the 418,497 bacterial sequences in Silva 115NR99 as reference configured to yield the Silva taxonomy [31,32]. Operational taxonomic units (OTUs) were produced by clustering sequences with identical taxonomic assignments. OTU counts were normalized between samples by dividing sequence counts by the total number of sequences generated per sample. Phylum-level and family-level OTU tables were generated by collapsing lower level OTUs into higher-level categories. OTUs with a relative abundance <0.0001 and with a prevalence <0.01 were culled and relative abundances then transformed using the square root function [33]. The software package Explicet (v2.9.4, www. explicet.org) was used to display OTU data and estimate alpha diversity indices (i.e., S Chao1 , Shannon complexity [H], and Shannon Evenness [H/H max ]) through 1000 replicate samplings of rarefied datasets [34].
Comparisons of individual phyla and families passing the initial filtering step, were conducted as follows. Because of the commonly observed over-dispersion in microbiome count data [35], the effects of APC Min/+ genotype on individual OTU abundances were examined using the negative binomial (NB) regression model as follows: Y ijk denotes the OTU k's observed count for mouse j in shipment i, μ ijk is the mean parameter for OTU k's count distribution of mouse j in shipment I and ϕ k is the dispersion parameter. Shipment refers to three different deliveries of mice (mice 1-20, mice 21-40 and mice 41-60), which has a zero-mean random coefficient. Coefficients b 0 and β 1 are fixed constant representing grand mean and APC genotype respectively. The log total sequence count for each mouse is considered as an offset. In addition, the skewness of the distribution of low abundance OTUs causes a large proportion of zero counts. Therefore, a zero-inflated version of the negative binomial (NB) model is also fitted to OTUs with zero counts in additional to the NB model: The NB and zero-inflated NB models are chosen based on AIC criterion [36]. The p values for the genotype effects on each OTU were then adjusted by the Benjamini-Hochberg [37] procedure to calculate the FDR. Significance was set as FDR<0.05. In addition, a 10-fold cross-validation was performed to validate the significant OTUs that were identified. The square root transformation was then applied to the relative abundances to correct for the skewness and to reduce the coefficient of variation. Comparisons of overall microbial composition between APC Min/+ mice versus the wild type mice were subsequently conducted using the permutation Hotelling T2 test with 10,000 permutations using the R package 'Hotelling' [38]. Principle coordinate analysis (PCoA) was conducted at the lowest taxonomic level (genus) using the wcmdscale function implemented by the vegan R package [39] and using Morisita-Horn dissimilarity scores.

Histological Analysis of the Intestinal Sections
Histological analysis was carried out in a subset of 10 APC Min/+ mice and 10 wild type mice (1 st cohort) by constructing "Swiss rolls" of intestinal segments. These segments were stained with hematoxylin and eosin and scored for adenomas and inflammation by a pathologist who was blinded with respect to the genotype of the mice (N.O.), as previously described [40]. Inspection of intestinal segments stained with 0.25%methylene blue was carried out on the second and third cohorts of mice with the aid of a Zeiss dissecting scope for detection of adenomas and aberrant crypt foci [20].

Comparison of Mouse Proximal Colonic mRNA Expression in APCMin/ + and Wild type Mice
IL-1β mRNA expression relative to actin mRNA, was measured in cecal, proximal colonic, and distal colonic intestinal tissue RNA samples in all 30 APC Min/+ and 30 WT mice as previously described [8]. RNA extracted from the proximal colon of 9 week-old WT mice treated with 3% DSS in water for 7 days was used as a positive control for the assay. The IL-1β ΔCt values (ΔCt -= Ct actin -Ct IL1β ) were compared between APC Min/+ and WT groups using the Mann-Whitney test. Significance was set at a threshold of p <0.05.
Aliquots (1 μg) of proximal colon RNA samples from 6 APC were subjected to paired-ends 100 bp Illumina sequencing. The RNA-Seq libraries were prepared and sequenced at the New York Genome Center. Between 81 and 314 million reads were generated for each of the RNA samples. The RNA-Seq data were deposited in NCBI's Gene Expression Omnibus database with accession number GSE67634. The short reads were aligned to the GRCm38 genome (http://useast.ensembl.org/Mus_musculus/Info/Annotation) using STAR (Spliced Transcripts Alignment to a Reference) [41], and then converted to raw gene counts using featureCounts [42]. The edgeR package [43] was used to identify differentially expressed (FDR<0.05) genes (DEGs) between the APC Min/+ and wild type mice, using additional cutoff of 2-fold differential expression between groups. Hierarchical clustering based on the reads per kilobase of exon per million mapped reads (RPKM) value of the 130 DEGs was carried out by using 1-r dissimilarity measurement and Ward linkage, and the cluster number (n = 7) was chosen based on inspection of the coefficient of determination (R2) plot as previously described [44]. Second, a negative binomial (NB) regression model was fit with gene clusters as following: x ijh is gene cluster h's expressions centroids (medians) of mouse j in shipment i. μ ijk is the mean parameter for phylum k's count distribution of mouse j in shipment i and ϕ k is the dispersion parameter. "Shipment" is as defined above. b 0 and β 1 are fixed coefficients representing grand mean and APC genotype respectively. The log total count of each mice is considered as an offset. Best subset model selection was conducted to choose the model with lowest AIC. All models were fitted with R package: glmmADMB [45].

Results
The relative abundance of Bacteroidetes spp is increased in APC Min/+ mice colonic mucosa and luminal content prior to the development of intestinal neoplasias.
In preliminary targeted qPCR studies, we observed that the relative abundance of Bacteroidetes spp. in fecal DNA was higher in 12-14 week-old APC Min/+ female mice compared to age-matched WT female mice (ΔCt = -2.4 vs. ΔCt = -5.1, p = 0.0004). To test the hypothesis that the increase in the relative abundance of Bacteroidetes spp. preceded polyposis, we compared the relative abundance of this clade in the ileal, cecal, proximal colonic, distal colonic mucosa and the distal colonic luminal content in 6 week-old mice. Because gender effects have been previously reported on the number and location of polyps, we restricted our analysis to female mice [46]. Previous studies conducted on 6 week APC Min/+ female mice had established the absence of detectable neoplasias at that age [18]. The absence of intestinal adenomas and aberrant crypt foci was confirmed in the mice included in this study by microscopic inspection of the intestinal segments. The mean histological scores for inflammation were 0.1 and 0 (p = 0.37) for APC Min/+ and WT colons, respectively (n = 10 in each category). We observed a significant increase in the relative abundance of Bacteroidetes spp. in proximal colonic, distal colonic and distal luminal contents between the APC Min/+ and WT-control mice, but no significant difference in the ileal or cecal mucosal samples (Table 1). In contrast we observed no significant difference is Lachnospiriceae spp., a prominent group of Firmicutes, except in the luminal content of the distal colon, which exhibited significantly higher loads in WT mice (p = 0.003).
16S rRNA sequence analysis of proximal colonic mucosal samples from APC Min/+ and WT mice Illumina sequencing of the 16S rRNA gene V1V2 hypervariable region was carried out for the proximal colonic DNA samples. A total of 13,248,412 high-quality sequences were generated (average sequence length: 317 nt; average sample size: 220,807 sequences/sample; minimum: 9,381 sequences; maximum: 411,636 sequences, exclusive of negative controls which were near zero). The median Good's coverage score was 99.9987% at the rarefaction point of 9,381 sequences, indicating deep sequence coverage of the intestinal microbiome.
The 16S rRNA sequencing results confirmed the targeted qPCR results in demonstrating a significant increase in the relative abundance of taxa within the Bacteroidetes phylum (FDR = 0.0009) in the APC Min/+ mice (Fig 1A and Table 2). Within the Bacteroidetes phylum, the predominant family was S24-7, whose relative abundance was also significantly increased  week old APC Min/+ and 30 wild type mice. The p-values were carried out using the Mann-Whitney U test. The Bonferroni correction was made to the pvalue, so that significance was set at p<0.025.
in APC Min/+ mice ( Fig 1A and Table 3, FDR = 0.0015). 16S rRNA sequencing detected a significant reduction in the relative abundance of taxa within the Tenericutes phylum (FDR< 0.0001), as well as a significant reduction in the relative abundance of the Anaeroplasmataceae family (FDR < 0.0001), which was the most prevalent family in the Tenericutes phylum. The results were also confirmed in all 10 cross-validations. The relative abundance of the Cyanobacteria phylum and the Chloroplast family, which is the most prevalent family in this phylum), was decreased in APC Min/+ mice (FDR = 0.047). However this observation was confirmed in only 4 out of 10 cross-validations. Although 16S rRNA sequence analysis did not detect a significant difference in the overall relative abundance of the Firmicutes phylum, it also confirmed the lack of any change in the relative abundance of taxa within Lachnospiriciae clade, thus confirming the qPCR data (see Table 1). Furthermore decreases as well as increased were observed in some of the individual families within the Firmicutes phylum (Table 3). Similarly, although 16S rRNA sequence analysis did not detect a significant difference in the overall relative abundance of the Proteobacteria phylum, increases as well as decreases were observed in some of the individual families within the Proteobacteria phylum.
Overall, the proximal colonic microbiomes differed between APC Min/+ and WT mice at both the phylum (p<0.0001) and family (p <0.0001) levels, as assessed by a 10000-permutation Hotelling T2 test of the 7 phyla and 42 families that passed the initial filtering step (Fig 1A). The proximal colonic microbiomes of APC Min/+ mice were also characterized by significantly lower OTU complexity (Shannon H; p = 0.009) and evenness (Shannon H/Hmax; p = 0.006) compared with WT mice (Fig 1B); OTU richness (S Chao1 ) was comparable in the two groups. Finally, principle coordinates analysis (PCoA) demonstrated significant associations of principle component axes 1 and 2 with APC genotype (see Fig 2).

Differentially expressed genes in 6 week-old APC Min/+ mice
We reasoned that differences in colonic microbial composition associated with the APC Min/+ genotype must be linked to alterations in mouse colon gene expression, despite the lack of histological evidence of neoplasias. To examine how the APC mutation could alter the host colon gene expression, we conducted parallel RNA-sequence analysis on 6 APC Min/+ mice and 6 WT mice, (sampling all three cohorts). A total of 130 host genes (fold change > 2 fold, FDR <0.05) were selected using edgeR [43]. The DEGs were grouped into seven clusters as described in Methods (see Table 4), with 106 upregulated genes distributed among three clusters and the 24 downregulated genes distributed among four clusters.
NB regression (see Methods) selected four (1, 4, 6, 7) out of seven gene clusters (see Table 5) significantly (p-value < 0.05) in addition to APC genotype, which were positively associated with the relative abundance of Bacteroidetes. Among those four clusters, cluster 1 (coefficient 0.051) is composed of 100 out of 106 upregulated genes. The downregulated gene clusters 4,6,7 covers 22 out of 24 downregulated genes. While APC genotype had a dominant effect on the relative abundance of Bacteroidetes, detection of additional associations with mouse colonic gene expression, suggest that alterations in host colonic gene expression play a role in influencing mucosal associated microbial composition. Because elevated IL-1β levels have previously been reported in 18-25 week old APC Min/+ compared with WT mice, RT-PCR assays were conducted on proximal colonic RNA samples in all of the mice as previously described [15]. No significant difference was observed in the ΔCt IL1β-actin values between APC Min/+ and WT mice (-13.5 vs. -16.2, p = 0.485). These values were both very low compared to that measured in DSS treated mice (-4.9), indicating that IL-1β was not highly expressed in the colons of either mouse group in our study.

Discussion
Alterations in the gut microbiome (dysbiosis) have been reported in human colonic neoplasia and in mouse models [1][2][3][4][5][6]47]. This study demonstrates that alterations in the gut microbiome, characterized by an increased relative abundance of Bacteroidetes spp. observed in association with intestinal neoplasias, actually precedes the development of microscopically    [47], and in some but not all studies of human colorectal neoplasia [1][2][3][4][5][6]. 16S rRNA sequence analysis revealed that the increased relative abundance of Bacteroidetes spp. corresponded primarily to an increased relative abundance of taxa within the uncultured family S24-7. Similar increases in S24-7 have also been reported in conventionally raised C57BL6 mice that were fed a high fat diet [48]. This association with intestinal dysbiosis is of interest, because increased dietary fat has previously been associated with increased number and /or size in both WT and APC Min/+ mice [49,50]. The relative abundance of the phylum Tenericutes observed in this study is higher than reported by some studies of C57Bl/6 mice [51], but similar to another study using C57Bl/6 mice purchased from the same vendor [52]. In this study, the APC Min/+ and WT mice were  [53,54], possibly related to coprophagic behavior. Alterations in gene expression have been previously reported in normal appearing mucosa of APC mutant mice after the development of intestinal polyposis [55]. We report differential expression of genes (DEG) in 6 week APC Min/+ mice prior to the detection of intestinal polyposis. In order to integrate host colonic gene expression with the microbial taxonomic data, we reduced the gene expression input variables by first selecting DEGs, reasoning that these genes would be most likely to be involved in disrupted colonic microbial interactions in the mutant mice. Variable dimensionality was further reduced by clustering the 130 DEGs into seven groups. The detection of significant associations between host colonic gene expression and the relative abundance of microbial taxa, after taking into consideration APC genotype, support the concept that host colonic microbial cross talk influences mucosal associated microbial composition. Cluster 4, which included downregulated genes encoding immunoglobulins and non--coding functional RNAs, demonstrated a significant linear relationship with the relative abundance of Bacteroidetes after controlling for APC genotype. The inverse correlation between immunoglobulin gene expression and the relative abundance of Bacteroidetes spp. is intriguing in light of previous reports linking an increased relative abundance of Bacteroidetes spp. with a reduction of immunoglobulin coated bacteria in humans [56,57]. The observation that some of the non-coding RNAs in this cluster may be located in the mitochondria, is intriguing in light of the observation that mutated APC proteins in contrast to WT APC proteins are detected in mitochondria [58] In summary, our results support the concept that APC haplo-insufficiency of the host colonic epithelial cell alters colonic microbial interactions prior to polyposis. It is thus conceivable that such microbiome changes contribute to the pathogenesis of colon cancer. An important corollary to such a notion would be that the colonic microbiome represents an important (and druggable) target for the prevention of colon cancer. Indeed, interventions directed at the microbiome (germ free and antibiotic treatment) have been reported to modulate tumor formation in mouse models of colon cancer [17]. However, it remains to be determined whether interventions directed at ameliorating dysbiosis in APC Min/+ mice, such as through probiotic, prebiotic or antibiotic interventions, could reduce tumor formation.