Inheriting microbiome variation In mice, the bacterial species that persist within the gut microbiome are maternally inherited. However, maternally inherited variations in mitochondrial DNA (mtDNA) sequence also correlate with gut microbiome diversity, as well as the production of reactive oxygen species (ROS). In mice with mtDNA variants associated with increased production of ROS, Yardeni et al. found reduced gut microbiome species diversity. When pups were cross-fostered to unlink inheritance of mtDNA and gut microbiota, after weaning, the gut microbiome species were reflective of inherited mtDNA variation. Both pharmacological and genetic reduction of mitochondrial ROS abundance increased microbiome diversity. These data suggest that microbiome diversity is genetically encoded, and further imply that antioxidants may improve the efficacy of cancer immunotherapy, which is sensitive to microbiome composition.

Abstract Changes in the gut microbiota and the mitochondrial genome are both linked with the development of disease. To investigate why, we examined the gut microbiota of mice harboring various mutations in genes that alter mitochondrial function. These studies revealed that mitochondrial genetic variations altered the composition of the gut microbiota community. In cross-fostering studies, we found that although the initial microbiota community of newborn mice was that obtained from the nursing mother, the microbiota community progressed toward that characteristic of the microbiome of unfostered pups of the same genotype within 2 months. Analysis of the mitochondrial DNA variants associated with altered gut microbiota suggested that microbiome species diversity correlated with host reactive oxygen species (ROS) production. To determine whether the abundance of ROS could alter the gut microbiota, mice were aged, treated with N-acetylcysteine, or engineered to express the ROS scavenger catalase specifically within the mitochondria. All three conditions altered the microbiota from that initially established. Thus, these data suggest that the mitochondrial genotype modulates both ROS production and the species diversity of the gut microbiome, implying that the connection between the gut microbiome and common disease phenotypes might be due to underlying changes in mitochondrial function.

