Altered intestinal microbiota-host mitochondria crosstalk in new onset Crohn's disease

Abstract

Intestinal microbial dysbiosis is associated with Crohn's disease (CD). However, the mechanisms leading to the chronic mucosal inflammation that characterizes this disease remain unclear. In this report, we use systems-level approaches to study the interactions between the gut microbiota and host in new-onset paediatric patients to evaluate causality and mechanisms of disease. We report an altered host proteome in CD patients indicative of impaired mitochondrial functions. In particular, mitochondrial proteins implicated in H₂ S detoxification are downregulated, while the relative abundance of H₂ S microbial producers is increased. Network correlation analysis reveals that Atopobium parvulum controls the central hub of H₂ S producers. A. parvulum induces pancolitis in colitis-susceptible interleukin-10-deficient mice and this phenotype requires the presence of the intestinal microbiota. Administrating the H₂ S scavenger bismuth mitigates A. parvulum-induced colitis in vivo. This study reveals that host-microbiota interactions are disturbed in CD and thus provides mechanistic insights into CD pathogenesis.

ARTICLE

Received 13 May 2016 | Accepted 27 Sep 2016 | Published 23 Nov 2016

DOI: 10.1038/ncomms13419 OPEN

Altered intestinal microbiotahost mitochondria crosstalk in new onset Crohns disease

Walid Mottawea1,2,3,*, Cheng-Kang Chiang1,2,*, Marcus Mhlbauer4,*, Amanda E. Starr1,2, James Butcher1,2, Turki Abujamel1,2, Shelley A. Deeke1,2, Annette Brandel1,2, Hu Zhou2,5, Shadi Shokralla6, Mehrdad Hajibabaei6, Ruth Singleton7, Eric I. Benchimol7,8,9, Christian Jobin10, David R. Mack7,8, Daniel Figeys1,2,11 & Alain Stintzi1,2

Intestinal microbial dysbiosis is associated with Crohns disease (CD). However, the mechanisms leading to the chronic mucosal inammation that characterizes this disease remain unclear. In this report, we use systems-level approaches to study the interactions between the gut microbiota and host in new-onset paediatric patients to evaluate causality and mechanisms of disease. We report an altered host proteome in CD patients indicative of impaired mitochondrial functions. In particular, mitochondrial proteins implicated in H2S

detoxication are downregulated, while the relative abundance of H2S microbial producers is increased. Network correlation analysis reveals that Atopobium parvulum controls the central hub of H2S producers. A. parvulum induces pancolitis in colitis-susceptible interleukin-10-decient mice and this phenotype requires the presence of the intestinal microbiota.

Administrating the H2S scavenger bismuth mitigates A. parvulum-induced colitis in vivo. This study reveals that hostmicrobiota interactions are disturbed in CD and thus provides mechanistic insights into CD pathogenesis.

1 Department of Biochemistry, Microbiology and Immunology, University of Ottawa, 451 Smyth Road, Ottawa, Ontario, Canada K1H 8M5. 2 Ottawa Institute of Systems Biology, University of Ottawa, 451 Smyth Road, Ottawa, Ontario, Canada K1H 8M5. 3 Department of Microbiology and Immunology, Faculty of Pharmacy, Mansoura University, Mansoura 35516, Egypt. 4 Department of Medicine, Division of Gastroenterology, Hepatology and Nutrition, University of Florida, Gainesville, Florida 32611, USA. 5 Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai 201203, China. 6 Biodiversity Institute of Ontario, Department of Integrative Biology, University of Guelph, Guelph, Ontario, Canada N1G 2W1. 7 Childrens Hospital of Eastern Ontario (CHEO) Inammatory Bowel Disease Centre and CHEO Research Institute, 451 Smyth Road, Ottawa, Ontario, Canada K1H 8L1. 8 Department of Pediatrics, University of Ottawa, 451 Smyth Road, Ottawa, Ontario, Canada K1H 8M5. 9 School of Epidemiology, Public Health and Preventive Medicine, University of Ottawa, 451 Smyth Road, Ottawa, Ontario, Canada K1H 8M5. 10 Department of Medicine, Department of Infectious Diseases and Pathology, University of Florida, Gainesville, Florida 32611, USA. 11 Department of Chemistry and Biomolecular Sciences, University of Ottawa, 451 Smyth Road, Ottawa, Ontario, Canada K1H 8M5. * These authors contributed equally to this work. Correspondence and requests for materials should be addressed to A.S. (email: mailto:[email protected]

Web End =astintzi@ mailto:[email protected]

Web End =uottawa.ca ) or to C.J. (email: mailto:[email protected]

Web End [email protected] ) or to D.R.M. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419

Inammatory bowel disease (IBD) is characterized by chronic and relapsing mucosal inammation of the gastrointestinal tract and comprises two main subtypes, Crohns disease (CD)

and ulcerative colitis (UC)1. A prevailing hypothesis is that IBD development is a consequence of aberrant interplay between an altered intestinal microbiota with the host2. The composition of the intestinal microbiota from IBD patients is typically characterized by decreased prevalence of protective microorganisms (that is, Clostridium XIVa and IV groups) and an expansion of detrimental bacteria (that is, Enterobacteriaceae and Fusobacteria)27. However, there is limited information on causality and mechanism of disease. Moreover, a fundamental unanswered question is whether the observed microbial dysbiosis is a cause or a consequence of inammation.

To address causality and microbial dysbiosis in CD pathogenesis, we conducted a systems-level study of the interaction between the intestinal microbiota and host at the mucosa luminal interface (MLI) in newly diagnosed paediatric patients. New-onset, paediatric subjects constitute an important population to study the role of the intestinal microbiota in IBD as there are no treatment effects on the disease process, co-morbidities are rare and there are fewer confounders as compared with adults. In addition, the inammatory phase of IBD characteristically predominates in younger subjects and thus, these patients are more likely to reveal the underlying mechanisms that promote early development of disease. To further identify key microbial drivers of inammation, we characterized the MLI microbial community at the time of diagnostic endoscopy. In contrast to stool microbiota, the MLI microbial community is in direct contact with the site of inammation and is therefore likely to be directly involved in the initiation and maintenance of the diseased state. To understand the interactions between the intestinal microbiome and the host, we characterized the host proteome and identied key associations between the host proteome and intestinal microbiota. We also validated the causative role of a microbe that predicts disease severity in a mouse model of colitis. The use of these approaches in the context of a phenotypically well-dened patient cohort has enabled us to elucidate a possible mechanism of disease and to identify key microbial drivers of inammation.

ResultsCharacteristics of the paediatric IBD gut microbiota. Although others have previously characterized the microbiota composition at the mucosa-luminal interface (MLI) in adult cohorts8,9, the MLI of new-onset paediatric IBD patients has not been previously reported. Given the intimate location of the MLI microbiota on the host intestinal epithelium, we hypothesized that studying the microbial composition at this location would provide direct information on how microbes affect inammation and gut physiological processes. To characterize the composition of the MLI microbiota in IBD, we collected the MLI from the ascending colon of new-onset paediatric IBD patients and controls by ushing sterile water onto the mucosa to dislodge the mucus layer from the mucosal epithelial cells followed by aspiration of these mucosal lavages through the colonoscope (Table 1 and Supplementary Table 1). We investigated the MLI microbiota using 16S rDNA-based Illumina sequencing on a total of 65 CD,21 UC patients and 42 control subjects with the remaining samples from our cohort (29 CD, 16 UC and 21 controls) used for independent validation with alternative quantication approaches (that is, Ion Torrent/454 sequencing and quantitative PCR (qPCR)). Microbial diversity was unchanged between controls and the CD/UC patients (Fig. 1a and Supplementary Figs 1 and 2) in contrast to the current model of decreased diversity in IBD

patients10, but in agreement with a previous study on new-onset paediatric IBD patients3. These results suggest that rst-onset paediatric patients have a distinct microbial prole as compared with adults and post-treatment cohorts.

Phylotypes differentially abundant between CD and controls or UC and controls were identied using a linear discriminant analysis effect size method11 and a zero-inated Gaussian mixture model (metagenomeSeq)12. These two statistical methods identied most of the same key microbes (Fig. 1b and Supplementary Fig. 3; and Supplementary Data 1 and 2). In addition, to ensure that our use of the V6 hypervariable region as opposed to the V4 region commonly used in microbiota studies did not unduly inuence our results, we sequenced the V4 and V6 regions in a subset of our samples using Ion Torrent sequencing. Importantly, while there are obviously unique taxa identied as being differentially abundant using each approach (Supplementary Fig. 4, Supplementary Data 3), many of the same key taxa are identied using each technique. Moreover, these results are highly concordant with those obtained by Gevers et al.3, who characterized the microbiota composition of new-onset paediatric CD patients using rectal mucosal biopsies instead of MLI samples. These include increased Clostridiales and Bacteroidales in controls and increased Enterobacteriaceae, Veillonella, Fusobacterium, Neisseria and Haemophilus in CD patients (Supplementary Data 2). The identication of the same key players highlights the usefulness of MLI as biological samples and the high quality of our inception cohort.

Numerous reports have cautioned that 16S rDNA results can be confounded by the presence of sequencing artifacts and/or contamination from the kits used during 16S library construction13. To quantify the levels of these artifacts and contaminants in our experimental setup, we puried genomic DNA from pure bacterial cultures and subsequently processed and analysed these samples as done for our MLI aspirates. These experiments revealed that our sequencing approach yields accurate results, with each pure bacterial sample primarily composed of a single dominant operational taxonomic unit (OTU) matching the genus of the sequenced bacteria with a median relative abundance of B99.1% (Supplementary Fig. 5).

Using this data set, we also determined that our median PCR/sequencing error rate was 0.007, which is comparable to those obtained in other studies14 (Supplementary Fig. 6). In addition, we determined that potential contamination from our sequencing reagents/kits accounted for a median of B0.6%

relative abundance. Importantly, the levels of sequencing artifacts and contamination quantied in these control experiments are low enough to be condently removed during our bioinformatic processing and statistical analyses (Methods section).

Key phylotypes correlating with CD inammation severity. To identify candidate causal microbes of CD pathogenesis, we examined changes in colonic microbiota composition at different stages of the disease from mild to moderate to severe inammation (as determined by the PCDAI (Pediatric Crohns Disease Activity Index) score). As taxonomy binning might mask important microbial dynamic behaviour, we identied key phylotypes correlating with CD inammation severity at the OTU level (dened by sequence clustering at 97% identity as a proxy of species-level resolution). Non-supervised multidimensional scaling analysis failed to separate the microbial communities as a function of disease severity likely due to high interpersonal microbial variation (Supplementary Fig. 7) and instead segregated by microbial diversity as reported by others3. Nevertheless, we assessed the association of individual OTUs with disease severity using generalized linear models and controlling for potential

2 NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419 ARTICLE

Table 1 | Subjects characteristics.

Control Crohns disease Ulcerative colitis Subjects 63 94 37

Gender (F/M) 33/30 34/60 18/19 Median age, y (IQ range) 14 (1016) 14 (11.7515) 15 (12.517) Disease activity (Inactive/mild/moderate/severe)* NA 3/19/25/47 1/4/16/16 Inammation locationw NA Upper GIT 1 E1 2

Ileum 28 E2 2 Colon 18 E3 5 Ileumcolon 47 E4 28

*Estimated based on PCDAI/PUCAI; Pediatric CD/UC Activity Index.wE1, ulcerative proctitis; E2, left-sided UC; E3, extensive UC; and E4, pancolitis.

