Microbiota contribute to the induction of type 2 diabetes by high-fat/high-sugar (HFHS) diet, but which organs/pathways are impacted by microbiota remain unknown. Using multiorgan network and transkingdom analyses, we found that microbiota-dependent impairment of OXPHOS/mitochondria in white adipose tissue (WAT) plays a primary role in regulating systemic glucose metabolism. The follow-up analysis established that Mmp12+ macrophages link microbiota-dependent inflammation and OXPHOS damage in WAT. Moreover, the molecular signature of Mmp12+ macrophages in WAT was associated with insulin resistance in obese patients. Next, we tested the functional effects of MMP12 and found that Mmp12 genetic deficiency or MMP12 inhibition improved glucose metabolism in conventional, but not in germ-free mice. MMP12 treatment induced insulin resistance in adipocytes. TLR2-ligands present in Oscillibacter valericigenes bacteria, which are expanded by HFHS, induce Mmp12 in WAT macrophages in a MYD88-ATF3–dependent manner. Thus, HFHS induces Mmp12+ macrophages and MMP12, representing a microbiota-dependent bridge between inflammation and mitochondrial damage in WAT and causing insulin resistance.
A combination of overconsumption of food high in simple sugars and animal fat diet and a sedentary lifestyle provoked dramatic increases in the rates of diabetes, obesity, and other metabolic diseases over recent decades. Type 2 diabetes (T2D) involves dysfunction in multiple organs, with critical events in the liver, white adipose tissue (WAT), and muscles, promoting the development of systemic insulin resistance (IR). Despite established knowledge of the disease processes occurring within each organ, their relative contribution to systemic metabolic dysregulation is less clear. Furthermore, gut microbiota are emerging as key players in multiple proposed mechanisms that affect different organs (Gurung et al., 2020; Hartstra et al., 2015; Sonnenburg and Backhed, 2016). Nevertheless, there is no consensus on which of these organs and pathways are most important for the development of IR and T2D. To recapitulate human diet–induced T2D in animal models, special diets were developed, including high-fat/high-sugar diet (HFHS), frequently referred to as western diet (WD; Small et al., 2018).
Several studies interrogating different omics data have been conducted attempting to address these questions (Hosseinkhani et al., 2021; Passaro et al., 2021). However, commonly used association analyses between diet, microbiome composition, and disease do not infer causality or attempt to establish master regulators/pathways that can be affected by inverse causality and other confounding factors. In this study, we used multiorgan network analysis and transkingdom causal inference employing statistical dependencies between differentially expressed elements (e.g., genes, microbes, and phenotypic markers) combined with causality principles (Pearl et al., 2009; Yambartsev et al., 2016), information theory (Shannon, 1953; Uda, 2020), and interventions (e.g., diet, antibiotics, etc.). This system’s approach provided us with two causal inferences: first, microbiota-dependent inflammation drives IR in WAT; second, this inflammation is caused by specific members of HFHS-enriched gut microbiota.
Specifically, we investigated pathways underlying the role of diet-modified microbiota in metabolic disease by comprehensively characterizing HFHS-induced changes in four metabolic organs in mice: intestine, liver, WAT, and muscle. We found that a novel predominant mechanism contributing to systemic IR is microbiota-dependent damage of mitochondrial oxidative phosphorylation (OXPHOS) in WAT. To establish how microbiota contributes to this process, we characterized IR-promoting microbiota-induced macrophages in WAT of mice and humans. Using network analysis coupled with loss- and gain-of-function studies, we established MMP12 as a critical effector of these cells. Finally, we identified a specific microbe (Oscillibacter valericigenes [OV]) that promotes Mmp12 expression through the production of TLR2 ligands. We demonstrate in vitro that Mmp12 is induced by stimulation of the TLR2/Myd88 signaling pathway and activation of the transcription factor (TF) ATF3.
HFHS diet drives alterations predominantly in a microbiota-dependent manner
To dissect the contribution of microbiota to the effects of WD (designated herein HFHS), specific pathogen–free (SPF) and germfree (GF) C57BL/6 mice were fed a HFHS or matched normal diet (ND) for 8 wk (Fig. 1 A). As expected (Backhed et al., 2007; Rabot et al., 2010), HFHS increased glucose intolerance in GF mice and the presence of microbiota further exacerbated IR and adiposity in SPF mice (Fig. 1 B and Table S1).
Next, using samples from SPF and GF mice, we analyzed whole-tissue transcriptomes of the intestine, as it harbors the most abundant microbial community in mammals, and of the three other major insulin-responsive organs implicated in T2D (liver, WAT, muscle). We found that WAT and liver transcriptomes were more sensitive to HFHS (2,427 and 794 differentially expressed genes [DEGs], respectively), compared to muscle and ileum (170 and 116 genes, respectively). The finding that more genes are regulated in WAT than in liver or muscle might be expected because mice that have consumed the HFHS for 8 wk represent the initial stages of T2D, with molecular changes progressing from adipose to the liver to the muscle (Kraegen et al., 1991; Perry et al., 2015). HFHS-regulated genes were classified as microbiota dependent (category 1) or microbiota independent (category 2; Fig. 1 C), while a small fraction of genes could not be assigned to either of the above categories. Microbiota-dependent genes were further classified as regulated by HFHS differently in SPF and GF mice (category 1A; e.g., Itgax) or regulated by diet in both types of mice, but whose expression was enhanced by the microbiota (category 1B; e.g., S100b; Fig. 1 C). We found that microbiota play a key role in HFHS-modulated gene expression in all four organs, with about 87% of genes in the liver, 76% in WAT, 57% in the ileum, and 44% in muscle belonging to category 1 (Fig. 1 D and Table S2).
Multiorgan network analysis identifies OXPHOS/mitochondria in WAT as a key driver of systemic metabolism
To model HFHS-induced alterations, we reconstructed a co-variation network, integrating gene expression with systemic metabolic parameters (Fig. 2 A). Approximately two thirds of DEGs in the four organs were retained in the network (Fig. 2 B and Table S3) after filtering out genes (nodes) and edges (connections between nodes of network) that did not fulfill statistical and causality criteria (Yambartsev et al., 2016). As expected, the network node degree distribution conformed with the power-law structure of regulatory biological networks (Fig. 2 C). Two subnetworks were identified in the liver and WAT (Fig. 2 A). In WAT, the larger subnetwork mainly reflected impaired mitochondrial functions, in particular OXPHOS (Fig. 2 D ). The smaller WAT subnetwork with mostly upregulated genes was enriched for inflammation and immune response genes (Fig. 2 D ). The larger liver subnetwork included mostly HFHS downregulated genes that were enriched for ribosomal components (Fig. 2 D ), while the other liver subnetwork, largely upregulated by HFHS, was represented by fatty acid metabolism genes (Fig. 2 D ). The genes controlling fatty acid metabolism in the liver and inflammation in fat were those most impacted by microbes (>85% of genes from the 1A and 1B categories; Fig. 2 A).
To disclose the regulatory relationships between the pathways and systemic metabolic parameters, we investigated network topology reflective of information flow in the network (including shortest path [Fig. 2 E] and bipartite betweenness centrality [BiBC; Fig. 2 F and Table S4]). This analysis showed (blue arrows in Fig. 2 A) that among the four subnetworks, mitochondrial OXPHOS in WAT directly affects systemic glucose metabolism, while fatty acid pathways in the liver and inflammation in WAT primarily act indirectly through microbiota-dependent damage of mitochondrial OXPHOS in the adipose (Fig. 2, A and G).
Three independent lines of evidence further supported the computational inference that IR is induced in WAT and is partially dependent on the gut microbiota. First, HFHS induced an IR gene signature (Jung et al., 2018) in WAT of SPF, but not of GF mice, indicating a significant contribution of microbiota (Fig. 3 A). Second, phosphorylated AKT and—to a lesser extent—p38 MAP kinase, both classical markers of insulin signaling (IS), negatively correlated with inflammation and positively with OXPHOS genes in WAT of SPF, but not GF mice (Fig. 3 B; Fig. S1, A and B; and Table S5). Third, genes with a high regulatory capacity (i.e., high BiBC [Dong et al., 2015] to systemic metabolic markers in the network) and upregulated by HFHS negatively correlated with the response of adipocytes to insulin in vitro (Figs. 3 C and S1 C; and Table S6). This last result indicated that a large proportion of genes inferred to mediate the negative effects of HFHS on systemic metabolism are involved in insulin sensitivity in adipocytes.
Thus far, our experiments coupled with network analysis predict that HFHS likely acts in WAT by driving microbiota-related OXPHOS damage, leading to systemic IR.
Mmp12-positive macrophages link microbiota-dependent inflammation and IR in WAT
Since our results indicated that WAT is the most critical tissue initially involved in HFHS-driven T2D, we focused further works on this tissue. Our multiorgan network model also inferred that microbiota-dependent activation of immune/inflammation pathways leads to OXPHOS damage and consequently systemic IR.
To identify microbiota-dependent genes that mediate this process, we interrogated the network for genes with high BiBC (Dong et al., 2015), reflecting their ability to mediate inflammatory damage on the OXPHOS pathway. Mmp12 was identified as the most prominent microbiota-dependent gene (category 1A) connecting inflammation and mitochondrial OXPHOS subnetworks in WAT (Fig. 4 A). Notably, HFHS induced Mmp12 expression in the presence of both mouse and human microbiota (Fig. 4 B). Mmp12 also negatively correlated with p-AKT/AKT (Fig. 4 C), further supporting its role in WAT-specific IR.
In WAT, Mmp12 is predominantly expressed in CD45+ hematopoietic cells (Fig. 4 D). Using meta-analysis of eight publicly available single-cell RNA sequencing (scRNA-seq) datasets (Table S7), we established a 95-gene signature of Mmp12+ CD45+ cells (Figs. 4 E and S2 A; and Table S8) and confirmed its expression in macrophages (Fig. S2 B; Lee et al., 2014; Martinez-Santibanez et al., 2015). This new Mmp12 meta-signature comprehensively defined metabolic macrophages described by three independent reports (Hill et al., 2018; Jaitin et al., 2019; Silva et al., 2019; Fig. 4 F). Although these studies used different sets of markers, they are all part of our newly discovered meta-signature together with other genes that were previously implicated in IR in WAT (Cd9 [Hill et al., 2018], Cd36 [Kennedy et al., 2011], Lgals3 [Li et al., 2016; Pejnovic et al., 2013], and Itgax [Patsouris et al., 2008]). This signature identifies a unique cell population—IR-associated adipose tissue macrophages (IR-ATMs; Fig. 4 E)—that was strongly increased by HFHS in SPF, but not in GF mice (Figs. 4 G and S2 C). These data demonstrate that HFHS increases not only Mmp12 expression but also IR-ATMs in WAT in a microbiota-dependent manner (Figs. 4 G and S2 C). Altogether, these results suggest that the Mmp12 meta-signature is characteristic of WAT-infiltrating, metabolically active macrophages with a negative impact on glucose metabolism.
MMP12 inhibits the OXPHOS/mitochondria pathway in WAT, worsening systemic glucose metabolism
While several genes in the IR-ATM signature were previously associated with ATM metabolic potential, our analysis predicted that MMP12 is an effector molecule contributing to impairment of mitochondrial OXPHOS in adipocytes and systemic glucose metabolism. Accordingly, we observed that genetic deletion of Mmp12 or inhibition of MMP12 protein activity by the specific inhibitor MMP408 (Li et al., 2012; Li et al., 2009; Long et al., 2015) in HFHS-fed SPF mice improved fasting glucose levels, glucose tolerance, and IR (Fig. 5, A and B, left panel; and Fig. S2 D). The most pronounced transcriptional effect of MMP408 was the restoration of mitochondria-related pathways (Fig. 5 C), including increased expression of genes responsible for various mitochondrial functions (Ndufaf5, Acacb, Maob, Lpin1, Nat8l, and Mrpl38; Table S9), and of antioxidative genes (Fig. 5 D). This was accompanied by a decrease in carbonyl protein, a marker of oxidative stress (Fig. 5 D, right panel). MMP12 inhibition had a significantly higher impact on microbiota-dependent than on microbiota-independent genes in the OXPHOS/mitochondria subnetwork (Fig. 5 E), and improved glucose metabolism in HFHS-fed SPF mice but not in GF mice (Figs. 5 B and S2 D). Finally, mitochondrial/OXPHOS-related gene expression that was decreased by HFHS and increased by MMP408 negatively correlated with markers of metabolic disease (Fig. 5 F and Table S10). These results confirm the negative impact of MMP12 on glucose metabolism and OXPHOS/mitochondria as well as its microbiota dependence.
Mitochondrial OXPHOS damage in adipocytes leads to an increase in the reactive oxygen species and is linked to systemic IR (Kusminski et al., 2012). Thus, we sought to identify whether adipocytes represent the target cells of MMP12. Treatment of adipocytes in vitro with activated recombinant MMP12 reduced insulin-stimulated glucose uptake by about 50% (Fig. 5 G) and reversed the insulin-induced transcriptional profile in adipocytes (Fig. 5 H and Table S11). On one hand, our in vivo data showed that MMP12 blockade leads to a reduction of oxidative markers (Fig. 5 D). On the other hand, mitochondrial OXPHOS dysregulation increases oxidative stress (Kusminski et al., 2012). Thus, we used an antioxidant (Tempol; Wilcox, 2010) to reduce oxidative stress in adipocytes treated with recombinant MMP12 (rMMP12) and observed partial restoration of insulin-stimulated glucose uptake in the presence of the antioxidant (Fig. 5 I).
Previous attempts to investigate the role of MMP12 in glucose metabolism using Mmp12 KO mice showed contradictory results, with some studies reporting phenotypes consistent with our results (Amor et al., 2016; Bauters et al., 2013; Lee et al., 2014; Martinez-Santibanez et al., 2015). In contrast to earlier research that used only KO mice, our study employed two “loss of function” (Mmp12 KO and chemical inhibitor) perturbations and one “gain of function” (rMMP12) perturbation, and showed full consistency across all three types of perturbations. The discrepancies between previous studies are likely to be explained by microbiota-dependent Mmp12 expression in WAT (Fig. 4, A and B) and the confounding effect of large differences between microbiota composition in mice of the same genetic background housed in different animal facilities (Rausch et al., 2016).
MMP12-dependent molecular signature in WAT is associated with IR in humans
To verify whether these findings in experimental animals translate into human pathology, we tested the Mmp12/IR-ATM meta-signature in patients with T2D. We first performed scRNA-seq analysis of the stromal vascular fraction (SVF) of WAT from patients with obesity and with or without T2D (Vijay et al., 2020). The IR-ATM meta-signature was overexpressed in WAT from obese individuals with T2D compared with those without T2D (Fig. 6 A), and predominantly in the MC1, a macrophage subtype with the highest increase in patients with T2D (Fig. 6, B–D; and Fig. S3 A). Meta-analysis of total WAT transcriptomes from five available patient cohorts (Hardy et al., 2011; Jain et al., 2013; Klimcakova et al., 2011; Soronen et al., 2012; Yang et al., 2019; three with similar body mass index in patients and controls) showed a consistent increase of Mmp12/IR-ATM signature (65/83 genes) in individuals with IR (Fig. 6, E and F; Fig. S3 B; and Table S12) and also highly overlapping with the results in mice (Fig. 6 G).
OV induces Mmp12/IR-ATM meta-signature in macrophages
To identify which bacterial taxa are responsible for the induction of Mmp12 in WAT, we performed transkingdom causal inference including a combination of standard techniques (such as instrumental variables and mediation analysis) incorporating additional aspects of principles of causality such as common cause (Pearl et al., 2009; Reichenbach et al., 1991), Simpson Paradox (Blyth, 1972), and correlation inequalities (Yambartsev et al., 2016). Specifically, we first perturbed microbiota of HFHS-fed SPF mice with six antibiotics for 4 wk. We next assessed the effects on the microbiome, Mmp12 expression in WAT and glucose metabolism markers (Fig. 7 A). Using these data, we inferred (see details of inference in Figs. S4 A and 7, B–D) that among 263 different taxa, three bacteria (Figs. 7 E and S4 B) were potential inducers of Mmp12 in HFHS-fed mice, with OV being the candidate with the largest impact (Fig. 7 E). To test this prediction, we first supplemented ND- and HFHS-fed mice for 4 wk with OV. We found that OV supplementation to mice increased the expression of several genes of the IR-ATM signature (Fig. 7 F), including Mmp12 in WAT (Fig. 7 G), with ∼70% similarity to human T2D (Fig. S4 C). OV supplementation also showed a trend of worsened glucose tolerance in mice with 2-h serum glucose levels being significantly increased (Fig. S4 D). Interestingly, another HFHS-associated bacterium (Romboutsia ilealis), reportedly a negative regulator of glucose metabolism (Rodrigues et al., 2021), was eliminated at the last step of the analysis (Fig. 7, D and E). Accordingly, supplementation of R. ilealis did not induce Mmp12 expression in mouse WAT (Fig. 7 G), indicating some level of specificity for OV.
To test whether OV can stimulate Mmp12 in macrophages directly or whether it requires cell–cell contact or circulating factors, we treated a macrophage cell line (RAW 264.7 cells) with heat-killed (HK) OV or cell-free supernatant (CFS). We found that the CFS from OV induced higher expression of Mmp12 (Fig. S5 A) and a large part of the IR-ATM meta-signature (Fig. 7 H) with significant overlap to observations in humans with IR (Fig. S5 B). These results offer another indication that macrophage-related inflammation in WAT of insulin-resistant patients is driven by microbiota.
TLR2-induced Mmp12 expression in macrophages is Atf3 dependent
To identify the mechanisms of Mmp12 induction by OV, we first tested (via knockdowns and/or knockouts) four major adaptor molecules (Myd88, Ticam-1, Mavs, Ripk2) that mediate signaling downstream of receptors of innate immunity and inflammasome activation. These experiments revealed that Mmp12 regulation by OV is Myd88 dependent (Figs. 8 A and S5 C). Next, we screened OV supernatant for TLR-ligands and found that this bacterium produces Tlr2 and Tlr5 agonists (Fig. 8 B). The follow-up studies demonstrated that Tlr2 agonist strongly induced Mmp12 in the macrophage cell line, while Tlr5 ligand shows only a trend in the same direction (Fig. 8 C).
While the classical view of Myd88-dependent signaling is focused on NF-κB activation, we investigated if any other factor might be involved in Mmp12 transcription in IR-ATMs. Analyzing scRNA-seq data (Aibar et al., 2017), we inferred ATF3 to be one of the top TFs regulating the IR-ATM signature (Fig. 8 D), including Mmp12 (Table S13). This computational inference was confirmed by ATF3 chromatin immunoprecipitation sequencing (ChIP-seq) results (Krebs et al., 2014), with a 3.4-fold difference (>2 is significant) between WT and Atf3 KO macrophages for the Mmp12 gene region in addition to 29 genes (out of 95) from the IR-ATM signature also targeted by ATF3 (Table S13). Both OV CFS and TLR2 ligands induced the expression of Atf3 in macrophages (Fig. 8 E). Although ATF3 dampens the classic inflammatory response (M1) in macrophages (De Nardo et al., 2014; Gilchrist et al., 2006), by treating Atf3 KO and control cell lines with OV CFS, we demonstrated that ATF3 is required for the complete induction of Mmp12 expression by macrophages (Fig. 8 F, left panel), despite its inhibitory effect on Tnf and Il6 expression (Fig. 8 F, right panel).
Using a robust multiorgan network and transkingdom causal inference analysis, we were able to infer, in mice and humans, and experimentally confirm novel pathways linking the deleterious effect of HFHS on glucose metabolism and diabetes via gut microbiome alterations. HFHS expands the gut commensal O. valericigenes and a few other bacterial species that, via the release of TLR2 ligands, may induce Mmp12 expression in adipose tissue macrophages in an ATF3-dependent manner. Mmp12 is part of a meta-signature expressed in a subset of WAT macrophages (IR-ATM) associated with T2D in mice and humans. MMP12 is a negative regulator of glucose metabolism that, through the damage of antioxidative and mitochondrial gene expression, causes OXPHOS dysfunction and IR in adipocytes, leading to systemic dysfunction of glucose metabolism in diabetes.
Dysregulation in the OXPHOS pathway in different organs, including WAT, is an important factor contributing to the pathogenesis of T2D (Lee et al., 2018; Rocha et al., 2020). Our study not only pointed to a central role of WAT mitochondrial OXPHOS damage in T2D, but more importantly showed for the first time the importance of gut microbiota in this process. We characterized a gene expression signature of microbiota-dependent MMP12-expressing macrophages (IR-ATM) that is associated with IR in humans. Notably, among six different clinical cohorts (one single cell and five whole tissue transcriptomes), where IR-ATM enrichment was observed, four cohorts represented patients with a similar level of obesity, differing only in IR. Our results agree with the association studies reporting that expression of MMP12 in WAT positively correlates with IR in humans (Lee et al., 2018) and another study that showed MMP12 plasma levels were elevated in patients with T2D and were associated with atherosclerosis (Goncalves et al., 2015).
Previous studies showed the critical role of diverse adipose macrophages, classified into several subtypes beyond classical M1/M2, in animal models of T2D (Russo and Lumeng, 2018). These macrophages can be identified by individual markers (e.g., LGALS3 [Li et al., 2016; Pejnovic et al., 2013], CD36 [Silva et al., 2019], and CD9 [Hill et al., 2018]) and have a negative impact on metabolism. Our analysis of multiple scRNA-seq datasets established that the previously identified subtypes converge into a single subset of metabolically harmful adipose macrophages (IR-ATM) characterized by a gene signature that includes both previously identified and novel marker genes. Moreover, we found that ATF3 (an inhibitor of proinflammatory cytokines such as TNF and IL-6) is a key TF for IR-ATMs and a positive regulator of Mmp12 expression. This is an additional strong indicator that a special type of macrophage (IR-ATMs), rather than classical inflammatory M1 (Chylikova et al., 2018; Fujisaka, 2021), promotes IR in WAT.
The negative effect of MMP12 on OXPHOS and consequently glucose uptake by adipocytes is a novel mechanism of IR, which could not be predicted based on the existing knowledge about MMP12. Thus, our study adds a new chapter to the variety of roles played by this protein ranging from a pathological factor in COPD (Hunninghake et al., 2009) to an antiviral response molecule (Marchant et al., 2014).
Our observation that TLR2 and MYD88 are involved in the control of Mmp12 is in line with reports demonstrating a role of TLR2 in T2D (Ehses et al., 2010) and of MYD88 in HFHS-induced inflammation in WAT (Tran et al., 2020). Accordingly, our findings suggest that microbial regulators of Mmp12 expression are not limited to one pathobiont (i.e., OV). Importantly, OV’s effect on the IR-ATM signature in adipose tissue in mice shows high similarity with gene expression alterations in humans, suggesting that microbiota-modifying treatments targeting IR-ATMs may offer a therapeutic strategy for T2D. Furthermore, our analysis predicted that at least two other bacteria (Barnesiella viscericola and Hydrogenoanaerobacterium saccharovorans) may impair systemic metabolism by inducing Mmp12 in WAT. Notably, these three bacteria (OV, B. viscericola, and H. saccharovorans) and/or their parent genera were negatively associated with metabolic health in some studies (Anhe et al., 2017; Jung et al., 2016; Wang et al., 2015). Our study, however, is the first to pinpoint a mechanistic role for these bacteria in T2D as well as to propose a molecular mechanism of their potential action. Although our study identified the Mmp12 pathway to be instrumental in mediating effects of HFHS on IR and T2D, alternative pathways are expected to exist, and possibly, to become prevalent in the presence of a different gut microbiota (Leshem et al., 2020). For example, among bacteria associated with a HFHS in the current data, OV acts through the induction of Mmp12 in macrophages, while we previously described that R. ilealis worsens glucose tolerance by inhibiting insulin levels, which may be relevant to more advanced stages of disease (Thévenod, 2008; Rodrigues et al., 2021).
In conclusion, our study developed a multiorgan model corresponding to initial stages of T2D and revealed a new mechanism of interaction between HFHS and microbiota in the pathogenesis of this illness. Specifically, we demonstrated that certain members of the gut microbial community (e.g., pathobionts such as O. valericigenes) exacerbate the detrimental effects of a HFHS on mitochondrial dysregulation in WAT by promoting IR-ATMs and their effector molecule MMP12. Our work sheds light on the one of most enigmatic phenomena in the field of T2D, potentially explaining why some patients with obesity develop diabetes while others do not.
Materials and methods
Animals and diets
7-wk-old SPF C57BL/6 male mice were purchased from The Jackson Laboratory and housed in a controlled environment (12-h daylight cycle) in groups of four to five per cage at the Laboratory Animal Resources Center at Oregon State University, with ad libitum access to food and water. After 1 wk of acclimatization, mice were either switched to a WD designated in this study as HFHS diet D12451 (40 kcal % lard and 17 kcal % sucrose) or to a matched ND D12450K (ND; 4 kcal % lard and 0% sucrose) produced by Research Diets (Table S14). For initial experiments (shown in Fig. 1), diets were γ-irradiated at 50 kGy at the Oregon State University Radiation Center. In follow-up experiments (shown in Figs. 2 and 4), diets irradiated by the manufacturer were used.
Mice were on these diets for 8 wk. C57BL/6 GF were bred at the Oregon State University Germ-Free Mouse Core at the Laboratory Animal Resources Center. At 8 wk of age, male mice were transferred to ventilated racks (Innovive) to fully sealed disposable cages with HEPA-filtered air supply and fed either irradiated ND or HFHS for 8 wk. Each group had five to six mice per experiment and repeated twice.
Mmp12 knockout (B6.129X-Mmp12tm1Sds/J) and C57BL/6J controls initially purchased from The Jackson Laboratory were bred to generate heterozygous Mmp12 KO mice. The heterozygous mice were then interbred to generate littermate KO Mmp12−/− and WT Mmp12+/+ mice. Genotyping was proceeded according to the protocol provided by The Jackson Laboratory with modifications. The following primers were used: mutant forward, 5′-GCTACTTCCATTTGTCACGTCC-3′; WT forward, 5′-ACTGGGCAACTGGACAACTC-3′; common, 5′-CTCCCATGTGCTGGGATTAC-3′. Briefly, a small piece of ear tissue was suspended in 50 μl lysis buffer containing 0.2% SDS, 100 mM Tris-HCl, 100 mM NaCl, 10 mM EDTA (pH = 8), and 1 μl proteinase K (Ambion, 20 mg/ml), and then digested at 55°C for 2 h. After being diluted with 450 μl distilled water, vortexed, and spun down, the supernatant was used for PCR (GoTaq Master Mixes) with the following program: initial denaturation at 94°C for 3 min, 10 cycles of 94°C for 30 s, 65°C (touch-down, 0.5°C/cycle) for 30 s, and 68°C for 30 s, 28 cycles of 94°C for 30 s, 60°C for 30 s, and 72°C for 30 s, and a final extension at 72°C for 2 min, hold at 4°C.
Adult male littermates (Mmp12 KO and littermate WT) were fed a HFHS for 5–6 wk before conducting the glucose tolerance test (GTT) and insulin tolerance tests. All mice were housed individually. The experiment was performed twice.
Work on mice was carried out in accordance with the Guide for the Care and Use of Laboratory Animals as adopted and promulgated by the U.S. National Institutes of Health and was approved by Oregon State University’s Animal Care and Use Committee (protocol IACUC-2021-0202).
Mmp12 inhibitor experiments
Adult male C57BL/6J SPF and GF mice were used. Mice were first fed a HFHS for 4 wk to induce glucose intolerance then divided into two groups. Mice were injected intraperitoneally with an MMP12 selective inhibitor MMP408 (3.5 mg/kg, dissolved in PBS containing 10% DMSO) or vehicle every day for 2–3 wk while continuing HFHS feeding. All mice were housed individually. The experiment was performed twice.
Microbe supplementation experiment
8-wk-old male C57BL/6J mice were fed with an ND and gavaged with 109 CFU of O. valericigenes DSM 18026 or R. ilealis DSM 25109 in PBS for 4–5 wk every day.
Mice were fasted for 6 h during the light phase with free access to water. A concentration of 2 or 1 mg kg−1 glucose (Sigma-Aldrich) was injected intraperitoneally. Blood glucose was measured at 0 (immediately before glucose injection), 15, 30, 60, and 120 min with a Freestyle Lite glucometer (Abbot Diabetes Care).
Insulin tolerance testing
Mice were fasted for 6 h during the light phase with free access to water. Insulin (Novolin R) was injected intraperitoneally using the dose of 1 U/kg body weight. Blood glucose was measured at 0 (immediately before insulin injection), 15, 30, 60, and 120 min with a Freestyle Lite glucometer (Abbot Diabetes Care).
Fasting insulin and fasting glucose
Mice were fasted for 6 h with free access to water. Fasting blood was collected either via submandibular bleed or from the tail vein. Insulin and glucose levels in fasting plasma or serum was measured with Mouse Insulin ELISA Kit (Crystal Chem) and Glucose Colorometric Assay Kit (Cayman Chemical), respectively, according to the manufacturer’s protocol. Homeostatic model assessment for insulin resistance (HOMA-IR) was calculated according to formula: HOMA-IR = glucose (mg/dl) × insulin (μU/ml)/405.
Antibiotic administration experiments
SPF mice were fed ND or HFHS for 4 wk, followed by the administration of antibiotics in drinking water for 4 wk. The following doses of antibiotic were used: ampicillin (Amp), 1 g/liter; neomycin, 1 g/liter; metronidazole (Met), 1 g/liter; vancomycin, 0.5 g/liter; norfloxacin, 1 g/liter; and cefoperazone (Cef) 0.5 g/liter. Body weight and feed intake were measured each week. GTT, insulin, and fasting glucose was measured every 4 wk. At the end of 8 wk of the experiment (i.e., 4 wk on antibiotics), Mmp12 gene expression was measured in WAT using quantitative PCR (qPCR) and stool microbiota was sequenced using 16S rRNA sequencing. Two independent experiments were conducted with n = 5 each.
Fecal transplantation experiments
GF mice were fed either irradiated ND or HFHS for 4 wk and then gavaged with 250 μl of inoculum prepared either from fecal contents from ND- or HFHS-fed mice. Gavaging was done on day 1 and day 3. The change in body weight, feed intake, glucose tolerance, and insulin were measured on day 4, 9, 16, and 30 after first gavage.
For humanized microbiota mice experiments, 7–8 wk old, male, GF mice were gavaged with 100 μl inoculum prepared from stool of a healthy human on day 1 and day 3 and fed with either irradiated ND or HFHS for 8 wk. Systemic measurements were done on week 3 and 8.
At the end of each experiment, mice were euthanized by cervical dislocation. The right lobe of the liver, terminal ileum, epididymal fat pad, subcutaneous fat, and gastrocnemius muscle were collected and snap frozen at −80°C until future use.
Protein carbonyl assay
Around 80 mg in epididymal adipose tissue was homogenated in 600 μl radioimmunoprecipitation assay buffer containing protease inhibitors using an OMNI Bead Ruptor and 2.8 mm ceramic beads (OMNI International). Protein carbonyl levels were measured with Protein Carbonyl Content Assay Kit (Sigma-Aldrich) according to the manufacturer’s protocol. The final values were normalized by total protein amount measure by the bicinchoninic acid assay.
Preparation of inoculum for fecal transfer
L-cystein at 0.05% solution was prepared in PBS in sterile condition. Stool from ND or HFHS fed mice was mixed with L-cystein solution to final concentration of 10% g/ml and vortexed vigorously for 5 min. This solution was centrifuged at 200 g for 5 min and aliquots from supernatant was prepared and stored at −80°C until gavaging. Inoculum from human sample was prepared in the same way.
O. valericigenes DSM 18026 was cultured in peptone yeast glucose (PYG) medium (AS 825; Anaerobe Systems) medium for 2–3 d at 30°C in an anaerobic GasPak jars (BD), then harvested, and stored in PBS containing 15% glycerol at −80°C. The R. ilealis DSM 25109 culture was provided by DSMZ in carboxymethyl cellulose medium after culture for 24 h at 37°C.
To prepare CFS from O. valericigenes, the bacterium was grown in PYG medium for 3 d, culture was centrifuged at 4,000 rpm for 20 min at 4°C, supernatant was filtered through 0.2 μM low protein binding filter and aliquoted in 2 ml tubes, then stored at −80°C. To prepare the HK O. valericigenes, the bacterial suspension was washed twice with broth, heated at 70°C for 10 min, aliquoted, and stored at −20°C. To confirm that there were no viable bacteria, one aliquot was thawed and plated in PYG agar plate, then incubated anaerobically in 30°C for 3 d.
3T3-L1 cells were grown in DMEM (high glucose) media supplemented with 10% FBS in a 5% CO2 incubator at 37°C. For differentiation, cells were seeded in a 96-well plate at a density of 5,000 cells/well. After 24 h, the medium was shifted to differentiation medium (90% DMEM, 10% FBS, 1.0 µM dexamethasone, 0.5 mM methylisobutylxanthine, 1.0 µg/ml insulin). After 48 h, the medium was removed and the cells were cultured in adipocyte maintenance medium (90% DMEM, 10% FBS, 1.0 µg/ml insulin) for another 48 h. Then the medium was changed to regular growth medium (90% DMEM, 10% FBS) and the cells were treated with activated rMmp12 or vehicle for 24 h. Afterward, the cells were incubated in glucose uptake testing medium (low glucose DMEM containing 0.2% BSA) with or without 100 nM insulin for 12 h. Then, 3 μl medium from each well was used for glucose assay using a glucose assay kit (Cayman). Levels of glucose uptake were calculated by subtracting the glucose readout from that of the fresh medium. Cells were frozen in RLT buffer and stored in −80°C for RNA extraction.
RAW cells were grown in DMEM (high glucose) media (Thermo Fisher Scientific) supplemented with 10% FBS in 5% CO2 at 37°C. Cells were than seeded in a 24-well plate with 250,000 cells per well in 500 μl complete media for 24 h. Cells were then treated with either 5% or 10% of CFS or HK for 24 h. For the control, respective concentration of PYG media was added to separate wells. Following 24 h treatment, media was removed, cells were scraped in RLT buffer (Qiagen), and stored at −80°C until RNA extraction.
For the siRNA knock down experiments, Myd88, Ticam-1, Ripk2, Mavs siRNAs, and scrambled controls were purchased from Thermo Fisher Scientific. Then 100,000 RAW cells per well in 500 μl complete medium were seeded in 24-well plate for 24 h. After 24 h, medium was replaced with fresh medium and cells were treated with siRNAs or scrambled control for 24 h. The final concentration of each siRNA was 20 pmol. After 24 h of siRNA transfection, 5% CFS or its control was added to respective wells and incubated for 24 h, followed by collection of cell lysates in RLT buffer.
WT, Ticam-1ko, and Myd88ko immortalized mouse macrophages (IMM) were originally generated in the Katherine A. Fitzgerald lab and grown in DMEM high glucose media supplemented with 10% FBS and 20 mM HEPES. Then 100,000 cells were seeded in 24-well plate in 500 μl complete medium for 24 h. All three types were seeded in the same plate. After 24 h of seeding, cells were treated with either 5% CFS or control for 24 h, then collected in RLT buffer and stored at −80°C.
TLR agonist experiments
WT IMM were originally generated in the Katherine A. Fitzgerald lab and grown in DMEM media (4.5 g/liter glucose) supplemented with 10% FBS, 1% antibiotics, and 20 mM HEPES. A total of 30,000 cells/well were seeded in triplicate in 96-well plates for 18–24 h in complete medium and incubated at 37°C with 5% CO2. Cells were then stimulated for 6 h with either 1% vol/vol O. valericigenes CFS or PYG media as control, TLR2 agonist (HKLM at 10E+8 cells/ml; InvivoGen), or TLR5 agonist (FLA-ST Ultrapure 1 µg/ml; InvivoGen). Cells were then collected in RLT lysis buffer by scraping before being stored at −80°C until RNA extraction.
Atf3 KO experiments
WT and Atf3ko IMM were originally generated by Dominic de Nardo (De Nardo et al., 2014) and grown in DMEM media (4.5 g/liter glucose) supplemented with 10% FBS, 1% antibiotics, and 20 mM HEPES. Then 30,000 cells/well were seeded in triplicate in 96-well plates for 18–24 h in complete medium and incubated at 37°C with 5% CO2. Cells were then stimulated for 6 h with 1% vol/vol O. valericigenes CFS or PYG media as control (O. valericigenes growth media), then collected in RLT lysis buffer by scraping before being stored at −80°C until RNA extraction.
RNA extraction and qPCR
RNA was extracted using an OMNI Bead Ruptor and 2.8-mm ceramic beads (OMNI International) followed by Qiashredder and an RNeasy kit (the right lobe of liver, terminal ileum, and cells), RNeasy Fibrous Tissue Mini Kit (for gastrocnemius muscle), and an RNeasy Lipid Tissue Mini Kit (for epididymal fat) using Qiacube (Qiagen) automated extraction according to the manufacturer’s specifications with DNase treatment (RNase-Free DNase Set; Qiagen). Total RNA was quantified using Quant-iT RNA Assay Kit (Thermo Fisher Scientific). Complementary DNA was prepared using qScript reverse transcription kit (QuantaBio). qPCR was performed using the Perfecta SYBR Green Fastmix mix (Quantabio) and the StepOne Plus Real Time PCR system and software (Applied Biosystems).
The primers used are listed in Table 1.
Normalization of qPCR data and statistics
Raw cycle threshold (CT) values of genes of interest from qPCR runs were normalized to CT values of the housekeeping gene (Polr2c) via delta CT method before calculation of relative expression using the 2−∆CT method. Relative expression levels of genes were median-normalized and/or also log2 transformed before departure from normality was assessed with a Shapiro–Wilk test. If data were normally distributed, one- or two-tailed parametric t tests were used, otherwise a Mann–Whitney test was used to compare treatments. Individual plot legends note how data are visualized specifically.
Sequencing of RNA (RNA-seq)
RNA libraries were prepared with the QuantSeq 3′mRNA-Seq Library Prep Kit (Lexogen) for the Apollo 324 NGS Library Prep System and sequenced using Illumina NextSeq. Sequences were processed to remove the adapter, polyA, and low-quality bases by BBTools (https://jgi.doe.gov/data-and-tools/bbtools/) using bbduk parameters of k = 13, ktrim = r, forcetrimleft = 12, useshortkmers = t, mink = 5, qtrim = r, trimq = 15, minlength = 20. Reads were aligned to mouse genome and transcriptome (ENSEMBL NCBIM37) using Tophat (v. 2.1.1) 70 with default parameters. Number of reads per million for mouse genes were counted using HTSeq (v. 0.6.0) 71 and quantile normalized. BRB-ArrayTools was used to identify DEGs.
Adipose tissue samples (n = 31) were blocked and randomized to mitigate any potential systematic experimental biases (Oberg and Vitek, 2009). The tissue samples were homogenized in cell lysis buffer (20 mM Tris-HCl pH 7.5, 150 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% vol/vol Triton X-100, 2.5 mM sodium pyrophosphate, 1 mM β-glycerophosphate, 1 mM Na3VO4, 1 μg/ml leupeptin, and 1 mM PMSF [added immediately before use]). Homogenization was achieved using a Bioruptor Plus sonicator and Bioruptor Protein Extraction Beads (Diagenode, Inc.) at 4°C. The samples were microcentrifuged (20 min; 4°C; 20,000 g) and the supernatants were retained. Bicinchoninic acid protein concentration assays were performed using a Micro BCA Protein Assay Kit (Thermo Fisher Scientific).
Multiplexed Western blots were performed using the following primary antibodies from Cell Signaling Technology: total AKT (C67E7; cat# 4691S), AKT phospho-Ser473 (D9E; cat# 4060S), total ERK1/2 (L34F12; cat# 4696S), ERK1/2 phospho-Thr202/Tyr204 (D13.14.4E; cat# 4370S), total p38 MAPK (D13E1; cat# 8690S), p38 MAPK phospho-Thr180/Tyr182 (D3F9; cat# 4511S), total PTEN (D4.3; cat# 9188S), PTEN phospho-Ser380/Thr382/Thr383 (cat# 9554S), total GSK-3β (D5C5Z; cat# 12456S), total JNK (cat# 9252S), and total cofilin (D3F9; cat# 5175S). For each tissue, a sample pool was prepared and probed using each primary antibody individually, and these results were used to ensure that the analytes would not be misidentified due to multiplexing. Between one and three primary antibodies were used per Western blot. Each Western blot was used to analyze either adipose or liver samples (not both). At least two lanes were used for a MW ladder standard and at least two lanes were used for the pooled tissue sample. Imaging was performed using a ChemiDoc MP imager, and densitometry was performed using Image Lab v. 6.0.1 (both Bio-Rad, Inc.).
DNA extraction and 16S rRNA gene libraries preparation
For microbial measurements, stool pellets were collected at T1 (4 wk) and T2 (8 wk). Ileum contents were collected at T2 (8 wk) only. To get microbial DNA, frozen fecal pellets were resuspended in 1.4 ml ASL buffer (Qiagen) and homogenized with 2.8 mm ceramic beads followed by 0.5-mm glass beads using an OMNI Bead Ruptor (OMNI International). DNA was extracted from the entire resulting suspension using a QiaAmp mini stool kit (Qiagen) according to the manufacturer’s protocol. DNA was quantified using a Qubit broad range DNA assay (Life Technologies). The V4 region of the 16S rRNA gene was amplified using universal primers (515f and 806r) as previously described. Individual samples were barcoded, pooled to construct the sequencing library, and then sequenced using an Illumina Miseq (Illumina) to generate pair-ended 250 bp reads.
16S rRNA gene sequencing data analysis
The samples were demultiplexed, and forward-end fastq files were analyzed using QIIME v. 1.9.1. The default quality filter parameters from QIIME’s split_libraries_fastq.py were applied to retain high quality reads (Phred quality score ≥20 and minimum read length = 75% of 250 nucleotides). To cluster 16S rRNA gene sequence reads into operational taxonomic units (OTUs) and assign taxonomy, a closed reference OTU picking with 97% sequence similarity was performed using UCLUST and the Greengenes reference database v. 13.8. The reference sequence of candidate OTUs from the Greengenes database was used to obtain species level taxonomic assignment using Megablast (top hit using default parameters). A threshold of 99% cumulative abundance across all samples in an experiment was used to retain abundant microbes, thus removing OTUs with approximately 1% abundance across all samples in that experiment. The read counts were normalized using cumulative sum scaling, accounted for bacterial DNA quantity based on qPCR, followed by quantile normalization. The principal component analysis for the 16S sequencing data was created using Clustvis and the GraphPad Prism software.
Isolation of SVF
Epididymal adipose tissue was minced in HBSS (Invitrogen) containing calcium, magnesium, and 0.5% BSA and 3 mg/ml collagenase type 1 (Rockland). After incubation at 37°C for 1 h, the cell suspension was filtered through 100-μm strainer and then washed with 20 ml wash buffer (DPBS containing 1 mM EDTA and 0.5% BSA). After centrifugation at 700 g for 10 min, the supernatant was removed and the cell pellet was re-suspended in 5 ml of RBC lysis buffer (Tonbo Biosciences) and incubated for 10 min at room temperature. A total of 10 ml wash buffer was added and the cell suspension was filtered through 40 μm and centrifuged again. CD45+ cells from the SVF were isolated using CD45 MicroBeads (Miltenyi Biotec) according to the manufacturer’s instruction.
The isolated cells were re-suspended in 1 ml FBS containing 10% DMSO and stored at liquid nitrogen tank until use.
The genes whose expressions are significantly changed by HFHS (false discovery rate [FDR] < 15%) in the conventional SPF mice were grouped into two main categories: category 1 (microbiota-dependent) and category 2 (microbiota-independent genes). Category 1A are genes whose expressions are significantly changed by the interaction effect of diet and microbiota (two-way ANOVA interaction effect FDR < 15%). Category 1B are genes whose expressions are not affected by the interaction effect (FDR > 15%) but are significantly regulated by microbiota under either ND or HFHS conditions (SPF vs. GF FDR < 15%). Genes whose expressions are not detectable in GF mice are grouped in Category 1 and cannot be categorized further. Category 2 are genes whose expression are independent of microbiota (interaction FDR > 15%, SPF vs. GF FDR > 15%).
Reconstructing the multiorgan network
First, Spearman rank correlations were calculated between all pairs of genes from four tissues (muscle, ileum, liver, and fat) and metabolic parameters (phenotypes) by pooling the samples per diet (HFHS, ND) from both experiments. Meta-analysis was performed to retain edges with the same sign of correlation coefficient in both the diets. Edges were further filtered by the following criteria: individual P value of correlation within each diet from pooled experiments <30%, combined Fisher’s P value over diets from pooled experiments <5%, and FDR cutoff of 5% for edges within tissues and 10% for phenotypes and between tissues, and edges needed to satisfy principles of causality (i.e., satisfied fold change relationship between the two partners in the HFHS vs. ND comparison). Next, correlations were calculated per diet per experiment. Finally, the edges obtained from pooling were retained if they had the same sign of correlation coefficient as in two groups (two experiments × two diets). False positive edges due to pooling experiments were removed in the creation of the network. The proportion of genes from each tissue type that made it to the final network (following statistical cutoffs) was determined by dividing the number of DEGs in each tissue after applying correlation cutoffs by the number of DEGs in each tissue prior to applying correlation cutoffs.
Detecting subnetworks essential for regulation of systemic parameters of metabolic disease
The networks were visualized in the Cytoscape Software 2.6.3. To identify subnetworks of correlated genes in fat and liver, we used the MCODE v. 1.2 (Molecular Complex Detection) plug-in for Cytoscape to identify clusters (subnetworks) of correlated genes. The largest two subnetworks in each tissue were selected for further analyses, as all other subnetworks were representing parts of them.
Finding the shortest paths between subnetworks
Dijkstra’s algorithm was used in the python module NetworkX v. 2.2 to calculate all pairwise shortest paths between the nodes of each subnetwork–subnetwork pair and each subnetwork–phenotype pair. The weighted average (reciprocal) of the shortest paths was calculated to determine the distance of each subnetwork from systemic parameters of metabolic disease. A chi-square test was then conducted on the shortest path distributions for each pair of subnetworks.
Identifying important nodes for the information flow between subnetworks
BiBC has previously been used to find nodes that control the information flow between two parts of a network (Morgun et al., 2015). Node BiBC was calculated between each subnetwork–phenotype pair of the multiorgan network. Per pair, we counted the number of nodes in the top 20% (ranked in decreasing order of BiBC) belonging to each of the identified subnetworks. A chi-square test was used to compare the distributions of the number of top BiBC nodes belonging to each subnetwork for each subnetwork–phenotype pair.
Establishing the information flow within network by incorporating the results of the shortest paths and BiBC analyses
To combine the results from the previous two analyses, we created a model that demonstrates how close each subnetwork is to systemic parameters of metabolic disease as well as how much one of the identified subnetworks relies on another subnetwork to influence those systemic parameters. Both the line weights and distance from systemic parameters indicate the average shortest path between the subnetwork and systemic parameters. Arrow weights indicate the number of nodes in the top 20% BiBC of the network indicated by the arrow direction (i.e., the OXPHOS/mitochondria subnetwork has more nodes in the top 20% BiBC than the inflammation subnetwork when using fatty acid metabolism and systemic parameters as the two groups of the BiBC analysis).
Analysis of association of IS proteins and gene expression
Spearman correlations between the expression levels of DEGs (ND vs. HFHS, FDR < 15%), and p-Akt/Akt protein levels in WAT were calculated using the following cutoffs: (1) correlation sign (±) is consistent between two experiments, (2) when the two experiments are pooled, the correlation P value is <0.05 (two-tailed). Spearman correlations between the expression levels of DEGs and p-p38/p38 protein levels in WAT were also calculated (two-tailed P < 0.05). To compare the relationship between the sign of the correlation (the levels of gene expression vs. the levels of IS protein) with the regulation of HFHS on this gene (up/down), a 2 × 2 contingency table was made and chi-square test was performed. Next, ToppGene was used to perform gene enrichment analysis on the two groups of genes: (1) upregulated by HFHS and negatively correlated with IS proteins; (2) downregulated by HFHS and positively correlated with IS proteins.
Analysis of IR signature in WAT
A gene expression meta-signature of IR in adipose tissue obtained from Jung et al. (2018) was used for this analysis. Dimension reduction using median summary metric was carried out to compare the presence of the signature in different groups using RNA-seq data. For this, we first normalized each gene expression using the median in all groups. For the genes downregulated in T2D, the reciprocal values were used. Then the sum of all genes in the signature was used to indicate the level of presence of the signature in each sample.
Analysis of the relationship between the effects of HFHS and insulin on gene expression
To analyze the relationship between the effects of HFHS and that of insulin, global gene expression was profiled by RNA-seq of the WAT of mice fed with ND or HFHS (n = 10) and 3T3-L1 cells treated with insulin or vehicle (n = 3), respectively. Spearman correlations were conducted between fold change of gene expression in mice (HFHS/ND, FDR < 15%) and fold change of gene expression in 3T3-L1 cells (insulin/vehicle) using GraphPad Prism7.
The raw gene expression matrix (UMI counts per gene per cell) was filtered, normalized, and clustered using R (https://www.R-project.org/). Cell and gene filtering were performed as follows: cells with a very small library size (<2,500) and a very high (>0.5) mitochondrial genome transcript ratio were removed. Genes detected (UMI count > 0) in less than three cells were removed. After log normalization, clustering was performed using standard Seurat package procedures. Principal component analysis was used to reduce the number of dimensions representing each cell. The number of components used for analysis was determined based on the elbow of a scree plot. Selection of a biologically relevant number of clusters was based on differential expression between neighboring clusters. Differential expression between clusters was calculated using a likelihood-ratio test for single-cell gene expression implemented in Seurat at a family-wise error rate of 5%. Neighboring clusters in principal component space were identified as the next-nearest cluster to each cell after the cell’s assigned cluster. Clusters were visualized using t-distributed stochastic neighbor embedding of the principal components (spectral t-SNE), as implemented in Seurat. The cell-type identities for each cluster were determined manually using a compiled panel of available known immune cell marker expression.
Identifying Mmp12-positive cells
The detectable expression of Mmp12 was identified in subset of cells in each of the single cell RNA-seq data set for further analysis.
The gene expression was correlated with the expression of Mmp12 in these cells in each dataset using the Pearson correlation coefficient and a t test was performed to determine the significance. The Bonferroni-Hochberg multiple comparisons procedure was used to compute the adjusted P value. Genes with mean Spearman correlation coefficients above 0.2 and a P value of 0.001 or less was selected for meta-analysis.
Mmp12 cluster cells expression was analyzed for differential expression across other clusters. The DEGs whose fold change is >2 in at least three or more datasets with P value <0.001 were selected. Significant genes were further shortlisted based on Bonferroni-Hochberg corrections.
The Mmp12 meta-signature list of genes was selected based on the meta-analysis of overlapped gene sets from both Mmp12-positive cell meta-correlations and Mmp12 cluster meta-differential lists of genes.
Dimension reduction of Mmp12 meta-signature using median summary metric was carried out to compare the metabolic macrophage cell types from different publicly available studies with data from the Gene Expression Omnibus (GEO; GSE126407, GSE128518 and GSE113595). For this, we normalized each gene expression and calculated the median value of all genes for each sample, obtaining a single value for each macrophage cell type in each of the study.
Human diabetic patient adipose tissue SVF scRNA-seq analysis
scRNA-seq profiling of SVF from subcutaneous and visceral adipose tissue of obese individuals with or without T2D available in GEO as Discovery cohort (GSE129363) were analyzed to check the relevance of mouse Mmp12 meta-signature and signature expressing macrophages to human diabetes. The analysis was performed using Seurat package as described previously.
The gene average expression of T2D and non-T2D patients’ samples were calculated separately. The percentage of genes in each condition for Mmp12 meta-signature was calculated and observed to be significantly higher in T2D condition. The single cell gene expression data for cluster specific markers were used to identify and classify the clusters. The cluster enriched in Mmp12 meta-signature identified using the expression of marker genes as part of myeloid cells were then labeled MC1 cluster. Percentage distribution of SVF of MC cluster cells as compared with other clusters in these single cell samples were estimated to show the enrichment of MC1 cells as compared with other clusters in T2D samples. This was also evident in the dot plot of average expression and percentage of cells with Mmp12 meta-signature enrichment in MC1 cluster cells among all clusters of the adipose tissue SVF. Additionally, for the specific cell type identification of MC1 cluster as macrophages, scQuery, a single cell database was queried with highly differential and are also part of Mmp12 meta-signature genes in MC1 cluster.
Analysis of TFs regulating IR-ATMs signature
The gene regulatory network (GRN) was inferred from 1,758 obese adipose tissue macrophages scRNA-seq data (GSE117176/GSM3272967) using standard SCENIC (Aibar et al., 2017) workflow. The normalized enrichment score for the GRN was estimated and 38 TFs were identified in this data. Atf3 was the only TF, with Mmp12 as part of predicted regulon with the normalized enrichment score. There were additionally 13 genes from the Mmp12 meta-gene signature identified to be part of Atf3 regulon. A single gene from the signature was predicted to be part of Maf regulon in this GRN. Expression matrix units: UMI counts.
Meta-analyses of Mmp12 meta-signature in human studies
GEO was searched for studies that compared the expression of one group (either obese or obese with T2D) to another group (either non-obese or obese with T2D, respectively). Studies were further filtered out if (1) either body mass index or fasting glucose were not listed, (2) the T2D group consisted of non-obese individuals, or (3) individuals in the study were not over the age of 18 (adults). After filtering out, five studies remained: GSE133786 (Yang et al., 2019), GSE24883 (Klimcakova et al., 2011), GSE20950 (Hardy et al., 2011), GSE29231 (Jain et al., 2013), and GSE26637 (Soronen et al., 2012).
Using BRB-ArrayTools (https://brb.nci.nih.gov/BRB-ArrayTools/), the gene expression data were then log2-transformed and quantile-normalized, then subset for the 87 genes in the human Mmp12 signature. Since these studies were all microarray data, duplicate arrays were managed by selecting the arrays in the disease group (the group with either obese insulin-resistant individuals or obese noninsulin-resistant individuals) that had the maximum geometric mean. Meta-Mar v2.7.0 (http://www.meta-mar.com/) was used to perform meta-analysis, while GraphPad Prism 8.2.1 was used to create the volcano plots. The online tool Morpheus (https://software.broadinstitute.org/morpheus) was used to create the heatmaps of genes.
Transkingdom causal inference analysis: Prediction of microbes that induce Mmp12 gene expression and worsen glucose metabolism (steps refer to Fig. S4 A)
Step 1. SPF mice were treated with six antibiotics, systemic metabolic parameters were measured, stool microbiota was sequenced, and Mmp12 gene expression in WAT was measured using qPCR at the end of 8 wk. Outliers were identified using Grubbs (α = 0.05). Statistical significance was calculated by comparing the ND and antibiotics groups against HFHS.
Step 2. Microbes were selected if they consistently showed the same (antibiotics vs. HFHS) fold-change direction in the three antibiotics group where Mmp12 was changed (Amp, Met, Cef) in both experiments and were significantly changed (Fisher’s P value <0.1) by at least two of the three antibiotics.
Step 3. OTUs passing the above criteria were correlated (Spearman correlation) with systemic measurements, and Mmp12 across all the samples treated with antibiotics. Meta-analysis was performed to retain correlations with same sign of coefficient in both experiments and accounting for correlation inequalities (Yambartsev et al., 2016; i.e., satisfied fold change relationship between the two partners in the antibiotics vs. HFHS comparison).
Step 4. Next, the OTUs were correlated with genes from Mmp12/IR-ATM meta-signature and direct neighbors of Mmp12 from the multi-organ network by pooling the HFHS samples from both experiments (i.e., samples from mice that did not receive antibiotics), therefore eliminating the common cause (Pearl et al., 2009; Reichenbach et al., 1991) and accepted edges that satisfied correlation inequalities (Yambartsev et al., 2016; i.e., satisfied fold change relationship between the two partners in the HFHS vs. ND comparison). Correlations were calculated in HFHS samples per experiment. Edges obtained from pooling were retained if they had the same sign of correlation coefficient as in HFHS group per experiment. False positive edges due to pooling experiments were removed. The cumulative scores obtained by adding the (OTU-gene) correlation coefficients were used to select top candidate OTUs for experimental validation.
Online supplemental material
Fig. S1 shows additional details of validation of network inferences revealing microbiota-dependent IR in WAT. Fig. S2 shows Mmp12/IR-ATM signature and Mmp12 blockade in mice. Fig. S3 shows Mmp12/IR-ATM signature in patients with T2D. Fig. S4 shows inference and validation of microbes inducing Mmp12/IR-ATMs. Fig. S5 shows testing innate adaptor molecules used by OV to stimulate Mmp12 expression and IR-ATM signature. Table S1 shows analysis of systemic metabolic parameters in SPF and GF mice fed with HFHS and control diets. Table S2 lists DEGs between HFHS and ND in conventional (SPF) mice (FDR < 15%) in four organs and their expression in GF mice. Table S3 shows the total number and percentage of DEGs and percentage of DEGs that were retained in the networks in the WAT, liver, muscle, and ileum. Table S4 shows chi-square values of each pairwise comparison between subnetworks for the shortest paths and BiBC analyses performed in Fig. S2 E. Table S5 shows the intensity values of each individual sample (WAT) from Western blot. Table S6 shows the effects of HFHS on gene expression in WAT of SPF mice and of insulin on gene expression in adipocytes. Table S7 shows Mmp12 and CD45(Ptprc) positive cell counts for eight publicly available scRNA-seq datasets. Table S8 shows Mmp12-IR-ATM meta-signature genes. Table S9 shows gene ontology enrichment analysis of genes regulated by MMP12 inhibitor (ex vivo). Table S10 shows genes significantly regulated in WAT in MMP12 inhibitor-treated mice and their correlation to metabolic parameters (ex vivo). Table S11 shows genes significantly regulated by rMMP12 and insulin in 3T3-L1 cells (in vitro). Table S12 shows Mmp12 meta-signature gene expression in human datasets and in macrophage cell lines stimulated with OV supernatant. Table S13 shows genes from IR-ATM signature targeted by Atf3 as predicted by SCENIC and identified via ChIP-seq. Table S14 shows composition of experimental (HFHS) and control (normal) diets.
The authors would like to thank Oksana Gavrilova, Mark Rochman, Allison K. Ehrlich, Donald B. Jump, Barbara Gvakharia, and Valerian Dolja for helpful discussions and feedback; Ian Fraser (National Institute of Allergy and Infectious Diseases, National Institutes of Health, Bethesda, MD) for kindly providing immortalized macrophage cell lines; Philip Proteau (Oregon State University, Corvallis, OR) for help with reagents to activate recombinant MMP12; and the Oregon State University Center for Quantitative Life Sciences for sequencing and computational support.
This work was supported by the National Institutes of Health R01 grant DK103761 (N. Shulzhenko), an Oregon Medical Research Foundation award (A. Morgun), and intramural National Institutes of Health funding (R.R. Rodrigues, A.K. Dzutsev, A. Nita-Lazar, and G. Trinchieri).
Author contributions: Conceptualization: Z. Li, M. Gurung, R.R. Rodrigues, N.K. Newman, J. Padiadpu, A.K. Dzutsev, A. Nita-Lazar, G. Trinchieri, N. Shulzhenko, A. Morgun; data curation: Z. Li, M. Gurung, R.R. Rodrigues, N.P. Manes, R.L. Greer, S. Vasquez-Perez; formal analysis: Z. Li, M. Gurung, R.R. Rodrigues, N.K. Newman, J. Padiadpu; funding acquisition: N. Shulzhenko, A. Morgun, A. Nita-Lazar, G. Trinchieri; investigation: Z. Li, M. Gurung, R.R. Rodrigues, N.K. Newman, J. Padiadpu, N.P. Manes, R.L. Greer, S. Vasquez-Perez, H. You, K.A. Hioki, A. Fel, A.K. Dzutsev, D. De Nardo, Z. Moulton, J.W. Pederson; methodology: Z. Li, M. Gurung, R.R. Rodrigues, N. Shulzhenko, A. Morgun; project administration: N. Shulzhenko, A. Morgun, G. Trinchieri; resources: A. Nita-Lazar, G. Trinchieri, D. De Nardo, N. Shulzhenko, A. Morgun; software: R.R. Rodrigues, N.K. Newman, J. Padiadpu; supervision: A. Nita-Lazar, G. Trinchieri, N. Shulzhenko, A. Morgun; validation: Z. Li, M. Gurung, R.R. Rodrigues, N.K. Newman, N. Shulzhenko, A. Morgun; visualization: Z. Li, M. Gurung, R.R. Rodrigues, N.K. Newman, J. Padiadpu, J.W. Pederson; N. Shulzhenko, A. Morgun; writing – original draft: Z. Li, M. Gurung, R.R. Rodrigues, N.K. Newman, J. Padiadpu, N. Shulzhenko, A. Morgun; writing – review & editing: Z. Li, M. Gurung, R.R. Rodrigues, N.K. Newman, J. Padiadpu, N.P. Manes, A.K. Dzutsev, A. Nita-Lazar, G. Trinchieri, N. Shulzhenko, A. Morgun.
Z. Li, M. Gurung, R.R. Rodrigues, and J. Padiadpu contributed equally to this paper.
N.K. Newman and N.P. Manes contributed equally to this paper.
Disclosures: The authors declare no competing interests exist.