INTRODUCTION The human gut contains trillions of microorganisms that coincide with an individual’s health status (1, 2). Changes in this gut microbiota community correlate with a variety of clinical phenotypes including diabetes mellitus (3), autism (4, 5), and Parkinson disease (6). A subset of these phenotypes is modified in animal models by gut microbiome transfer (5, 6). Concurrently, a wide range of human (7) and mouse (8) metabolic and degenerative diseases also correlates with mitochondrial genotypes, including diabetes mellitus (9–11), autism (9, 12, 13), and Parkinson disease (14). Moreover, changes in the microbiota correlate with changes in mitochondrial metabolism (15), human mitochondrial DNA (mtDNA) haplogroup variants (16), and a mouse mtDNA variant and haplogroup (17). The question then becomes does the microbiome cause the disease or do mitochondrial alterations determine both the microbiome composition and the animal phenotype? Both mtDNA (18) and the microbiome community (19) are vertically inherited from the mother. In wild-trapped mice, the microbiome is stable over multiple generations, suggesting that vertical inheritance is the dominant mode of gut bacterial transmission (19). Because newborn pups receive their initial gut inoculum from their nursing mother (20), vertical transmission may result from repeated maternal microbiota inoculation or from the transmission of the mother’s mtDNA. Each cell contains hundreds of mitochondria and maternally inherited mtDNAs. Mitochondria not only produce most of the cellular energy by the process of oxidative phosphorylation (OXPHOS) but also are integral to most of the cell’s metabolic pathways, regulate cytosolic Ca++, modulate cellular redox status, and generate much of the cell’s reactive oxygen species (ROS) (21). The mtDNA encodes 13 critical subunits of the OXPHOS complexes plus the 12S and 16S ribosomal RNAs (rRNAs) and the 22 transfer RNAs (tRNAs) for mitochondrial protein synthesis. Of the 13 mtDNA subunits, 7 encode polypeptides for complex I [reduced form of nicotinamide adenine dinucleotide (NAD+) (NADH) dehydrogenase], 1 for complex III (bc 1 complex), 3 for complex IV (cytochrome c oxidase), and 2 for complex V [adenosine 5′-triphosphate (ATP) synthase]. All other mitochondrial polypeptide genes are encoded within the nuclear DNA (nDNA) (22, 23). There are three phenotypically relevant classes of mtDNA mutations: maternally transmitted ancient adaptive polymorphisms, materially inherited recent deleterious mutations, and somatic mtDNA mutations that accumulate with age. The adaptive polymorphisms arose over thousands of years on radiating maternal mtDNA lineages. Those functional variants that were beneficial in a particular environment were enriched along with the linked mtDNA variants to generate a group of related haplotypes called a haplogroup. For example, a functional mutation that increased mitochondrial ROS could enhance innate immunity and be selected because of resistance to certain infections (23–25). This situation has been observed in the mouse by comparing the mitochondrial ROS production of mice with the same nucleus that harbors mtDNAs from NZB and 129 or Balb/c mice lineages (26–28). Recent deleterious mutations arise along material lineages and can result in disease. Deleterious mutations can either be homoplasmic (pure mutant) or heteroplasmic (mixed mutant and normal) (23). As an example, a heteroplasmic missense mutation in the mtDNA ND6 gene at m.14600G>A (P25L) results in neurological disease (29), which is recapitulated in the equivalent mouse ND6 nt 13997 m.G A (P25L) mutation (30). mtDNA damage and mutations accumulate throughout life, eroding energy metabolism and increasing ROS production (31–33). Last, mutations in nDNA-coded mitochondrial genes can also result in mitochondrial dysfunction and phenotypic manifestations (23). Deleterious mtDNA or nDNA variants that partially impede mitochondrial OXPHOS result in redirection of electrons from the electron transport chain that would normally reduce O 2 to 2H 2 O to add a single electron to O 2 , thus generating superoxide anion. Superoxide anion is the first of the ROS species, which can be converted to hydrogen peroxide (H 2 O 2 ) by Mn superoxide dismutase located in the mitochondrial matrix (21). Addition of mitochondrially targeted catalase (mCAT) can remove the mitochondrial H 2 O 2 , protect the mtDNA, and extend the lifespan (31). The mtDNAs of 129 and NZB mice differ in multiple nucleotide positions, but the most notable difference has been reported to be a polymorphism in the mtDNA DHU loop of the mtDNA tRNAArg gene. This polymorphism at nucleotide (nt) 9821 involves a homopolymer of A’s, whose length ranges from 8 to 10 A’s. The NZB mtDNA tRNAArg has 10 A’s, but 129 mtDNA has 9 A’s. The addition of extra A’s has been observed to increase mitochondrial ROS production (26, 27). We now report the use of these mice to determine the role of the mitochondrial function on the gut microbiome composition. We found that the mitochondrial genotype of mice correlates with the gut microbiota species diversity. Even when the initial microbiota community was obtained from a foster mother with different microbiota community by cross-fostering, the microbiota community progressed toward that characteristic of the pup’s genotype within 2 months. Last, our data suggest that one important factor in the mitochondrial modulation of the gut microbiome is mitochondrial redox status and associated ROS production.