confounders (gender and inammation status of the sampled site). A total of 161 OTUs were found to strongly correlate,(53 OTUs positively and 108 OTUs negatively), with disease severity (Fig. 1c; Supplementary Data 4). The OTUs negatively correlating with severity include major butyrate producers such as Blautia, Lachnospiraceae, Roseburia, Eubacterium rectale, Ruminococcus, Clostridium and Faecalibacterium. This decrease in butyrate producers was further conrmed using a qPCR assay targeting the butyryl-CoA:acetate CoA-transferase (BCoAT) gene, which encodes a key enzyme of the major pathway for butyrate production in the gut environment. As expected, CD patients had signicantly fewer copies of the BCoAT gene as compared with the control subjects indicating a decreased microbial capacity to produce butyrate (Fig. 1d). The OTUs found to positively correlate with severity include members of the order Bacteriodales, the families Enterobacteriaceae and Veillonellaceae, and the genera Atopobium, Fusobacterium, Leptotrichia, Sutterella, Vibrio, Parabacteroides, Prevotella, Peptostreptococcus, Peptococcus and Streptococcus. It is noteworthy that one-fourth of these OTUs (Atopobium, Fusobacterium, Veillonella, Prevotella, Streptoccocus and Leptotrichia) are known to produce H2S through the fermentation of sulfur-containing amino-acids15. We therefore evaluated whether the relative abundance of these H2S producers could predict CD severity using receiver-operating characteristic analysis. These H2S producers predicted severe inammation with an AUC of 0.74 (95% condence intervals 0.61 to 0.87; Supplementary Fig. 8) indicating that this group of microbes could be used to classify CD patients with severe inammation. Indeed, the relative abundance of these H2S producers correlated positively with the severity of inammation (Fig. 1e). H2S can also be produced by sulfate-reducing bacteria (SRB) but we found no evidence of a link between the SRB and CD patients as shown by a qPCR assay targeting the dissimilatory sulte reductase gene (Fig. 1f). All together, these results suggest a possible role for the H2S-producing bacteria in CD pathogenesis through the fermentation of sulfur-containing compounds and not through dissimilatory sulfate reduction.

Atopobium parvulum is a key network hub. To predict interactions among individual OTUs, we constructed a network of correlations between the abundance of the OTUs found to be associated with disease severity (Fig. 2a). The resulting network included 89 nodes (OTUs) and 341 edges (representing 312 co-occurrence and 29 co-exclusion interactions). Two major co-excluded modules of OTUs appeared in the network analysis separating OTUs on the basis of their relative abundance with respect to inammation (increased or decreased in CD as compared with control subjects). One module of co-occurring OTUs consisted of Vagoccocus (OTU-570199), Streptococcus (OTU-218959), OTU-535825 (Veillonellaceae) with A. parvulum

(OTU-529659) as the central hub of this module. The OTUs from this module were all found to be increased in CD. TheA. parvulum-module anti-correlated with the other major module, which is mostly composed of OTUs that are decreased in CD and belong to the class Clostridia. The relative abundance ofA. parvulum was validated by qPCR and found to be positively correlated with CD severity (Fig. 2b). The correlation betweenA. parvulum and CD was also conrmed by 454 pyrosequencing of a subset of samples followed by linear discriminant analysis effect size method analysis (Supplementary Fig. 9).

Functional alterations in the paediatric CD host proteome. To gain mechanistic insights into the role of microbes in CD severity, we conducted an unbiased, shotgun, quantitative proteomic analysis of mucosal biopsies from controls (n 10) and CD

(n 21) subjects with various levels of disease severity and tested

for associations between these hostproteomes and the intestinal microbes. A total 3,323 proteins were quantied using the superSILAC proteomics approach described. The majority of proteins were quantiable within a 10-fold tissue/super-SILAC ratio (85.7%: non-normalized L/H ratios, 89.3%: normalized L/H ratios) in all the samples (Supplementary Fig. 10a). Within the 475 proteins quantied by ratios outside the range of 10-fold in at least one sample (non-normalized L/H ratios), only 91 proteins had a median ratio above 10-fold (Supplementary Fig. 10b). This is comparable to a recently published superSILAC proteome study16 that showed that 83.7% of identied proteins were within a 10-fold ratio. There were 39 proteins identied in the light samples that did not have a heavy counterpart (Supplementary Fig. 10c). Of these, only one protein (Ribonuclease pancreatic; Uniprot ID: P07998) was detected in more than half of the samples without a corresponding SILAC-labelled peak from super-SILAC standard. Taken together, these results highlight the applicability of the reference proteome for use with the human biopsy sample proteomes.

Three hundred and twenty of the 3,323 quantied proteins were identied as differentially expressed by comparing the CD patients versus control (t-test with qo0.05). PCA of the proteomic proles grouped the control subjects away from the CD patients with the rst two components explaining 61.7% of the variation (Fig. 2c). KEGG pathway analysis, through the DAVID bioinformatics resources, revealed the top 10 enriched pathways (Fishers exact test Po0.01; Supplementary Fig. 11).

Gene ontology analysis identied the mitochondrial proteins as the major discriminant feature of all differentially expressed proteins (Fig. 2d). Most mitochondrial proteins (93.5%; 100/107 mitochondrial proteins) were found to be signicantly down-regulated in the CD patients as compared with controls (Fishers exact test Po0.05) with the extent of downregulation rising with increased severity (Supplementary Fig. 12). Interestingly, components of the mitochondrial hydrogen sulde detoxication

NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419

a b

8

Lachnospiraceae (f)

Bacteroidales (o) Bacteroidales (o)

Shannon index

6

Lachnospiraceae (f)

Lachnospira (g)

Neisseria (g) Neisseriales (o) Neisseriaceae (f)

Haemophilus (g) Photobacterium (g)

Ferrimonas (g) Ferrimonadaceae (f) Alteromonadales (o)

Pasteurellaceae (g) Pasteurellales (o)

Eggerthella (g) Peptostreptococcaceae (f)

Peptostreptococcus (g) Vibrio (g) Veillonella (g) Vibrionales (o) Vibrionaceae (f) Clostridiaceae (f)

Clostridium (g) Fusobacterium (g)

Fusobacteria (c) Fusobacteriaceae (f)

Fusobacteriales (o)

Fusobacteria (p) Veillonellaceae (f)

Erwinia (g) Gemellales (o)

Gemella (g) Gemellaceae (f) Enterobacteriaceae (f)

Enterobacteriales (o) Enterobacteriaceae (f) Gammaproteobacteria (c)

Proteobacteria (p)

4 2 0LDA score (Log10)

Control

Erysipelotrichi (c) Erysipelotrichaceae (f)

Erysipelotrichaceae (f)

Erysipelotrichales (o) Tenericutes (p)

4

2

0 Control

Mild

Moderate

Severe

c

Logfold change

3 2 1

0

1

2

3

2 4

CD

d e f

30

20

**

25

CT (BCoAT Ct- Univ 16S Ct)

**

**

20

10

Log 2(H 2S-producing bacteria)

CT (DSR Ct- Univ 16S Ct)

15

20

15

10

10

5

5

0 0 0

Control Control

CD

Mild

Moderate

Severe

Control

Mild

Moderate Severe

Figure 1 | Assessment of the microbiota composition at the mucosa-luminal interface of new-onset paediatric IBD. (a) Diversity of the intestinal microbiota in control and CD paediatric patients as a function of disease severity. Shannon index was calculated using 500,000 reads per sample. Crosses indicate the mean while horizontal lines indicate the median. (b) Histogram of linear discriminant analysis (LDA) effect size score for CD differentially abundant taxa compared with controls (n 65 and 42 for CD and controls, respectively); only OTUs meeting an LDA signicant threshold Z2 with a

Po0.05 (pairwise Wilcoxon test) are shown and are denoted with their lowest dened taxonomy. (c) Phylogenetic tree of the differentially abundant OTUs identied by metagenomeSeq analysis (fold change Z2 and Po0.05); an increasing red intensity indicates OTUs whose relative abundance increased, whereas an increasing blue intensity indicates OTUs whose relative abundance decreased in CD patients with severe inammation as compared with mild (outer circle) or severe as compared with moderate (inner circle); the maximum colour output for this gure was set at a Log2 value of 3. (d) The qPCR quantication of butyrate producing bacteria in control and CD microbiota. Butyryl-CoA: acetyl-CoA transferase gene (BCoAT); n 20 and 34 for control

and CD, respectively; statistical comparison by MannWhitney two-tailed test; **Po0.01). (e) Log2 sum of rareed reads assigned to taxa known to produce H2S through amino acid fermentation were plotted for the control and as a function of CD severity (n 39, 11, 11, 43 for control subjects and mild,

moderate and severe CD, respectively). (f) The qPCR quantication of sulfate reducing bacteria in control and CD patient microbiota as a function of disease severity. Dissimilatory sulde reductase gene (DSR); n 16, 9, 9 and 8 for control, mild, moderate and severe CD, respectively). (a,e,f) Statistical

comparison by KruskalWallis test using Dunns post hoc test and followed by a Bonferroni correction for the signicance level; *Po0.05; **Po0.01. Crosses indicate the mean while horizontal lines indicate the median.

4 NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419 ARTICLE

a

b

0125 ** **

Control

25

20

Ct (Aparv Ct-Univ 16S Ct)

10

15

5

0

Mild

Moderate

Severe

c

5

4

3

2

d

P value

PC2 (7.4%)

Mitochondrion Mitochondrial part Mitochondrial membrane

Mitochondrial envelope Mitochondrial inner membrane

Organelle inner membrane Organelle membrane

Organelle envelope

Envelope Mitochondrial lumen

25.08

18.1818.8116.3016.61

27.59

20.6920.69

12.23

< 0.001

1

< 0.001

< 0.001

0

< 0.001 < 0.001 < 0.001 < 0.001 < 0.001 < 0.001 < 0.001

1

2

3

415

10 5 0 5 10PC1 (54.3%)

Control Mild Moderate Severe

6.00

0.00 12.00 18.00 24.00 30.00

33.86

36.00

Percentage (%)

Figure 2 | Enrichment of A. parvulum and altered mitochondrial proteome dene the severity of Crohns disease. (a) Interaction network of differentially abundant OTUs; each node represents an OTU and are sized according to their number of interactions; each edge denotes a signicant co-exclusion (red) or co-occurence (grey) relationship between OTUs. Nodes are coloured by their number of signicant co-exclusions. (b) The qPCR quantication ofA. parvulum in control and CD microbiota as a function of disease severity. DCt was calculated by subtracting average Ct values of universal 16S rDNA from average Ct values of A. parvulum specic primers (n 18 for controls, nine for mild CD and seven for moderate and severe CD). **Po0.01 estimated using

KruskalWallis followed by Dunns post hoc test. Crosses indicate the mean while horizontal lines indicate the median. (c) Principal component analysis of the differentially expressed proteins among controls and CD patients categorized as a function of disease severity. (d) Functional annotation (cellular component) analysis of the differentially expressed proteins; the 10 most signicantly enriched functional groups (GO terms) are shown (Fishers exact test Po10 13). The illustrated P values are for classications that were signicantly enriched compared to the whole proteomic data set.

complex17, namely the sulfur dioxygenase (ETHE1), the thiosulfate sulfurtransferase (TST) and the components of complexes III and IV of the mitochondrial respiratory chain, were downregulated in the CD patients compared with controls (two-tailed MannWhitney test Po0.05). Secondary validation by qRTPCR conrmed the repression of tst, cytochrome c oxidase subunit IV (hcox41) and the sulde dehydrogenase genes (SQRDL) transcripts, all of which detoxify H2S, in CD and/or UC patients (Supplementary Fig. 13). These ndings indicate that the decreased abundance of these H2S-detoxication proteins is a hallmark of CD disease activity and possibly UC as well.

We next performed transkingdom correlation analysis between the levels of differentially expressed mitochondrial proteins and the OTUs that correlated with CD inammation severity (Fig. 3). In total, 46 OTUs displayed signicant associations with 96 host mitochondrial proteins with a Po0.05 in using both Kendall and Spearman correlations. Interestingly, OTUs that were found to be decreased in CD displayed strong

positive correlations with mitochondrial proteins while OTUs found to be increased mainly displayed negative correlations. Of particular note, OTUs displaying positive correlations are known butyrate producers. OTU-580521 (Ruminococcaceae), OTU-182190 (Lachnospiraceae), OTU-470114 (Clostridium leptum), OTU-182190 (Clostridium) and OTU-175682/109018 (Eubacterium rectale) were the top six OTUs positively associated with mitochondrial proteins expression. Among the OTUs displaying negative correlations, OTU-529659 (A. parvulum), OTU-570199 (Vagococcus), OTU-158660 (Bacteroides), OTU-353985 (Parabacteroides distasonis) and OTU-64396 (Fusobacterium) were identied as the top candidates exhibiting reverse correlations with up to 31 mitochondrial proteins including proteins from the respiratory chain and the H2S

detoxication complex.

A. parvulum causes pancolitis in Il10 / mice. Altogether these observations suggest that A. parvulum might play a key role

NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419

Figure 3 | Microbiotamitochondria correlation analysis. Clustered heatmap of SpearmanKendall correlation analysis between differentially abundant OTUs and mitochondrial proteins; red and blue colours indicate a positive or negative Spearman correlation respectively (P value o0.05 for both the

Kendall Tau and the Spearmans rank statistics was used to dene signicance denoted by ). Note: only OTUs and mitochondrial proteins with at least

one signicant correlation are shown. Whether the OTU was increased or decreased in mild versus severe CD patients is denoted on the bar to the left of the OTU name (blue decreased, red increased).

in the development or progression of inammation severity in the CD patients. To evaluate the colitogenic potential of A. parvulum, we utilized colitis-susceptible Il10 / mice18,19. Germ-free

Il10 / ;NF-kBEGFP mice were transferred to specic pathogen free (SPF) housing and after 2 weeks acclimatization, gavaged with A. parvulum (108 colony-forming units (CFUs) per mouse). Thereafter the mice were gavaged with A. parvulum weekly for 6 weeks and killed 6 weeks later. Compared with control uninfected Il10 / ;NF-kBEGFP mice, A. parvulum-colonized Il10 / mice displayed colitis (crypt hyperplasia, ulcers, goblet cell depletion and immune cell inltration) with histologic inammation scores signicantly higher compared with controls (MannWhitney U-test Po0.05; Fig. 4a,b). Colonoscopy imaging revealed mucosal erythema, friability and mucosal ulceration in A. parvulum-colonized Il10 / mice compared with the healthy mucosa observed in control mice (Fig. 4c). At the molecular level, the colon of A. parvulum-colonized Il10 / mice showed increased expression of Cxcl1 and Il17 compared with controls (MannWhitney U-test Po0.01; Supplementary Fig. 14). To determine whether the colitogenic effect of A. parvulum was dependent on the presence of commensal bacteria, we mono-associated germ-free Il10 / mice with A. parvulum for 8 weeks.

Interestingly, A. parvulum mono-associated GF Il10 / mice did not cause signicant colitis compared with controls suggesting that this bacterium cannot induce disease on its own (Fig. 4d; Histological score: 1.6 versus 0.67 respectively; not signicant).

Bismuth (III)-subsalicylate (hereafter denoted as bismuth) is a known H2S scavenger20. To test the potential role of H2S in

A. parvulum-mediated colitis in Il10 / mice, gnotobiotic Il10 / mice were transferred to SPF conditions and randomized into four groups; SPF only, SPF plus A. parvulum,

SPF plus bismuth and nally, SPF plus A. parvulum and bismuth.

As expected A. parvulum worsened the development of colitis (Fig. 5ac). Interestingly, treatment with bismuth signicantly improved the colitis scores in A. parvulum-colonized mice (Fig. 5ac) and prevented expression of Il1b, Il12p40 and Cxcl1 messenger RNA (mRNA; Supplementary Fig. 15). Although these ndings suggest a role for H2S production in A. parvulum-induced colitis, the reduced colitis observed with bismuth administration could also be due to bismuths potential antimicrobial activity on A. parvulum and/or the intestinal microbiota.

DiscussionIn this study, we found key alterations in the microbiotahost interplay in newly diagnosed children with IBD providing insight into the mechanism of IBD pathogenesis and potentially opening new avenues for treatment interventions. Our study found that mitochondrial proteins were the primary proteins downregulated in CD indicating a central role for the mitochondria in CD pathogenesis. The signicant downregulation of these proteins was associated with a depletion of butyrate producers together with a bloom of pathobionts many of which are known to be potent H2S producers. These results suggest a regulatory relationship between the mitochondria and microbiota and a disturbance of this relationship in the CD patients.

Accumulating data support a role for mitochondrial dysfunction in IBD pathogenesis. Reduced mitochondrial functions have been reported in UC patients with an up to 60% decreased activity of the respiratory chain complexes II, III and IV (ref. 21). Morphological changes in the mitochondria have been observed in the IBD patients2224 and functional defects at complexes III and IV have been described in a paediatric CD patient25. Variants in mitochondrial DNA, which result in increased metabolic activities, protect mice from colitis26. The decreased expression

6 NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419 ARTICLE

a c

A. parvulum

A. parvulum

+A. parvulum

b

20

*

CecumDistal colon

15

A. parvulum

Inflammation score

10

500 m

500 m

5

+A. parvulum

0

A. parvulum

+A. parvulum

500 m

500 m

d

+A. parvulum

GF

500 m 100 m

20X

Figure 4 | Atopobium parvulum induces rapid and severe pan-colitis in Il10 / mice. (a) Representative histological sections of the cecum and colon of Il10 / mice under specic pathogen free (SPF) conditions that were either associated with ( ) or without ( ) A. parvulum. (b) Blinded histological

score of inammation (n 6 to 7 per group; crosses indicate the mean while horizontal lines indicate the median; comparison by MannWhitney two-tailed

test; *Po0.05). (c) Macroscopic monitoring of inammation with a murine endoscope (data is representative of two animals). (d) Representative histological sections of the distal colon of germ-free (GF) Il10 / mice mono-associated with or without A. parvulum.

of the mitochondrial proteins observed here in paediatric CD patients would perturb mitochondrial functions. Such mitochondrial perturbations have been shown to result in increased epithelial permeability through ROS production and promote transcytosis of bacteria across the epithelial layer27. Increased intestinal permeability has been reported in active CD further suggesting a role for mitochondrial dysfunction in the pathogenesis of IBD28,29.

Butyrate is the main energy source of colonocytes and is produced by the bacterial fermentation of unabsorbed carbohydrate. Here, we demonstrate that the intestinal microbiota of paediatric CD patients is characterized by a signicant decrease in the relative abundance of butyrate producers as has been reported previously in the context of paediatric IBD3. Our data also revealed a positive correlation between the relative abundance of the butyrate producers and mitochondrial proteins, suggesting a signalling role for butyrate in host mitochondrial gene expression. In support of this observation, butyrate is known to activate the expression of the genes encoding the host mitochondrial H2S

detoxication components30 and our proteomic analyses indicate a diminished capacity for H2S detoxication by paediatric

CD patients. Therefore, we postulate that the depletion of butyrate-producing microbes from the colonic microbiota may dampen host H2S defense systems. This impaired host would be highly susceptible to further damage caused by enhanced H2S production, resulting in additional metabolic stress and subsequently increased mucosal inammation. Notably, this decrease in butyrate producers is associated with a signicant

increase in pathobionts, many of which are known to produce H2S. These include Fusobacterium, Veillonella and Atopobium parvulum whose relative abundance correlated negatively with the expression of the mitochondrial proteins. Butyrate exhibits anti-inammatory properties, whereas H2S is an important mediator of many physiological and pathological processes and has been previously associated with UC and colorectal cancer17,31. Interestingly, A. parvulum was identied by network analysis as the most prominent microbe associated with mitochondrial dysfunction. Importantly, the increased abundance of A. parvulum as a function of inammation severity in paediatric CD patients was not observed in UC and therefore is not simply a consequence of inammation, suggesting a causal effect in CD.

Although causal relationships would be better identied by monitoring changes in microbiota composition before the onset of the disease while following disease establishment, logistical constraints preclude this experiment. Instead we tested the hypothesis that changes in microbial composition between distinct and progressive stages of the disease would unravel causal interactions among the host and specic members of the gut microbes. Indeed, the observation that the relative abundance of A. parvulum correlates positively with inammation severity and negatively with both the mitochondrial protein expression and the level of butyrate producers prompted us to assess the colitogenic potential of A. parvulum. Interestingly, we found thatA. parvulum triggers the development of colitis in Il10 / mice under specic pathogen-free conditions but not under mono-

NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419

a

A. parvulum +A. parvulum

+BismuthBismuth

500 m 100 m

500 m 100 m

20X

20X

b c

A. parvulum +A. parvulum

10

*** *

8

+BismuthBismuth

Inflammation score

6

4

2

0

A. parv

Bismuth

+ +

+ +

Figure 5 | Induction of A. parvulum-associated colitis requires the gut microbiota and hydrogen sulde. (a) Representative histological Swiss-rolled sections of the colon of Il10 / mice associated ( ) or not associated ( ) with A. parvulum, treated ( ) or not treated ( ) with bismuth and kept

under specic pathogen free (SPF) conditions for 14 weeks. (b) Representative murine endoscopies of SPF Il10 / mice associated ( ) or not-associated

( ) with A. parvulum, treated ( ) or not-treated ( ) with bismuth. (c) Blinded inammation scores (n 7 to 8 per group) for Il10 / mice under SPF

conditions; crosses indicate the mean while horizontal lines indicate the median; two-tailed MannWhitney test was used for statistical comparison; *Po0.05, ***Po0.001.

association conditions. These results indicate a key role for bothA. parvulum and the gut microbiota in colitis development. In support to our ndings, a recent report on CD patient response to exclusive enteral nutrition (EEN) therapy has reported that a putative A. parvulum OTU was enriched in CD patients before treatment and was negatively associated with days on EEN32. Although a role for H2S production by A. parvulum in colitis development remains to be proven directly, excess H2S has been recently shown to act as an autocrine T-cell activator33, potentially contributing to unwanted T-cell responses against commensal bacteria, consistent with our observation that the intestinal microbiota is required for A. parvulum-induced experimental colitis. Accordingly, potential H2S trapping by bismuth may have contributed to the mitigation of Atopobium-induced colitis development in our mouse model. This nding is in agreement with a prospective mechanistic role for H2S in inammation. A role for H2S in colitis development is supported by the recent nding that microbial H2S formation contributes to mucus degradation opening the intestinal barrier to toxic compounds and pathobionts34.