DISCUSSION The gut microbiome composition correlates with clinical presentations such as diabetes mellitus (22, 23), autism (4, 5), and Parkinson disease, as does mtDNA variation in these same diseases (9, 12–14, 37). We now show that the mtDNA genotype correlates with the gut microbiota and by cross-fostering experiments that the offspring’s mitochondrial redox state and ROS abundance can restructure the gut microbiome, even when initially derived from a fostering mother who transmits a different gut microbiota community. From these observations, we can conclude that modification of the mitochondrial redox status and associated ROS production is associated with changes in the gut microbiota community, implicating the mitochondrial function in controlling the composition of the microbiome community. The question remains: How does the mitochondrial redox state and ROS production modulate the gut microbiome community? One possibility is that the mitochondrial redox status might be modulating the gut microbiome through mitochondrial metabolites. A marked difference was observed in the NADH lifetime of C57BL/6 versus C57BL/6J, demonstrating a clear difference in the NAD+/NADH ratios in their intestinal tissues. Because the NAD+/NADH ratio is central to regulating metabolic flux through the mitochondrion, this could alter mitochondrial intermediates that could diffuse into the gut and select for altered microbiome populations. We have found that changes in the nuclear and mitochondrial NADH lifetimes correlate with changes in mitochondrial metabolites, which, in turn, modulate the epigenome (38). Another possibility is that altered mitochondrial redox signaling could affect the function of the intestinal epithelial, neuronal, and smooth muscle cells. Because mitochondrial genetic defects are commonly associated with intestinal dysmotility, this could also affect the gut microbiota community. Similarly, mitochondrial function is central to the function of the immune system (24). The mild mitochondrial dysfunction associated with the ND6P25L mtDNA mutation preferentially impairs the more oxidative regulatory T cells, leading to loss of inhibition of the more glycolytic effector T cells (39). Reductions in T regulatory cells would increase the inflammatory response, which could affect the gut microbial community. Mitochondria are also known to be central to the innate immune system as well. Release of mtDNA in stressed cells can activate the cGAS-Sting pathway, which regulates interferon production (24). Similarly, mtDNA oxidation within macrophages can activate the NLRP3 inflammasome, causing nuclear factor κB (NFκB) activation and stimulation of systemic inflammation (25). Perhaps mitochondrial variation might mediate the differential effects of gut microbiota inoculation on the therapeutic efficacy in metastatic melanoma patients treated with anti–PD-1 (programmed cell death protein-1) therapy (40–43). Mitochondrial NADH also regulates mitochondrial H 2 O 2 . NADH is converted to NADPH (reduced form of nicotinamide adenine dinucleotide phosphate) in the mitochondrion by the nicotinamide nucleotide transhydrogenase using the mitochondrial inner membrane potential as the added source of energy. NADPH is used to reduce oxidized glutathione (GSSG) to 2GSH, and GSH is used to reduce H 2 O 2 to 2 H 2 O by glutathione peroxidase. Hence, changes in NAD+/NADH ratio could also modulate H 2 O 2 amounts. Conversely, mitochondrial H 2 O 2 might act directly on specific gut microbiota. This possibility is supported by the fact that simply expressing mCAT at the mitochondrion had the single greatest impact on the gut microbiome. The function of mCAT is to convert mitochondrial H 2 O 2 to H 2 O and O 2 . Because H 2 O 2 tends to react with other molecules rapidly, if it is acting as a mitochondria-to-microbiota messenger, it would need to be generated in close proximity to the gut microbiota. This expectation is consistent with our observation that the expression of mCAT mRNA was 10 times higher in the gut than in other somatic tissues in the mCAT mice. Because mitochondrial H 2 O 2 detoxification requires both NADPH and GSH, alterations in H 2 O 2 detoxification by mCAT could also account for the observed changes in gut NADH lifetimes and the effect of adding NAC to generate more GSH. Undoubtedly, there are other potential mitochondrial–gut microbiota signal pathways, perhaps even harkening back to the use of peptides and small molecules in quorum sensing carried over from the bacterial origin of the mitochondrion. Regardless of the molecular basis of the specific mitochondrial molecular signals that are regulating the gut microbiome that correlate with H 2 O 2 production and NAD+/NADH redox status, our data provide substantial support for the hypothesis that the correlation between the gut microbiome composition and the range of clinical manifestations reported is because both the clinical phenotypes and the gut microbiota are regulated by the mitochondrial genotype and associated functions. This conclusion suggests new approaches for treating both complex diseases and gut microbiota dysbiosis.

MATERIALS AND METHODS Mice The Institutional Animal Care and Use Committee from the Children’s Hospital of Philadelphia approved all protocols, and the protocols comply with all relevant ethical regulations regarding animal research. The mice were fed a 5LOD diet from PicoLab Laboratory and were maintained on a 13-hour/11-hour light-dark cycle. Seven mouse models were used for this study (Table 1); six were on C57BL/6J background (44–46): B6 mtDNA, ND6 m.13997 G>A (ND6P25L) (30), 129 mtDNA homoplasmic, NZB mtDNA homoplasmic, NZB/129 mtDNA heteroplasmic (35), and transgenic (Tg)-mCAT (31). The last strain was C57BL/6EiJ. All mice, except for the Tg-mCAT, were bred and maintained in the Wallace laboratory colony for more than 10 years. The Tg-mCAT strain was purchased from The Jackson Laboratory and bred to the Wallace laboratory C57BL/6J strain. To determine the effect on NAC, 2-month-old C57BL/6J mice were divided into two groups: one group was treated with NAC (2 g/kg body weight per day) (Sigma-Aldrich, St. Louis, MO) in the drinking water for 10 weeks, and the second group served as control. Feces collection Fecal samples were collected from male mice at 2 to 3 months or 5 months of age in the morning. Individual mice were moved to an empty clean cage, and fresh fecal pellets were collected into a 1.5-ml Eppendorf tube and placed on dry ice. All samples were stored at −80°C until DNA extraction. Cross-fostering Cross-fostering experiments were performed using C57BL/6EiJ and C57BL/6J mice. Breeding pairs of C57BL/6EiJ and C57BL/6J mice were set up at the same time. Fourteen days after introduction, females were monitored daily for pregnancy stage and the males were removed from pregnant females. The entire litter of the newly born pups was transferred during the first 24 hours of life to a nursing mother from the opposite strain. After 3 weeks, all pups were weaned. Feces were collected when pups were 2 months old (20). DNA extraction DNA was extracted from frozen mouse pellets using the MO BIO PowerSoil HTP DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA, USA). Sequencing libraries were generated by polymerase chain reaction (PCR) amplifying the V1 and V2 regions of the 16S rRNA gene using barcoded universal primers 27F (5′-AGAGTTTGATCCTGGCTCAG-3′) and 338R (5′-3TGCTGCCTCCCGTAGGAGT-3′) (47, 48). For each sample, PCRs were performed in quadruplicate with AccuPrime Taq High Fidelity (Invitrogen, Carlsbad, CA, USA). Individual reactions were combined, purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA), and then pooled with other samples in equimolar quantities. Extraction blanks and DNA-free water were subjected to the same extraction and amplification procedures to allow for empirical assessment of environmental and reagent contamination. Positive controls, consisting of eight artificial 16S gene fragments synthesized in gene blocks (Integrated DNA Technologies, Coralville, IA, USA) and combined in known abundances, were also included. Libraries were sequenced on the Illumina MiSeq to obtain 250-base pair (bp) paired-end reads. Bioinformatics processing Sequence data were processed using QIIME version 1.9.1 (49). Read pairs were quality-filtered and joined using default parameters. Sequences were grouped into OTUs using UCLUST, with a 97% sequence identity threshold setting (50). Taxonomic assignments were generated using the default method in QIIME, which selects among top matches in the Greengenes reference database (51). Representative sequences from each OTU were aligned using PyNAST (49), and a phylogenetic tree was inferred from the multiple sequence alignment using FastTree (52). The Shannon diversity of each sample was calculated with the formula H ′ = − ∑ i = 1 N p i ln p i where p i is the relative abundance of the ith organism in the sample. Similarity between samples was assessed by weighted and unweighted UniFrac distance (53, 54). Statistical analysis Data files from QIIME were analyzed in the R environment. Because mice sharing the same cage harbor similar microbiota (55), the cage of each mouse was modeled as a random effect in the statistical analysis. The alpha diversity differences between the genotypes were assessed with linear mixed-effects models. The PERMANOVA method was used to test for differences in bacterial community composition, as quantified by unweighted and weighted UniFrac distance between samples (56). Relative abundance was assessed using the most specific taxonomic assignment available for each OTU. Bacteroidetes/Firmicutes ratios were calculated by summing up the relative abundance at the phylum amounts. Differences between groups were tested using linear mixed-effects models on log-transformed ratios. Taxa were selected for testing if the mean abundance exceeded 1% among the samples to be analyzed. Taxon abundance was tested using generalized linear mixed-effects models. To account for multiple comparisons of taxon relative abundance, P values were adjusted to control for a false discovery rate of 5%. Each filled or empty circle in a graph represents an individual mouse. Measuring ROS production from mouse strain livers Mitochondrial H 2 O 2 production was assayed by Amplex Red. Whole liver homogenates, clarified at 20,000g for 10 min in radioimmunoprecipitation assay buffer and protease inhibitor, were assayed for comparison between C57BL/6J mice with and without NAC treatment and with and without the mCAT transgene. Isolated liver mitochondria were used and assayed for the C57BL/6EiJ versus C57BL/6J comparison. Mice were euthanized by cervical dislocation, and the livers were quickly removed and placed on ice. Livers were washed with the isolation buffer [215 mM mannitol, 75 mM sucrose, 0.1% bovine serum albumin (BSA), 1 mM EGTA, and 20 mM Hepes (Na+) (pH 7.2)] and then transferred to 10 ml of isolation buffer for trimming and glass Dounce homogenization. For mitochondrial isolation, unbroken cells and nuclei were removed by centrifugation at 1000g for 10 min, and the mitochondria were pelleted from the supernatant by centrifugation at 20,000g for 10 min. The mitochondrial pellet was washed in isolation buffer without EGTA (10,000g for 10 min) and resuspended in isolation buffer without EGTA and BSA. Protein amounts were determined using bicinchoninic acid (BCA) assay (Pierce). One hundred micrograms of isolated liver homogenate or isolated mitochondria was assayed using the Amplex Red Hydrogen Peroxide/Peroxidase Assay Kit (Invitrogen, catalog no. A22188). Amplex Red fluorescence was measured with an excitation wavelength of 530 nm and a 590-nm emission filter. FLIM of NADH autofluorescence from mouse small intestine Cross sections of freshly isolated small intestine were imaged in Tyrodes buffer (pH 7.4) using LSM710 (Zeiss) equipped with a time-correlated single-photon counting module (HPM-100-40 and SPCM 9.81, Becker and Hickl) (57). NADH was excited by a femtosecond-pulsed two-photon laser (Coherent) at 730 nm, and its autofluorescence signal was detected through a 680-nm short-pass and 460/50-nm band-pass emission filter. To cover the whole cross section, a 20× lens (numerical aperture, 0.8) was used in combination with the tile imaging function (4 × 4 tiles, 2121.2 μm by 2121.2 μm). FLIM images were analyzed in SPCImage 7.4 using a biexponential decay model with T1, T2, and the shift not specified. Non-NADH autofluorescence was reduced by setting the minimum lifetime to 200 ps, and the maximum Chi2 to 5. The average mean NADH lifetime (T mean ) of the image sections was quantified. RNA isolation and real-time qRT-PCR of human catalase mRNA Total RNA was isolated using TRIzol reagent (Thermo Fisher Scientific) according to the manufacturer’s protocols. Briefly, mouse was sacrificed by cervical dislocation, and tissues (small intestinal and liver) were flash-frozen on liquid nitrogen and kept at −80°C. The tissues were homogenized using a motorized blade homogenizer (Polytron) in 1.0 ml of TRIzol. Total RNA was resuspended in 50 μl of ribonuclease-free H 2 O. Contaminating DNA was removed from RNA using the TURBO DNA-free Kit from Ambion. Approximately 3.0 μl of total RNA was treated in a 10-μl reaction according to the manufacturer’s protocol. First-strand complementary DNA (cDNA) was created from total RNA using SuperScript IV First-Strand Synthesis System (Thermo Fisher Scientific). The entire deoxyribonuclease-treated RNA (10 μl) was used for first-strand synthesis with oligo-dT primers. TaqMan-based qRT-PCR was performed using TaqMan Expression Assays (Thermo Fisher Scientific) on the Viaa7 platform. Assay was Hs00937387_m1 specific human catalase mRNA, and the control (housekeeping) assay was mouse HPRT Mm00446968_m1. qRT-PCRs were 5 μl of 2× Universal Master Mix (no UNG), 1 μl of first-strand cDNA, 0.5 μl of TaqMan assay, and 3.5 μl of H 2 O. All reactions were run in triplicate on a 384-well plate.