H2S is an important mediator of many biological processes (angiogenesis, cytoprotection, metabolism, inammation) and can be pro- or anti-inammatory depending on its concentration and the particular circumstance35,36. In light of these opposite effects, a recent study reported that paediatric CD patients undergoing EEN treatment and presenting with reduced colonic inammation had increased levels of fecal H2S37. The ndings from the latter study challenge the concept of a pro-inammatory role for H2S in CD. However, as EEN feed is enriched in sulfur-containing amino-acids, the observed increase in H2S levels might simply be a consequence of a change in diet. Clearly, complex modes of action of H2S need to be further explored to reconcile these apparent controversial observations. Furthermore, the impact of bacterial-derived exogenous H2S on the intestine may be quite different than endogenous cell-based H2S.

Altogether, our results emphasize the importance of the microbial community and its interaction with the host in the pathogenesis of IBD. We propose a central role for mitochondrial dysfunction in colitis development. Although it is difcult to identify the cause underpinning observed downregulation of

8 NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419 ARTICLE

mitochondrial proteins, the loss of butyrate producers in paediatric CD patients will lead to decreased butyrate production which in turn will impair mitochondrial functions likely resulting in ROS production, increased epithelial permeability and translocation of commensal microbes across the epithelial barrier. Intriguingly, antibiotic exposure has been shown to reduce the abundance of butyrate producers3 and also markedly increases the risk of CD in children (odds ratio 2.75, 95% condence interval 1.724.38; ref. 38). We have shown that A. parvulum causes colitis in a susceptible mouse model and that the host microbiota is required for colitis development. Therefore, the increased abundance of A. parvulum in paediatric CD patients with severe inammation suggests a key role for this microbe in driving mucosal inammation and disease exacerbation. Our results point toward a mechanism of CD pathogenesis involving a disruption of the mitochondriamicrobiota relationship leading to dysfunctional mitochondria in the presence of microbes themselves involved in a complex interactive microenvironment at the MLI required for mucosal inammation.

Methods

Ethics statement. The collection of samples from paediatric patients was approved by the Research Ethics Board of the Childrens Hospital of Eastern Ontario (CHEO). The protocol for the use of Il10 / mice was approved by the

Institutional Animal Care and Use Committee of University of North Carolina at Chapel Hill and the University of Florida. Written informed consent was obtained from the participants or parents.

Study design and sampling. We conducted a cross-sectional study of all eligible patients under 18 years of age scheduled to undergo colonoscopy for their initial diagnostic work-up for suspected IBD. Those who were not diagnosed with IBD acted as controls. Exclusion criteria relating to known conditions affecting the intestinal microbiota composition included: (1) a body mass index greater than the 95th percentile for age; (2) diabetes mellitus (insulin and non-insulin dependent);(3) infectious gastroenteritis within the preceding 2 months; and (4) use of any antibiotics or probiotics within 4 weeks before colonoscopy, or (5) irritable bowel syndrome. These same exclusion criteria were applied to the non-IBD control group. All IBD cases met the standard diagnostic criteria following thorough clinical, endoscopic, histologic and radiological evaluation39. Disease location was based on endoscopic and radiologic appearance and characterized using the Paris modication of the Montreal Classication for IBD40. Clinical disease activity of CD was determined using the PCDAI41 and of UC using the PUCAI (Pediatric Ulcerative Colitis Activity Index)42. All the controls had macroscopically and histologically normal mucosa, and did not carry a diagnosis for any known chronic intestinal disorder (for example, celiac disease, eosinophilic enterocolitis, irritable bowel syndrome). Subject clinical data were collected and managed using REDCap (Research Electronic Data Capture) hosted at the CHEO Research Institute. REDCap is a secure, web-based application designed to support data capture for research studies43.

Mucosalluminal interface samples were collected from the right (ascending) colon (RC); n 94, 63 and 37 for CD, control and UC, respectively at the time of

diagnostic colonoscopy. Colonoscopy preparation was done as per standard protocol modied to one day44. During colonoscopy, once the cecum and proximal ascending colon was reached, any loose uid and debris was aspirated. Thereafter, sterile water was ushed onto the mucosa to dislodge the mucus layer from mucosal epithelial cells and the mixture of water, mucus and intestinal cells was aspirated into a sterile container through the colonoscope. These samples were immediately placed on ice in the endoscopy suite and immediately transferred to the lab for processing and storage at 80 C. In addition, mucosal biopsies (n 21

and 10 for CD and control, respectively) were collected from macroscopically involved and non-involved (if the latter were identied) areas of the mid right colon. The biopsies were ash-frozen on dry ice in the endoscopy suite and immediately stored at 80 C until further processing.

Metagenomic DNA extraction. Metagenomic DNA was extracted using the Fast DNA Spin Kit (MP Biomedicals, Solon, OH, USA) and a FastPrep-24 instrument (MP Biomedicals) with two mechanical lysis cycles at speed 6.0 m s 1 for 40 s, with 5 min cooling on ice between the two cycles. Next, the DNA was isolated and puried as per the Fast DNA Spin Kit protocol. The extracted DNA was quantied by Qubit uorometer (Invitrogen, Carlsbad, CA, USA) and stored at 20 C

until use.

16S rDNA-V6 library construction for Illumina sequencing. The V6 region of 16S rDNA was amplied using two successive PCR reactions. The rst PCR added

the Illumina paired-end sequencing adaptors and the barcode sequences using modied universal 16S rDNA-V6 primers45 (Supplementary Table 2, Supplementary Fig. 16). Each reaction was prepared in a total volume of 50 ml using 50 ng of the extracted DNA, 0.5 mM of each primer and 1 Phusion Flash

High-Fidelity PCR Master Mix (Thermo Scientic, Vilnius, Lithuania). The reaction was heated to 98 C for 30 s and then subjected to 10 cycles of 98 C for 5 s, 61 C for 15 s with 1 C drop each cycle and 72 C for 15 s, followed by additional 15 cycles with an annealing temperature of 51 C for 15 s and a nal extension at 72 C for 5 min. The second PCR was carried out using 10 ml of the rst PCR products in a nal volume of 50 ml using the primers PCRFWD1/PCRRVS1 (Supplementary Table 2). The PCRFWD1/PCRRVS1 primers are complementary to the ow cell primers at the 50 end and Illumina paired-end sequencing adaptor at their 30 end. The second PCR conditions were 1 min at 98 C, 15 cycles of 10 s at 98 C, 30 s at 65 C and 30 s at 72 C followed by a nal extension step at 72 C for 5 min. The amplicons of each sample were visualized on 1.5% agarose gel and puried using a Montage PCR96 Cleanup Kit (Millipore, Billerica, MA, USA).

Afterwards, the DNA concentration in each reaction was quantied usingthe Qubit uorometer (Invitrogen, Carlsbad, CA, USA) according to the manufacturers instructions. Finally, equimolar quantities of the amplicon from all the samples were pooled, gel puried using a QIAquick Gel Purication Kit (Qiagen, Hilden, Germany) and sent to The Center for Applied Genomics, Toronto, for 100 bp paired-end Illumina HiSeq2500 sequencing.

Illumina sequencing of pure culture bacterial cultures. Pure cultures of various bacterial species (Actinomyces odontolyticus ATCC 17424, Campylobacter jejuni NCTC 11168, Enterococcus faecalis ATCC 14433, Fusobacterium nucleatum ATCC 25586, Lactobacillus rhamnosus ATCC 7464, Streptococcus cristatus clinical isolate) were grown as required under either anaerobic, microaerophilic (10% O2, 5% CO2, 85% N2) or aerobic conditions at 37 C. DNA for each sample was extracted using the Fast DNA Spin Kit as described above (hereafter referred to as pure samples) and Sanger sequenced to conrm the identity of each bacterial species. In addition, aliquots of the extracted DNA were repuried a second time using the FastDNA kits to assess the impact of potential kit contamination (hereafter referred to as re-extract samples). 16S rDNA-V6 libraries for all the samples were constructed, with triplicate libraries constructed for the re-extract samples and sent for Illumina sequencing as described above.