http://www.sciencemag.org/about/science-licenses-journal-article-reuse This is an article distributed under the terms of the Science Journals Default License.

Acknowledgments: We thank A. Butic, S. L. Royston, R. Morrow, K. L. Mitchell, and J. A. Tintos (Wallace laboratory) and C. Hofstaedter, D. Kim, M. Moraskie, and H. Zhang (CHOP Microbiome Center, Children’s Hospital of Philadelphia) for technical assistance and B. Boursi (Center for Clinical Epidemiology and Biostatistics, Perelman School of Medicine at the University of Pennsylvania, Philadelphia) for fruitful discussions. All authors acknowledge the essential contribution of mouse models. Funding: This work was supported by the PennCHOP Microbiome Pilot Grant and NIH grants MH108592, NS021328, MN110285, and DO W81XWH-16-1-0401 (to D.C.W.) and the PennCHOP Microbiome program (to G.D.W.). Author contributions: T.Y. and D.C.W. directed the project. T.Y., D.G.M., and D.C.W. wrote the manuscript. T.Y., C.E.T., K.B., D.G.M., and D.C.W. were responsible for experimental design. T.Y., L.M.M., and C.E.T. performed the experiments. T.Y., C.E.T., K.B., P.M.S., L.N.S., G.D.W., and D.C.W. performed data analysis. Competing interests: The authors declare that they have no competing financial interests. Data and materials availability: The DNA sequence datasets generated during and/or analyzed during the current study are available in the NCBI SRA repository (BioProject ID PRJNA423318). All other data needed to evaluate the conclusions in the paper are present in the paper.