Illumina microbiota sequencing data analysis. Raw data were rst analysed for decreasing quality along the read. As expected, we noted decreasing quality near the end of the reads. Given that our paired end reads are expected to overlap by B2125 nt (Supplementary Fig. 16), we trimmed four nucleotides from the 30 end of the sequences using the fastx_trimmer script (http://hannonlab.cshl.edu/

Web End =http://hannonlab.cshl.edu/) to remove error-prone bases that could result in incorrect read merging. To note, this trimming step reduced the total number of singletons/doubletons in the nal rareed OTU table by 5% (data not shown). The trimmed paired-end sequences were merged using Fast Length Adjustment of Short reads46. More than 95% of the reads were successfully merged, while the sequences that failed to merge were discarded. The merged reads were then quality ltered with a minimum quality score of 20 over 90% of the sequence using the fastq_quality_lter command from the Fastx toolkit (http://hannonlab.cshl.edu/

Web End =http://hannonlab.cshl.edu/). High-quality reads were demultiplexed according to the forward and the reverse barcode sequences and the barcode trimmed using NovoBarCode (Novocraft.com). Afterwards, the reads were converted from fastq to fasta and were fed to QIIME 1.8.0 (ref. 47) to determine the taxonomic and diversity proles of the samples. First, reads were clustered into OTUs using a closed-reference OTU picking workow against the Greengenes reference set (release 4 February 2011) based on an average percentage of identity of 97%. Note that the same Greengenes database was used for the Illumina, Ion Torrent and 454 pyrosequencing data sets discussed below. Singletons and doubletons were removed and a table of OTU counts per sample was generated. Next, the OTU table was randomly subsampled to a total number of reads per sample of 200,000. The resulting rareed OTU table was used to analyse the microbiota structure and diversity using the microbial ecology tools available in the QIIME package and as input for all other downstream analyses. Any further analysis-specic ltering of the input OTU table is described within each individual Methods section. The phylogenetic tree was constructed using the FastTree method of QIIME 1.8.0 and the generated le was exported to the Interactive Tree of Life (iTOL) software48,49 to view and format the constructed tree.

Illumina pure culture sequencing data analysis. The raw reads obtained for the pure culture sequencing experiments were processed as described above for the Illumina metagenome analysis to generate an OTU table rareed to 200,000 reads per sample and analysed in phyloseq50. Each pure culture consisted of single dominate OTU representing with a median of 499% of the total reads and with a taxonomic classication matching the expected bacteria at the genus level (Supplementary Fig. 5, Supplementary Data 5). Any OTUs present with a taxonomic classication matching the expected bacteria at the family level were considered to be potential sequencing errors. The reads that were matched against these erroneous OTUs were extracted and re-aligned using BLASTN against the dominant OTU to determine the PCR/sequencing error rate (Supplementary Fig. 6,

NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419

Supplementary Data 5). OTUs present in 490% of the samples (in either the pure or re-extract sample sets) were considered to be background contamination originating from either the DNA extraction kits or reagents used during PCR amplication (Supplementary Fig. 5).

Library construction for Ion Torrent sequencing. The V4 or V6 region of 16S rRNA was amplied in a single PCR reaction that incorporated both the Ion Torrent adaptors and an 11 bp barcode sequence using modied universal 16S rDNA-V4 (ref. 51) and V6 (ref. 45) primers (Supplementary Table 3). Note that the conserved V6 sequence used for the Ion Torrent sequencing is identical to those used in the 16S rDNA-V6 Illumina library construction described above. Each reaction was prepared in a total volume of 50 ml using 50 ng extracted metagenomic

DNA, 0.5 mM of each primer, and 1 Phusion Flash High-Fidelity PCR Master

Mix. For the V6 reactions, the PCR cycling conditions were identical to the rst PCR reaction described in the 16S rDNA-V6 Illumina library construction Methods section. For the V4 reactions, the reaction was heated to 98 C for 30 s and then subjected to 34 cycles of 98 C for 30 s, 55 C for 30 s and 72 C for 90 s, followed by a nal extension at 72 C for 5 min as previously described51. The samples were puried using Purelink PCR purication columns and quantied using a Qubit uorometer (Life Technologies). PCR products were visualized on a1.5% agarose gel to ensure successful amplication and that the products were the expected size. Two hundred nanograms of each sample was subsequently pooled together for both the V4 and V6 libraries. The pooled libraries were size-selected using Agencourt AMPure XP DNA purication beads to remove primer dimers. The molar concentrations of each library pool were determined using an Agilent Bioanalyzer and the V4/V6 libraries were combined in an equimolar ratio. This equimolar V4/V6 pool was subjected to emulsion PCR using the Ion OneTouch HiQ Template Kit and sequencing of the templated beads was completed on a single Proton Chip on an Ion Torrent Proton using an Ion Proton HiQ sequencing kit according to the manufacturers instructions (Life Technologies).

Ion Torrent sequencing data analysis. Ion Torrent sequencing was done for both the V4/V6 regions from 36 right colon samples (n 13, 12, 11 for control, CD and

UC). The unaligned BAM le generated by the Ion Torrent Proton (note that adaptor sequences/primer dimers are automatically removed by the Proton analysis pipeline) was downloaded, converted to fastq format and reads split according to the expected sizes for the V4/V6 regions plus the barcodes (234266 and 162194, respectively, representing a 16 bp window), using samtools52 in conjunction with bioawk (https://github.com/lh3/bioawk

Web End =https://github.com/lh3/bioawk). The reads were quality trimmed using the modied Mott trimming algorithm implemented in seqtk (https://github.com/lh3/seqtk

Web End =https://github.com/lh3/ https://github.com/lh3/seqtk

Web End =seqtk ) and reads that were shorter than the expected size for either the V4 or V6 regions (as described above) after trimming were discarded. The V4/V6 reads were then demultiplexed using FastqMultX53, converted from fastq to fasta and fed to QIIME 1.8.047 to determine the taxonomic proles of the samples. The reads were clustered into OTUs using a closed-reference OTU picking workow against the Greengenes reference set (release 4 February 2011) based on an average percentage of identity of 97%. The singletons/doubletons were removed from both the V4/V6 OTU tables and the number of reads/sample was normalized across sample groups (100,000 reads per sample for V4 and 200,000 reads per sample for V6). The rareed tables were used to calculate alpha/beta diversities in phyloseq50 and to identify the microbial biomarkers for each tested category using the linear discriminant analysis effect size algorithm11 with the default parameters. To compare the Ion Torrent sequencing results to the Illumina results, the Illumina V6 data from the samples subjected to Ion Torrent sequencing was extracted from the Illumina rareed OTU table (note that the Illumina V6 and Ion Torrent V6 samples were rareed to the same depth) and analysed as described above for the Ion Torrent data.

16S rDNA-V6 454-pyrosequencing. The 454-pyrosequencing sequencing library of the 16S rRNA-V6 region was constructed using a two-step PCR strategy using the primers; 16SF 50-AAACTCAAAKGAATTGACGG-30 and 16SR 50-ACGAGCTGACGACARCCATG-30 (ref. 45). Ten replicates from each sample were amplied in a 25 ml reaction containing 50 ng DNA template, 1 PCR buffer,

2 mM MgCl2, 0.2 mM dNTPs mix, 0.2 mM of each primer and 2.5 U Platinum Taq polymerase (Invitrogen, Carlsbad, CA, USA). The reactions were heated to 95 C for 5 min, followed by 15 cycles of 94 C for 40 s, 46 C for 1 min and 72 C for 30 s, and a nal extension at 72 C for 5 min. The amplicons from each sample were pooled and puried with Qiagens MiniElute PCR purication columns (Qiagen, Hilden, Germany) following its standard protocol and eluted in 25 ml molecular biology-grade water. Two microlitres of the pooled puried amplicons were then used as a template for the second PCR following the same PCR conditions except using 30 amplication cycles instead of 15 cycles. The second PCR added the titanium tails to the 16S-V6 amplicons, which are essential for the 454 sequencing procedure. The amplicons from the second PCR were pooled and puried again using Qiagen MiniElute PCR purication columns. Negative control reactions (no DNA template) were included in all the experiments. PCR products were visualized on a 1.5% agarose gel to check the amplication success. The puried amplicons from different samples were normalized to the same concentration (100 ng ml 1), captured to streptavidin coated sepharose beads and exposed to emulsion PCR.

Afterwards, the immobilized double stranded amplicons were denatured into single-stranded DNA, which was then annealed to the sequencing primer by heating at 65 C for 5 min. The single-stranded DNA-containing beads were sequenced on a 454 Genome Sequencer FLX System (Roche Diagnostics GmbH) using GS Titanium chemistry according to the standard amplicon sequencing protocol. Samples covering the three tested phenotypes (two to four samples each) were sequenced on the same run with each being sequenced in a 1/16 section of 70 75 picotiter plate. Note that as the samples were physically separated from

each other during the sequencing reaction, there was no need to use barcoded primers.

454-pyrosequencing data analysis. A total of 346,160 reads were generated from the 454 pyrosequencing of 16S rDNA-V6 region from 26 right colon samples. The raw pyrosequencing reads were de-noised and processed to remove low quality and short reads using Quantitative Insights Into Microbial Ecology pipeline release1.8.0 (QIIME 1.8.0; ref. 47) according to the following parameters: (1) minimum read length of 100 bp, (2) exact matching to the sequencing primers, (3) no ambiguous nucleotides and (4) a minimum average quality score of 20. Next, sequences were de novo clustered using UCLUST54 based on average percentage of identity of 97%. The most abundant read from each OTU was picked as a representative sequence for that cluster. The representative sequences were then aligned using PyNAST with a minimum alignment length of 100 and a minimum percentage identity of 75%, followed by checking the chimeric OTUs with the blast_fragments approach. Only six representative sequences were identied as chimeric and therefore were removed. Taxonomy assignments were made with BLAST55 by searching the representative sequences against the Greengenes database (release 4 February 2011) with an e value o1e 8 and a condence score

of Z0.5. Next, singletons (OTUs that had only one matching sequence) were ltered out from the resulting OTU table. The OTU table was then used to determine the alpha and beta diversity within and between the samples using the QIIMES default criteria. To identify the microbial biomarkers for each tested category, the relative abundance of different phylogenetic levels computed by QIIME was analysed by linear discriminant analysis effect size algorithm11 using its default parameters. The 454-pyrosequencing reads assigned as Atopobium by QIIME analysis were retrieved and found to match to A. parvulum following alignment of the reads against the RDP56 and NCBI databases (the aligned region covered the entire 454 sequence length with 499% sequence identity to A.

parvulum).

Statistical analysis of the microbiota data. Several statistical approaches were used to identify taxa signicantly associated with disease status and severity. Identication of OTUs exhibiting differential abundance between patients with different disease severity (mild versus moderate versus severe) was performed using the metagenomeSeq R packagea statistical approach that controls for confounding factors and accounts for undersampling12. Briey, the OTU counts were normalized using a cumulative sum scaling approach and OTUs with a normalized count o5 in at least one of the groups were discarded. OTUs exhibiting differential abundance as a function of disease severity were identied using a zero-inated

Gaussian mixture model by incorporating gender and the inammation status of the ascending colon (inamed or not inamed) as confounding covariates. Differentially abundant OTUs were selected according to the following algorithm:(1) OTU present in more than 50% of the samples from at least one group;(2) a fold change in relative abundance Z2 and (3) a statistical signicance with an adjusted P value (calculated using the BenjaminiHochberg method) o0.05.

A KruskalWallis test with a post hoc Dunns test was performed to compare the relative abundance of specic taxa or groups of bacteria as a function of disease status (CD versus UC versus control) and disease severity (mild versus moderate versus severe). The relative abundance of the taxa identied was also analysed by principal coordinate analysis (PCoA). All the statistical analyses were performed using XLSTAT (Addinsoft, NY, USA) and/or GraphPad Prism version 6 (GraphPad, La Jolla, CA, USA).

Validation of V6-16S sequencing. The quantication of A. parvulum and sulfate reducing bacteria relative to total bacteria was determined by conducting quantitative PCR (qPCR) on the extracted metagenomic DNA using an Applied Bio-systems 7,300 and primers listed in Supplementary Table 4. Each sample was tested in duplicate in a total volume of 25 ml per reaction. Fifty nanograms of template

DNA was added to a reaction mixture containing 0.3 mM of each primer (0.5 mM for 16S rDNA universal primers) and 1 QuantiFast SYBR Green PCR master

mix (Qiagen, Hilden, Germany). The amplication conditions were 5 min at 95 C followed by 40 cycles of 95 C for 10 s and 66 C for 1 min with data collection at the second step of each cycle. Ct values were then extracted using the Applied Biosystems 7300 software version 1.3.1 and the relative abundance of each taxon was calculated as a DCt value (taxon Ctuniversal 16S rRNA Ct). To validate the specicity of Apar-711F and Aparv-881R, fresh PCR amplicons from total DNA extracted from two different mucosal aspirates was cloned using the TOPO TA cloning vector (Invitrogen, Carlsbad, CA, USA) according to the manufacturers instructions. Next, the plasmid containing the 16S rDNA gene fragment was extracted from six different clones by the QIAprep Spin Miniprep kit (Qiagen,

10 NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419 ARTICLE

Hilden, Germany) followed by capillary sequencing using M13F and M13R primers.

Quantication of butyrate-producers by qPCR. The overall abundance of butyrate-producing bacteria was determined by quantifying the amount of the butyryl CoA:acetate CoA-transferase (BCoAT) gene using the primers BCoATscrF/ BCoATscrR (Supplementary Table 4) as described elsewhere57,58. The BCoAT gene was amplied from 50 ng of metagenomic DNA in a 25 ml qPCR reaction containing 1 QuantiTect SYBR Green PCR master mix (Qiagen) and 2.5 mM of

BCoATscrF/BCoATscrR primers. The amplication conditions were as follows: one cycle of 95 C for 15 min; 40 cycles of 94 C for 15 s, 53 C and 72 C each for 30 s with data acquisition at 72 C. The 16S rDNA gene was used to normalize total bacterial content in each sample and was amplied using adapted universal primers UniF/UniR (Supplementary Table 4)59. The assay was done in duplicate for each sample. Delta-Ct (DCt) for each target was calculated by subtracting the Ct of the 16S rDNA from the BCoAT Ct. Then, the DCt values were compared between the CD and control groups using a MannWhitney two-tailed test with a Dunns post hoc test and P values below 0.05 were considered signicant.

Microbial correlation network analysis. A correlation network was constructed to identify co-occurrence and mutual exclusion among OTUs from the control subjects and CD patients. To focus on bacteria potentially associated with disease, only the OTUs exhibiting differential abundance as a function of disease severity were chosen for this network analysis. In addition, OTUs that were found to be absent in more than two-third of the samples or with a minimum occurrence o15 across all the samples were removed from the data set, leaving a matrix of 97 OTUs. The network was built as previously described60 using the CoNet Cytoscape plug-in. The correlation scores were calculated for each OTU pair using a combination of two correlation (Pearson and Spearman) and two dissimilarity measures (BrayCurtis and KullbackLeibler). The 400 top- and bottom-ranking edges from each method were retrieved. Edge- and method-specic permutations (100) were used to control for potential false-positive correlations. The resulting distribution was used to calculate statistical signicance for each edge using 100 bootstrap iterations. The P value was adjusted for multiple comparisons using a BenjaminiHochberg correction and a value o0.05 was deemed signicant. The

P value was obtained for each correlation or dissimilarity measure and only edges supported by at least two signicant measures were included. The nal network consisted of 89 nodes (OTUs) and 341 edges (correlations). The resulting network was visualized using Cytoscape 3.2.1 (ref. 61).

Stable isotope labelling by amino acids in cell culture. The human hepatic HuH7 cells (HuH-7) and human embryonic kidney 293 cells (HEK-293) were originally acquired from ATCC and the human colorectal cancer 116 cells (HCT-116) were obtained from the Japanese Collection of Research Bioresources Cell Bank. HuH-7, HEK-293 and HCT-116 were individually grown at 37 C in a 5% CO2 humidied incubator (to note, all three cell lines used in this study were found to be mycoplasma positive). Stable isotope labelling by amino acids in cell culture (SILAC) medium was prepared as follows: DMEM lacking lysine, arginine and methionine was custom prepared by AthenaES (Baltimore, MD, USA) and supplemented with 30 mg l 1 methionine (Sigma-Aldrich; Oakville, ON, Canada), 10% (v/v) dialysed FBS (GIBCO-Invitrogen; Burlington, ON,Canada), 1 mM sodium pyruvate (Gibco-Invitrogen), 28 g ml 1 gentamicin (Gibco-Invitrogen), and [13C6,15N2]-L-lysine, [13C6,15N4]-L-arginine (heavy isotopic form of amino acids; Heavy Media) from Sigma-Aldrich (Oakville, ON, Canada) at nal concentrations of 42 and 146 mg l 1 for arginine and lysine, respectively.

For HCT-116, the concentration of arginine was increased to 84 mg l 1. The cells were grown for at least 10 doublings in SILAC media to allow for complete incorporation of the isotopically labelled amino acids into the cells.

Determination of the rate of SILAC amino acids incorporation. The cells were grown to 80% conuency in SILAC medium (5 106 cells were plated in 10-cm

dish). Next, the cells were washed twice with ice-cold phosphate-buffered saline and lysed by addition of 1 ml of 1X RIPA buffer (50 mM Tris (pH 7.6), 150 mM NaCl, 1% (v/v) NP-40, 0.5% (w/v) deoxycholate, 0.1% (w/v) SDS with protease inhibitor cocktail (Complete Mini Roche; Mississauga, ON,Canada) and phosphatase inhibitor (PhosStop Roche tablet)). The lysates were then transferred to 15 ml conical tubes and the proteins were precipitated by addition of 5 ml ice-cold acetone followed by incubation at 20 C overnight. The proteins were collected

by centrifugation (3,000g, 10 min, 4 C), washed with ice-cold acetone two times and the protein pellets were re-solubilized in 300 ml of a 50 mM NH4HCO3 solution containing 8 M urea. The protein concentrations were determined by the Bradford dye-binding method using Bio-Rads Protein Assay Kit (Mississauga, ON, Canada). For the general in-solution digestion, 200 mg of protein lysates were reconstituted in 50 mM NH4HCO3 (200 ml) and the proteins were reduced by mixing with 5 ml of 400 mM DTT at 56 C for 15 min. The proteins were then subjected to alkylation by mixing with 20 ml of 400 mM iodoacetamide in darkness (15 min at room temperature) followed by addition of 800 ml of 50 mM NH4HCO3 to reduce the urea concentration to B0.8 M. Next, the proteins were digested with TPCK-trypsin solution (nal ratio of 1:20 (w/w, trypsin:protein) at 37 C for 18 h. Finally, the

digested peptides were desalted using C18 Sep-Pack cartridges (Waters), dried down in a speed-vac, and reconstituted in 0.5% formic acid before mass spectro-metric analysis (as described below) and the determination of labelling efciency. The incorporation efciency was calculated according to the following equation: (1-1/Ratio(H/L)); where H and L represents the intensity of heavy and light peptides detected by mass-spectrometry, respectively. Labelling was considered complete when values reached at least 95% for each cell type.

Proteomic analysis of biopsies. Biopsies were lysed in 4% SDS, 50 mM Tris-HCl (pH 8.0) supplemented with proteinase inhibitor cocktail (Roche) and homogenized with a Pellet pestle. The lysates were sonicated three times with 10 s pulses each with at least 30 s on ice between each pulse. The protein concentrations were determined using the Bio-Rad DC Protein Assay. The proteins were processed using the Filter Aided Sample Preparation method as previously described with some modications62. Biopsy lysates (45 mg of proteins) and heavy SILAC-labelled cell lysates (15 mg from each HuH-7, HEK- 293 and HCT-116 cells) were mixed at a 1:1 weight ratio and transferred into the lter. The samples were centrifuged (16,000g, 10 min), followed by two washes of 200 ml 8 M urea, 50 mM Tris-HCl pH8.0. The samples were then reduced by incubation in 200 ml of 8 M urea, 50 mM Tris-HCl (pH 8.0) supplemented with 20 mM dithiothreitol. After centrifugation, the samples were subjected to alkylation by adding 200 ml of 8 M urea, 50 mM

Tris-HCl pH 8.0, containing 20 mM iodoacetamide (30 min at room temperature protected from light). The samples were washed twice using 200 ml 8 M urea,50 mM Tris-HCl pH 8.0 to remove excess SDS. To further dilute urea, two washes of 200 ml 50 mM Tris-HCl pH 8.0 were performed. For the trypsin digest, the samples were incubated in 200 ml of 50 mM Tris-HCl pH 8.0, containing 5 mg of Trypsin (TPCK Treated, Worthington) on a shaker (250 rpm) at 37 C overnight. Finally, 200 ml of 50 mM Tris-HCl pH 8.0 was added to elute the peptides by centrifugation (twice). The peptides were fractionated, using an in-house constructed SCX column with ve pH fractions (pH 4.0, 6.0, 8.0, 10.0 and 12.0). The buffer composition was 20 mM boric acid, 20 mM phosphoric acid and 20 mM acetic acid, with the pH adjusted by using 1 M NaOH. Finally, the fractionated samples were desalted using in-house C18 desalting cartridges and dried in a speed-vac before LC-MS analysis.

Mass-spectrometry analyses. All resulting peptide mixtures were analysed by high-performance liquid chromatography/electrospray ionization tandem mass spectrometry (HPLC-ESI-MS/MS). The HPLC-ESI-MS/MS consisted of an automated ekspert nanoLC 400 system (Eksigent, Dublin, CA, USA) coupled with an LTQ Velos Pro Orbitrap Elite mass spectrometer (ThermoFisher Scientic,San Jose, CA, USA) equipped with a nanoelectrospray interface operated in positive ion mode. Briey, each peptide mixture was reconstituted in 20 ml of 0.5%

(v/v) formic acid and 12 ml was loaded on a 200 mm 50 mm fritted fused silica

pre-column packed in-house with reverse phase Magic C18AQ resins (5 mm; 200 pore size; Dr Maisch GmbH, Ammerbuch, Germany). The separation of peptides was performed on an analytical column (75 mm 10 cm) packed with reverse phase

beads (3 mm; 120 pore size; Dr Maisch GmbH, Ammerbuch, Germany) using a 120 min gradient of 530% acetonitrile (v/v) containing 0.1% formic acid (v/v) (JT Baker, Phillipsburg NJ, USA) at an eluent ow rate of 300 nl min 1. The spray voltage was set to 2.2 kV and the temperature of heated capillary was 300 C. The instrument method consisted of one full MS scan from 400 to 2,000 m/z followed by data-dependent MS/MS scan of the 20 most intense ions, a dynamic exclusion repeat count of 2 and a repeat duration of 90 s. The full mass was scanned in an Orbitrap analyzer with R 60,000 (dened at m/z 400), and the subsequent

MS/MS analyses were performed in LTQ analyzer. To improve the mass accuracy, all the measurements in the Orbitrap mass analyzer were performed with on-they internal recalibration (Lock Mass). The charge state rejection function was enabled with charge states unassigned and single states rejected. All the data were recorded with Xcalibur software (ThermoFisher Scientic, San Jose, CA, USA).

Mass-spectrometry database search and bioinformatic analysis. All the raw les were processed and analysed by MaxQuant, Version 1.2.2.5 against the decoy Uniprot-human database (86,749 entries), including commonly observed contaminants. The following parameters were used: cysteine carbamidomethylation was selected as a xed modication, with methionine oxidation, protein amino(N)-terminal acetylation and heavy proline set as variable modications. Enzyme specicity was set to trypsin. Up to two missing cleavages of trypsin were allowed. SILAC double labelling (light: K0R0; heavy: K8R10) was set as the search parameter to assess the conversion efciency. The precursor ion mass tolerances were 7 p.p.m. and the fragment ion mass tolerance was 0.5 Da for MS/MS spectra. The false discovery rate for peptides and proteins was set at 1% and a minimum length of six amino acids was used for peptide identication. The peptides le was imported into Perseus (version 1.2.0.17) to extract the lysine- and arginine-containing peptides, respectively. The protein-group le was imported into Perseus (version 1.3.0.4) for statistical analysis. Differences between control and CD protein expression were determined using a t-test with an adjustment for multiple comparisons (BenjaminiHochberg method; false discovery rate 0.05) and q values

of less than 0.05 were considered signicant. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was achieved using the DAVID Bioinformatics

NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419

Resources (http://david.abcc.ncifcrf.gov/

Web End =http://david.abcc.ncifcrf.gov/). DAVID statistical analyses were performed against the whole genome. Proteomics has a tendency to oversample proteins from the cytosol and nucleus while under-sampling membrane-associated proteins. Therefore, the results from DAVID were re-tested against the set of proteins that were not changing in our data set to eliminate any pathway/GO enrichment biases.

Total RNA extraction from intestinal mucosal aspirates. RNA integrity was preserved by adding an equal volume of RNAlater (Ambion) to the mucosal aspirates before freezing at 80 C. The frozen aliquot (2 ml) was thawed on ice

and the total RNA was extracted following the hot phenol protocol described previously63. Briey, 4 ml of each sample in RNAlater were pelleted by centrifugation at 13,000 g for 5 min at 4 C. The pellets were washed twice by

resuspension in 50% RNAlater/PBS buffer and centrifugation at 13,000g for 5 min at 4 C. Next, the pellets were resuspended and lysed in 2 ml denaturing buffer (4 M guanidium thiocyanate, 25 mM sodium citrate, 0.5% N-laurylsarcosine and 1% N-acetyl cysteine, 0.1 M 2-mercaptoethanol). The lysate was divided into 500 ml aliquots, to which 4 ml of 1 M sodium acetate (pH 5.2) was added. Each aliquot was then incubated with 500 ml of buffer saturated phenol (pH 4.3) at 64 C for 10 min with intermittent mixing. The tube was placed on ice for 10 min followed by centrifugation at maximum speed for 30 min at 4 C. One millilitre of chloroform was added to the aqueous layer and incubated for 15 min on ice, followed by centrifugation at 18,000g for 30 min at 4 C. Afterwards, RNA was precipitated from the aqueous layer by the addition of one-tenth volume 3 M sodium acetate, 500 mM DEPC-treated EDTA and two volumes of cold ethanol followed by overnight incubation at 80 C. Later, the RNA was pelleted by centrifugation at

maximum speed for 30 min at 4 C, aspirated with 80% cold ethanol and resuspended in 100 ml of RNAse free dH2O. The extracted RNA was treated twice with Dnase I (Epicentre) followed by PCR amplication using the 16S rDNA universal primers; Bact-8F and 1391-R (ref. 64), to conrm the absence of genomic DNA. The quality and the quantity of the extracted RNA were determined by NanoDrop 2000 spectrophotometer (ThermoScientic) and conrmed by Bio-Rads Experion StdSens RNA system according to the manufacturers instructions and stored at 80 C until use.

Quantication of H2S detoxication genes expression level. The quantication of the expression level of TST, SQRDL (Sulde Quinone Reductase Like) and COX4-1 (Cytochrome C oxidase subunit IV isoform 1) relative to hGAPDH (Glyceraldehyde-3-Phosphate Dehydrogenase) was determined using an Applied Biosystems 7300 and Quantitect SYBR Green RT-PCR kit (Qiagen, Hilden, Germany). The primers used were either designed using the NCBI Primer-BLAST tool65 or extracted from the literature (Supplementary Table 5). The specicity of the primers was conrmed by TOPO TA cloning and capillary DNA sequencing as previously mentioned in the qPCR validation of V6-16S sequencing methods section. Each reaction contained 100 ng RNA template, 0.5 mM of each primer, 1

Quantitect SYBR Green RT-PCR master mix and 0.25 ml Quantitect RT-mix. The one-step qRTPCR conditions were 50 C for 30 min, 95 C for 15 min followed by 40 cycles of 15 s at 94 C, 30 s at 60 C and 30 s at 72 C with data collection at the third step of each cycle. The amplication specicity was checked by the melting prole of the amplicon and 2% agarose gel electrophoresis. The Ct values were then extracted using the Applied Biosystems 7300 software version 1.3.1. The Ct values of TST, SQRDL or COX4 were normalized to the Ct values of hGAPDH generating DCt. Next, DDCt was calculated by subtracting the average DCt of the control group from the DCt of each sample. The relative quantication was then calculated as 2 DDCt as described previously66.

Microbiotamitochondrial proteins correlation analysis. Pairwise correlations between the relative abundance of the signicant OTUs and mitochondrial proteins identied via MS were calculated using the non-parametric Kendall Tau and Spearmans rank correlation coefcients. To avoid spurious correlations, OTUs and mitochondrial proteins absent in more than half of the samples in each group were excluded from each pairwise comparison. As a result, pairwise correlations were calculated between 62 OTUs and 106 mitochondrial proteins from the 21 patients for which we had both data sets. To further avoid spurious correlations, correlations were considered signicant only if the upper-tail probability occurrence was below 0.05 for both the Kendall Tau and the Spearmans rank statistics.

Il10 / mouse experiments and tissue processing. Germ-free SvEv129/ C57BL6 Il10 / ;NF-kBEGFP male/female mice (812 weeks old, n 12) were

transferred to specic pathogen free (SPF) conditions and after 2 weeks mice from a randomly selected cohort (n 6) were gavaged once weekly with A. parvulum

(1 108 CFUs) for 6 weeks and sacriced 6 weeks later (14 weeks total).

A. parvulum ATCC 33793 was grown in fastidious anaerobic broth (Lab M, Canada). To investigate the involvement of complex microbiota in the development of colitis, we performed two subsequent experiments using C57BL6 Il10 /

male/female mice. In the rst experimental setting, gnotobiotic Il10 / mice (n 16) were randomized into two groups; 1: GF only (n 6) and 2: GF A.

parvulum (1 108 CFUs; n 10). The mice were killed after 8 weeks of mono-

association. The mice were fed normal diet (Teklan Global 18% Protein Rodent

Diet). In the second experimental setting, gnotobiotic Il10 / male/female mice (n 31) were transferred to SPF conditions for 2 weeks and then randomized into

two groups; 1: SPF only (n 7), 2: SPF plus A. parvulum (1 108 CFUs; n 8), 3:

SPF plus bismuth (III) diet (n 8) and 4: SPF plus A. parvulum and bismuth (III)

diet (n 8). As described for the rst cohort, the mice were gavaged weekly with

A. parvulum (1 108 CFUs) for 6 weeks and then killed after 6 weeks of coloni

zation (14 weeks total). Bismuth (III) subsalicylate (Sigma-Aldrich, Saint Louis, MO, USA) was incorporated to the chow (Teklan Global 18% Protein Rodent Diet) at a concentration of 7 g kg 1 (Harlan Laboratories, Madison, WI, USA) and then irradiated for gnotobiotic experiments. The mice were fed with this diet 1 week before the colonization with A. parvulum. The tissue samples from the colon were collected for RNA and histology as previously described67. All the animal procedures were approved by the University of North Carolina at Chapel Hill Animal Care and Use Committee. The histological images were acquired using a DFC310 FX (LEICA) coupled to LEICA Application Suite AFv4.6 software. Intestinal inammation was blindly scored as previously described18. The tissue was divided into four quarters to allow individual scoring of the mid-section, which was affected by A. parvulum more so than in other colitis models. A score was given to each quarter separately and then added to generate a nal colitis score on a scale of 016.

Statistical analyses of Il10 / mice results. Unless specically noted, the statistical analyses were performed using GraphPad Prism version 6 (GraphPad,

La Jolla, CA, USA). The comparisons of mouse studies were made with a nonparametric analysis of variance, and then a MannWhitney U-test. The experiments were considered statistically signicant if Po0.05.

Mouse endoscopy. Colonoscopy was performed using a Coloview System (Karl Storz Veterinary Endoscopy) as described previously68. The mice were anaesthetized using 1.5 to 2% isourane and B4 cm of the colon from the anal verge from the splenic exure was visualized. The procedures were digitally recorded on an AIDA Compaq PC.

qRTPCR on mouse intestinal samples. Total RNA from intestinal tissues was extracted from the distal part of the colon using TRIzol (Invitrogen, Carlsbad, CA, USA) following the manufacturers protocol. The complementary DNA was reverse-transcribed using M-MLV (Invitrogen, Carlsbad, CA, USA) and mRNA expression levels were measured using SYBR Green PCR Master mix (Applied Biosystems) on an ABI 7900HT Fast Real-Time PCR System and normalized to b-actin. The primers used are listed in Supplementary Table 5. The PCR reactions were performed for 40 cycles according to the manufacturers recommendations, and RNA fold changes were calculated using the DDCt method66.

Data availability. The demultiplexed reads for all V6-16S experiments that support the ndings of this study have been deposited to the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra

Web End =http://www.ncbi.nlm.nih.gov/sra). The pure culture experiments were submitted under accession number SRP075131, the Illumina microbiota reads under accession numbers SRP034595/SRP056939, the 454 microbiota reads under accession number SRP034632 and the Ion Torrent reads under accession number SRP075131. All the raw data les for the mass spectrometry that support the ndings of this study have been submitted to ProteomeXchange under accession number PXD002882. The remaining data are available within the Article and Supplementary Information les or available from the author on request.

References

1. Rigoli, L. & Caruso, R. A. Inammatory bowel disease in pediatric and adolescent patients: a biomolecular and histopathological review. World J. Gastroenterol. 20, 1026210278 (2014).

2. Li, J., Butcher, J., Mack, D. & Stintzi, A. Functional impacts of the intestinal microbiome in the pathogenesis of inammatory bowel disease. Inamm. Bowel Dis. 21, 139153 (2015).

3. Gevers, D. et al. The treatment-naive microbiome in new-onset Crohns disease. Cell Host Microbe 15, 382392 (2014).

4. Hansen, R. et al. Microbiota of de novo pediatric IBD: increased Faecalibacterium prausnitzii and reduced bacterial diversity in Crohns but not in ulcerative colitis. Am. J. Gastroenterol. 107, 19131922 (2012).

5. Kolho, K. L. et al. Fecal microbiota in pediatric inammatory bowel disease and its relation to inammation. Am. J. Gastroenterol. 110, 921930 (2015).

6. Morgan, X. C. et al. Dysfunction of the intestinal microbiome in inammatory bowel disease and treatment. Genome Biol. 13, R79 (2012).

7. Papa, E. et al. Non-invasive mapping of the gastrointestinal microbiota identies children with inammatory bowel disease. PLoS ONE 7, e39242 (2012).

8. Presley, L. L. et al. Host-microbe relationships in inammatory bowel disease detected by bacterial and metaproteomic analysis of the mucosal-luminal interface. Inamm. Bowel Dis. 18, 409417 (2012).

12 NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419 ARTICLE

9. McHardy, I. H. et al. Integrative analysis of the microbiome and metabolome of the human intestinal mucosal surface reveals exquisite inter-relationships. Microbiome 1, 17 (2013).

10. Walters, W. A., Xu, Z. & Knight, R. Meta-analyses of human gut microbes associated with obesity and IBD. FEBS Lett. 588, 42234233 (2014).

11. Segata, N. et al. Metagenomic biomarker discovery and explanation. Genome Biol. 12, R60 (2011).

12. Paulson, J. N., Stine, O. C., Bravo, H. C. & Pop, M. Differential abundance analysis for microbial marker-gene surveys. Nat. Methods 10, 12001202 (2013).

13. Salter, S. J. et al. Reagent and laboratory contamination can critically impact sequence-based microbiome analyses. BMC Biol. 12, 87 (2014).

14. Schloss, P. D., Gevers, D. & Westcott, S. L. Reducing the effects of PCR amplication and sequencing artifacts on 16S rRNA-based studies. PLoS ONE 6, e27310 (2011).

15. Carbonero, F., Beneel, A. C., Alizadeh-Ghamsari, A. H. & Gaskins, H. R. Microbial pathways in colonic sulfur metabolism and links with health and disease. Front. Physiol. 3, 448 (2012).

16. Tyanova, S. et al. Proteomic maps of breast cancer subtypes. Nat. Commun. 7, 10259 (2016).

17. Blachier, F. et al. Luminal sulde and large intestine mucosa: friend or foe? Amino Acids 39, 335347 (2010).

18. Karrasch, T., Kim, J. S., Muhlbauer, M., Magness, S. T. & Jobin, C. Gnotobiotic IL-10 / ;NF-kappa B(EGFP) mice reveal the critical role of TLR/NF-kappa B signaling in commensal bacteria-induced colitis. J. Immunol. 178, 65226532 (2007).

19. Kim, S. C., Tonkonogy, S. L., Karrasch, T., Jobin, C. & Sartor, R. B. Dual-association of gnotobiotic IL-10 / mice with 2 nonpathogenic commensal bacteria induces aggressive pancolitis. Inamm. Bowel Dis. 13, 14571466 (2007).

20. Suarez, F. L., Furne, J. K., Springeld, J. & Levitt, M. D. Bismuth subsalicylate markedly decreases hydrogen sulde release in the human colon. Gastroenterology 114, 923929 (1998).

21. Sifroni, K. G. et al. Mitochondrial respiratory chain in the colonic mucosal of patients with ulcerative colitis. Mol. Cell. Biochem. 342, 111115 (2010).

22. OMorain, C., Smethurst, P., Levi, J. & Peters, T. J. Subcellular fractionation of rectal biopsy homogenates from patients with inammatory bowel disease. Scand. J. Gastroenterol. 20, 209214 (1985).

23. Soderholm, J. D. et al. Augmented increase in tight junction permeability by luminal stimuli in the non-inamed ileum of Crohns disease. Gut 50, 307313 (2002).

24. Beltran, B. et al. Mitochondrial dysfunction, persistent oxidative damage, and catalase inhibition in immune cells of naive and treated Crohns disease. Inamm. Bowel Dis. 16, 7686 (2010).

25. Restivo, N. L., Srivastava, M. D., Schafer, I. A. & Hoppel, C. L. Mitochondrial dysfunction in a patient with crohn disease: possible role in pathogenesis.J. Pediatr. Gastroenterol. Nutr. 38, 534538 (2004).26. Bar, F. et al. Mitochondrial gene polymorphisms that protect mice from colitis. Gastroenterology 145, 10551063 (2013).

27. Wang, A. et al. Targeting mitochondria-derived reactive oxygen species to reduce epithelial barrier dysfunction and colitis. Am. J. Pathol. 184, 25162527 (2014).

28. Soderholm, J. D. et al. Epithelial permeability to proteins in the noninamed ileum of Crohns disease? Gastroenterology 117, 6572 (1999).

29. Zeissig, S. et al. Changes in expression and distribution of claudin 2, 5 and 8 lead to discontinuous tight junctions and barrier dysfunction in active Crohns disease. Gut 56, 6172 (2007).

30. Ramasamy, S., Singh, S., Taniere, P., Langman, M. J. & Eggo, M. C. Suldedetoxifying enzymes in the human colon are decreased in cancer and upregulated in differentiation. Am. J. Physiol. Gastrointest. Liver Physiol. 291, G288G296 (2006).

31. Linden, D. R. Hydrogen sulde signaling in the gastrointestinal tract. Antioxid. Redox Signal. 20, 818830 (2014).

32. Quince, C. et al. Extensive modulation of the fecal metagenome in children with Crohns disease during exclusive enteral nutrition. Am. J. Gastroenterol. 110, 17181729 (2015).

33. Miller, T. W. et al. Hydrogen sulde is an endogenous potentiator of T cell activation. J. Biol. Chem. 287, 42114221 (2012).

34. Ijssennagger, N. et al. Gut microbiota facilitates dietary heme-induced epithelial hyperproliferation by opening the mucus barrier in colon. Proc. Natl Acad. Sci. USA 112, 1003810043 (2015).

35. Guo, F. F., Yu, T. C., Hong, J. & Fang, J. Y. Emerging roles of hydrogen sulde in inammatory and neoplastic colonic diseases. Front. Physiol. 7, 156 (2016).

36. Wang, R. Physiological implications of hydrogen sulde: a whiff exploration that blossomed. Physiol. Rev. 92, 791896 (2012).

37. Gerasimidis, K. et al. Decline in presumptively protective gut bacterial species and metabolites are paradoxically associated with disease improvement in pediatric Crohns disease during enteral nutrition. Inamm. Bowel Dis. 20, 861871 (2014).

38. Ungaro, R. et al. Antibiotics associated with increased risk of new-onset Crohns disease but not ulcerative colitis: a meta-analysis. Am. J. Gastroenterol. 109, 17281738 (2014).

39. North American Society for Pediatric Gastroenterology, Hepatology and Nutrition et al. Differentiating ulcerative colitis from Crohn disease in children and young adults: report of a working group of the North American Society for Pediatric Gastroenterology, Hepatology, and Nutrition and the Crohns and Colitis Foundation of America. J. Pediatr. Gastroenterol. Nutr. 44, 653674 (2007).

40. Levine, A. et al. Pediatric modication of the Montreal classication for inammatory bowel disease: the Paris classication. Inamm. Bowel Dis. 17, 13141321 (2011).

41. Hyams, J. et al. Evaluation of the pediatric crohn disease activity index: a prospective multicenter experience. J. Pediatr. Gastroenterol. Nutr. 41, 416421 (2005).

42. Turner, D. et al. Development, validation, and evaluation of a pediatric ulcerative colitis activity index: a prospective multicenter study. Gastroenterology 133, 423432 (2007).

43. Harris, P. A. et al. Research electronic data capture (REDCap)a metadata-driven methodology and workow process for providing translational research informatics support. J. Biomed. Inform. 42, 377381 (2009).

44. Jimenez-Rivera, C., Haas, D., Boland, M., Barkey, J. L. & Mack, D. R. Comparison of two common outpatient preparations for colonoscopy in children and youth. Gastroenterol. Res. Pract. 2009, 518932 (2009).

45. Sundquist, A. et al. Bacterial ora-typing with targeted, chip-based Pyrosequencing. BMC Microbiol. 7, 108 (2007).

46. Magoc, T. & Salzberg, S. L. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27, 29572963 (2011).

47. Caporaso, J. G. et al. QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335336 (2010).

48. Letunic, I. & Bork, P. Interactive Tree Of Life (iTOL): an online tool for phylogenetic tree display and annotation. Bioinformatics 23, 127128 (2007).

49. Letunic, I. & Bork, P. Interactive Tree Of Life v2: online annotation and display of phylogenetic trees made easy. Nucleic Acids Res. 39, W475W478 (2011).

50. McMurdie, P. J. & Holmes, S. Phyloseq: a bioconductor package for handling and analysis of high-throughput phylogenetic sequence data. Pac. Symp. Biocomput. 235246 (2012).

51. Milani, C. et al. Assessing the fecal microbiota: an optimized ion torrent 16S rRNA gene-based analysis protocol. PLoS ONE 8, e68739 (2013).

52. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 20782079 (2009).

53. Aronesty, E. ea-utils: Command-line tools for processing biological sequencing data. https://github.com/ExpressionAnalysis/ea-utils (2011).

54. Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 24602461 (2010).

55. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403410 (1990).

56. Cole, J. R. et al. The Ribosomal Database Project: improved alignments and new tools for rRNA analysis. Nucleic Acids Res. 37, D141D145 (2009).

57. Hippe, B. et al. Quantication of butyryl CoA:acetate CoA-transferase genes reveals different butyrate production capacity in individuals according to diet and age. FEMS Microbiol. Lett. 316, 130135 (2011).

58. Louis, P. & Flint, H. J. Development of a semiquantitative degenerate real-time pcr-based assay for estimation of numbers of butyryl-coenzyme A (CoA) CoA transferase genes in complex bacterial samples. Appl. Environ. Microbiol. 73, 20092012 (2007).

59. Walker, A. W. et al. Dominant and diet-responsive groups of bacteria within the human colonic microbiota. ISME J. 5, 220230 (2011).

60. Faust, K. et al. Microbial co-occurrence relationships in the human microbiome. PLoS Comput. Biol. 8, e1002606 (2012).

61. Shannon, P. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 24982504 (2003).62. Wisniewski, J. R., Zougman, A., Nagaraj, N. & Mann, M. Universal sample preparation method for proteome analysis. Nat. Methods 6, 359362 (2009).

63. Stahl, M. et al. L-fucose utilization provides Campylobacter jejuni with a competitive advantage. Proc. Natl Acad. Sci. USA 108, 71947199 (2011).

64. Palmer, C., Bik, E. M., DiGiulio, D. B., Relman, D. A. & Brown, P. O. Development of the human infant intestinal microbiota. PLoS Biol. 5, e177 (2007).

65. Ye, J. et al. Primer-BLAST: a tool to design target-specic primers for polymerase chain reaction. BMC Bioinformatics 13, 134 (2012).

66. Arthur, J. C. et al. Intestinal inammation targets cancer-inducing activity of the microbiota. Science 338, 120123 (2012).

67. Sun, X., Threadgill, D. & Jobin, C. Campylobacter jejuni induces colitis through activation of mammalian target of rapamycin signaling. Gastroenterology 142, 8695 (2012).

68. Lippert, E. et al. Gnotobiotic IL-10; NF-kappaB mice develop rapid and severe colitis following Campylobacter jejuni infection. PLoS ONE 4, e7413 (2009).

NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 13

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13419

Acknowledgements

This work was funded by the Government of Canada through Genome Canada and the Ontario Genomics Institute (OGI-067), CIHR grant number GPH-129340, CIHR grant number MOP-114872 and the Ontario Ministry of Economic Development and Innovation (REG1-4450). We also acknowledge funding from the IBD foundation, Crohns and Colitis Canada (CCC), the CHEO foundation, the Childrens Hospital of Eastern Ontario Research Institute and the Faculty of Medicine of the University of Ottawa. D.F. acknowledges a Canada Research Chair in Proteomics and Systems Biology and funding from CIHR grant number 117960-1. C.J. acknowledges funding from NIH (R01 DK73338, R01 DK47700 and R21 CA195226). We thank Maureen Bower of the National Gnotobiotic Rodent Resource Center at the University of North Carolina for assistance with gnotobiotic mice (NIH 5P40OD010995, CGIBD grant P30 DK034987). We thank A. Mack, M. Boland, C. Jimenez-Rivera and J. Barkey for their assistance in clinical data collection and making this research project possible. W.M. acknowledges funding from Mansoura University, Egypt. T.A. was supported by a scholarship from King Abdulaziz University, through the Saudi Arabian Cultural Bureau in Canada. M.M. acknowledges funding from CCFA. E.I.B. was supported by a Career Development Award from the Canadian Child Health Clinician Scientist Program, and a New Investigator Award from the Canadian Institutes of Health Research, Canadian Association of Gastroenterology, and Crohns and Colitis Canada. J.B. is the recipient of a CIHR/CAG/CCC Postdoctoral Award.

Author contributions

W.M., C.-K.C. and M.M. were involved in the experimental plan and data analysis. W.M. processed the biological samples and generated the microbiota data, C.-K.C. and A.E.S. generated and analysed the proteomic data, M.M. performed the mice experiments. S.A.D. and H.Z. helped perform the proteomic analysis. T.A., J.B., S.S. and M.H. participated in characterizing the microbiota. A.B. performed the pure culture experiments. E.B. and R.S. participated in patient diagnosis and enrolment, sample isolation and collection of clinical data. A.S., D.R.M., D.F. and C.J. contributed equally to the work. A.S. was involved in all the aspects of the study from experimental design to data analysis. D.M. recruited the patients, performed samples collection and clinical study. D.F. was involved in the experimental design and analysis of the proteomes. C.J. was

involved in the experimental design and analysis of the mice experiments. W.M., J.B. and A.S. wrote the paper. All the authors have revised and approved the manuscript submission.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: A.S., D.M., D.F., W.M., T.A. and C.-K.C. have a corresponding patent application Methods for the diagnosis and treatment of inammatory bowel disease led in 2014 (WO/2014/138999). A.S., D.M. and D.F. have co-founded Biotagenics, a clinical microbiomics company. The remaining authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Mottawea, W. et al. Altered intestinal microbiotahost mitochondria crosstalk in new onset Crohns disease. Nat. Commun. 7, 13419 doi: 10.1038/ncomms13419 (2016).

Publishers note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afliations.

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

14 NATURE COMMUNICATIONS | 7:13419 | DOI: 10.1038/ncomms13419 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications