Lower brown adipose tissue activity is associated with non-alcoholic fatty liver disease but not changes in the gut microbiota.

In rodents, lower brown adipose tissue (BAT) activity is associated with greater liver steatosis and changes in the gut microbiome. However, little is known about these relationships in humans. In adults (n = 60), we assessed hepatic fat and cold-stimulated BAT activity using magnetic resonance imaging and the gut microbiota with 16S sequencing. We transplanted gnotobiotic mice with feces from humans to assess the transferability of BAT activity through the microbiota. Individuals with NAFLD (n = 29) have lower BAT activity than those without, and BAT activity is inversely related to hepatic fat content. BAT activity is not related to the characteristics of the fecal microbiota and is not transmissible through fecal transplantation to mice. Thus, low BAT activity is associated with higher hepatic fat accumulation in human adults, but this does not appear to have been mediated through the gut microbiota.

Introduction

Non-alcoholic fatty liver disease (NAFLD) is the leading cause of chronic liver disease and is an important risk factor for type 2 diabetes and cardiovascular disease. NAFLD includes hepatic steatosis and non-alcoholic fatty steatohepatitis (NASH), which may progress to end-stage liver disease and hepatocarcinoma. Although numerous pharmacotherapies have been effective for treating NAFLD in rodents, none have been approved for use in humans. The development of NAFLD occurs in part because of increased adipose tissue lipolysis and hepatic de novo lipogenesis. The gut microbiome also plays an important role in the pathogenesis of NAFLD and its progression to more severe disease.4, 5, 6, 7, 8 Despite these important observations, therapies targeting lipolysis, lipogenesis, or the gut microbiome have had limited efficacy as stand-alone therapies for reducing NASH in humans, suggesting that additional therapeutic strategies may be required.

In numerous preclinical studies in rodents, enhancing the metabolic activity of brown adipose tissue (BAT) potently reduces liver lipids.9, 10, 11, 12 In mice housed at room temperature, BAT is a major contributor to whole-body energy expenditure because of thermogenesis mediated through uncoupling protein 1 (UCP1) and other futile cycles., Recent studies have found that the gut microbiome plays an important role in regulating BAT activity and the browning of white fat15, 16, 17 in rodents. However, it is not known if these findings are also present in humans.

In humans, BAT presents mainly in the neck and supraclavicular (SCV) region.18, 19, 20 Although BAT thermogenesis is a relatively minor contributor to whole-body energy expenditure,21, 22, 23 cold exposure increases human BAT activity and the clearance of circulating glucose, non-esterified fatty acids (NEFA), and triglycerides (TG).24, 25, 26, 27, 28 These data suggest that activating BAT may have positive effects on NAFLD. Consistent with this idea, retrospective studies have noted a relationship between 18F-fluorodeoxyglucose (18F-FDG) uptake into the SCV area and NAFLD status.29, 30, 31, 32 Although interesting, it is important to note that in order to quantitate BAT activity it is essential that cold exposure be standardized between participants and to recognize that TG and not glucose are likely the primary fuel source for activated BAT., Therefore, given previous studies in rodents and the potential influence of BAT activity on metabolic health in humans, it is important to determine if enhanced BAT activity may reduce hepatic fat content and, if so, to determine whether the gut microbiome has a role in regulating this process.

In the present study, we have investigated if BAT activity is related to hepatic fat content in adult humans. Compared with previous studies, the use of standardized cold exposure and magnetic resonance imaging (MRI) allows detailed quantification of BAT activity. In addition, given previous studies in rodents, we investigated if the gut microbiota could be an intermediary factor regulating BAT activity. We find that low BAT activity is associated with increased liver fat, but this is not related to differences in the gut microbiota.

Results and discussion

Participant characteristics

This cross-sectional two-visit study was conducted in women and men 18–57 years of age. Of the 483 people that inquired about the study, because of strict inclusion criteria described in STAR Methods, 73 (n = 43 men) participants were enrolled, and 64 participants (n = 39 men) completed both study visits (Figure S1). Of the 9 participants who did not complete visit 2, 5 were identified as ineligible at visit 1 and 4 withdrew because of symptoms of claustrophobia during the MRI scan. Of the 64 participants with both visits, post-cold images for 4 participants were removed because of motion artifacts, leaving 60 participants (n = 37 men) with usable pre- and post-cold MRI scans of the SCV region. Visits 1 and 2 occurred on average 18.5 days apart (range 14–29 days). For premenopausal women, visit 2 was scheduled in relation to their menstrual cycle as described in STAR Methods. According to the American Diabetes Association guidelines for classification of glycemic status on the basis of hemoglobin A1c (HbA1c), 6 participants had prediabetes and one had diabetes (HbA1c 8.0%). Diabetes had been previously identified in this individual, but they remained eligible as they were on no medications.

Proton density fat fraction declined in SCV but not in subcutaneous adipose tissue

After a 3 h cold exposure delivered through a water-perfused suit, set at 18°C, the absolute decline in SCV proton density fat fraction (PDFF) was 3.0% ± 2.4%, and the percentage decline relative to baseline was 4.4% ± 3.9% (Table 1). In contrast, the PDFF in abdominal subcutaneous adipose tissue (SAT) remained unchanged after cold exposure (Figure S2), consistent with expectations based on studies reported by us and others in which PDFF of BAT and white adipose tissue (WAT) had different responses to cold. The mean absolute decline in SCV PDFF in our participants was slightly greater than the previously reported range of 0.4%–2.2% (reviewed in Ong et al.). This may be attributed to the more intense cold exposure (longer duration and lower temperature) and/or the application of cold to the whole body. The shivering intensity was <3.0% of maximal voluntary contraction for all participants (Table 1), consistent with previous studies using a similar standardized cold exposure protocol.,,

Table 1

Participant characteristics and study conditions of those with and without NAFLD

	All (n = 60)	Normal liver fat (n = 30)	Hepatic steatosis (n = 29)	p
Age (years)	25.9 (22.9–36.4)	23.8 (21.2–27.2)	31.02 (25.1–43.4)	<0.001†
Gender (male)	37 (62%)	19 (63%)	18 (62%)	0.920
Weight (kg)	81.9 (18.6)	70.7 (13.6)	91.7 (14.8)	<0.001∗
BMI (kg/m2)	27.9 (6.8)	22.4 (21.0–25.0)	30.4 (26.0–35.9)	<0.001†
Waist circumference (cm)	90.4 (16.9)	76.8 (72.5–83.8)	101.4 (89.8–110.2)	<0.001†
Body fat (%)a	30.6 (12.6)	23.9 (9.6)	37.4 (11.8)	<0.001∗
Hepatic fat (%)b	5.5 (4.6–7.1)	4.6 (4.0–5.1)	7.1 (6.2–8.5)	<0.001†
Pre-cold SCV PDFF (%)	72.4 (8.0)	67.9 (7.0)	76.7 (6.2)	<0.001∗
Absolute cold-induced decline in SCV PDFF (%)	3.0 (2.4)	3.6 (2.2)	2.5 (2.4)	0.082
Cold-induced percentage decline in SCV PDFF (%)	4.4 (3.9)	4.7 (2.7–8.0)	2.5 (1.2–5.7)	0.019†
VAT (cm2)c	30.3 (26.1)	8.7 (5.2,17.2)	58.6 (22.7–64.6)	<0.001†
SAT (cm2)c	99.3 (68.1)	43.2 (29.3–68.5)	150.1 (87.2–188.3)	<0.001†
Fasting plasma glucose (mmol/L)	4.7 (4.4–4.9)	4.6 (4.1–4.8)	4.8 (4.6–5.1)	0.005†
2 h glucose (mmol/L)d	4.4 (3.8–5.8)	4.4 (1.3)	5.00 (1.4)	0.299
HbA1c (%)	5.2 (4.9–5.5)	5.1 (4.9–5.3)	5.3 (4.9–5.6)	0.112
Total cholesterol (mmol/L)	4.4 (0.8)	4.1 (0.6)	4.6 (0.8)	0.005∗
Triglycerides (mmol/L)	0.8 (0.6–1.2)	0.8 (0.2)	1.3 (0.7)	<0.001∗
HDL-C (mmol/L)	1.3 (0.3)	1.4 (0.3)	1.2 (0.3)	0.052
LDL-C (mmol/L)	2.6 (0.7)	2.3 (0.6)	2.8 (0.6)	0.004∗
Non-HDL-C (mmol/L)	3.1 (0.8)	2.7 (0.6)	3.4 (0.8)	<0.001∗
Systolic blood pressure (mm Hg)	111.4 (12.1)	105.6 (10.7)	116.2 (9.9)	<0.001∗
Diastolic blood pressure (mm Hg)	74.4 (9.8)	70.7 (9.4)	77.8 (8.9)	0.004∗
AST (U/L)	18.5 (16.0–21.0)	18.0 (16.0–21.0)	19.0 (16.0–25.5)	0.326
ALT (U/L)	20.0 (15–25.8)	16.5 (13.0–21.5)	25.0 (19.0–35.5)	<0.001†
GGT (U/L)	18.5 (14.0–25.0)	14.5 (11.0–19.0)	20.0 (18.0–29.5)	0.001†
Outdoor temperature 1 h before visit 2 (°C)	7.18 (10.0)	4.3 (7.7)	10.5 (11.2)	0.015∗
Shivering intensity (% MVC)e	1.95 (1.0–2.9)	2.3 (1.3–3.5)	1.2 (0.8–2.6)	0.017†
ΔOutlet-inlet of the suit (°C)f	1.6 (0.2)	1.6 (0.2)	1.5 (0.2)	0.671

Categorical variables are presented as n (%), and continuous variables are presented as mean (SD) for normally distributed variables and median (interquartile range) for non-normally distributed variables. BMI, body mass index; SCV, supraclavicular; PDFF, proton density fat fraction; VAT, visceral adipose tissue; SAT, subcutaneous adipose tissue; HbA1c, hemoglobin A1c; HDL-C, high-density lipoprotein cholesterol; LDL-C, low-density lipoprotein cholesterol; AST, aspartate aminotransferase; ALT, alanine aminotransferase; GGT, gamma-glutamyl transferase.

p < 0.05, independent-samples t test; †p < 0.05, Mann-Whitney U test. Differences in gender were assessed using the chi-square test.

One participant was above the weight threshold for body composition analysis via dual-energy X-ray absorptiometry (DEXA).

No liver scan was performed for one participant who did not fit inside the scanner.

No abdominal volume data were available for four participants (n = 2, acquisition error; n = 1, did not fit inside the scanner; n = 1, motion issues).

Two-hour glucose was not performed in 29 participants (n = 2, fainted during baseline bloodwork; n = 1, known type 2 diabetes; n = 26, glucose solution was later recalled, so the test was not reliable).

Motion artifacts were present for one participant.

No outlet-inlet data were available for five participants (n = 4, data were not acquired; n = 1, error in the data acquisition).

Given the variability in the baseline SCV PDFF, the decline in SCV PDFF that was relative to the baseline value was used as our primary measure of BAT activity (i.e., cold-induced percentage decline in SCV PDFF). In studies in which BAT activity was measured using 18F-FDG uptake, individuals were classified as BAT positive or negative. Using our measure, those who are BAT negative would be expected to have no cold-induced decline in SCV PDFF. Thus, we classified those with a cold-induced decline in SCV PDFF within one SD of zero as having low BAT activity (n = 33). A comparison of those with low versus high BAT activity is presented in Table S1. Neither shivering nor the change in outlet-inlet temperature were different between BAT groups (Table S1), suggesting that both groups were exposed to similar cold intensity., Consistent with previous literature, those with low BAT activity were older, were more likely to be female, and had higher total body fat. No differences in macronutrient intake or dietary fiber were identified over the 24 h prior to measurement of BAT activity (Table S2).

Characteristics of those who were NAFLD positive

The median hepatic fat fraction in the 60 participants was 5.5% (4.6%–7.1%). Hepatic fat exceeded 5.6%, the classification for hepatic steatosis (NAFLD positive), in 29 participants (48%) (Table 1). Those who were NAFLD positive were older and had higher body mass index (BMI), body fat percentage, visceral adipose tissue (VAT), and abdominal SAT (Table 1). They also had higher serum alanine aminotransferase (ALT), gamma-glutamyl transferase (GGT), non-high-density lipoprotein (HDL) cholesterol, fasting blood glucose, and blood pressure (Table 1). The participant with diabetes and 4 of the 6 participants with prediabetes were NAFLD positive.

Individuals with hepatic steatosis had lower cold-induced percentage decline in SCV PDFF

Study participants who were NAFLD positive had lower BAT activity compared with NAFLD-negative participants (2.5% [1.2%–5.7%] versus 4.7% [2.7%–8.1%], p = 0.019) (Figure 1A; Table 1). Consistent with that, a lower cold-induced percentage decline in SCV PDFF was associated with higher hepatic fat (rho = −0.41, p = 0.001) (Figure 1B). The relationship of the decline in SCV PDFF with hepatic fat remained significant, independent of age, gender, and either total body fat or VAT area (Table 2). Our data suggest that in humans, as in rodents,,, higher BAT activity is associated with lower hepatic fat accumulation, and this relationship is not fully explained by age, gender, or adiposity.

Figure 1

Individuals with hepatic steatosis have lower cold-induced percentage decline in supraclavicular proton density fat fraction

(A) Comparison of the cold-induced percentage decline in supraclavicular (SCV) proton density fat fraction (PDFF) (%) between those without (non-alcoholic fatty liver disease [NAFLD] negative, n = 30) and with (NAFLD positive, n = 29) hepatic steatosis; data are presented as mean ± standard error of the mean (SEM) for each group.

(B) Relationship of the cold-induced percentage decline in SCV PDFF (%) with hepatic PDFF (%) in men (blue circles) and women (red triangles).

†p < 0.05, Mann-Whitney U test; §p < 0.05, Spearman’s correlation.

Table 2

Hierarchical multivariate analysis for the predictors of hepatic fat (age, gender, cold-induced percentage decline in SCV PDFF, total body fat, and visceral adipose tissue)

	Coefficients							Model
	N	B	SE	β	95% CI		p	R2	adjusted R2	p
					lower bound	upper bound
Model 1A (+ cold-induced percentage decline in SCV PDFF)	59							0.338	0.289	< 0.001∗
Age (years)		0.085	0.049	0.221	−0.013	0.183	0.087
Gender (male)		0.430	0.899	0.056	−1.373	2.233	0.635
Cold-induced percentage decline in SCV PDFF (%)		−1.153	0.318	−1.172	−1.790	−0.516	0.001∗
Cold-induced percentage decline in SCV PDFF (%)2		0.067	0.023	0.896	0.020	0.113	0.006∗
Model 1B (+ body fat percentage)	59							0.433	0.380	< 0.001∗
Age (years)		0.051	0.047	0.132	−0.043	0.145	0.284
Gender (male)		1.256	0.885	0.164	−0.519	3.031	0.162
Cold-induced percentage decline in SCV PDFF (%)		−0.885	0.310	−0.900	−1.507	−0.264	0.006∗
Cold-induced percentage decline in SCV PDFF (%)2		0.061	0.022	0.822	0.018	0.105	0.007∗
Body fat (%)		0.129	0.043	0.435	0.042	0.216	0.004∗
Model 2B (+ VAT)	56							0.504	0.454	< 0.001∗
Age (years)		−0.013	0.061	−0.030	−0.135	0.109	0.836
Gender (male)		−0.336	0.845	−0.043	−2.033	1.360	0.692
Cold-induced percentage decline in SCV PDFF (%)		−0.760	0.309	−0.769	−1.380	−0.139	0.017∗
Cold-induced percentage decline in SCV PDFF (%)2		0.051	0.022	0.690	0.008	0.095	0.021∗
VAT (cm2)		0.081	0.022	0.556	0.036	0.126	0.001∗

Quadratic function of cold-induced percentage decline in SCV PDFF (%) was used, as it fits better than the linear model (linear R2 = 0.173 versus quadratic R2 = 0.301). SCV, supraclavicular; PDFF, proton density fat fraction; VAT, visceral adipose tissue.

p < 0.05.

Cold-induced percentage decline in SCV PDFF, age, and adiposity

We also confirmed that the cold-induced percentage decline in SCV PDFF was related to demographic characteristics and measures of metabolic health previously reported to be related to 18F-FDG uptake measured using positron emission tomography (PET)/computed tomography (CT).,,, Lower age (rho = −0.48, p < 0.001) (Figure 2A), total body fat (r = −0.62, p < 0.001) (Figure 2B), waist circumference (r = −0.47, p < 0.001) (Figure 2C), BMI (r = −0.46, p < 0.001), and VAT (rho = −0.60, p < 0.001) were associated with greater cold-induced percentage decline in SCV PDFF. In a multivariate regression analysis, the relationship between cold-induced percentage decline in SCV PDFF and total body fat remained significant independent of age and gender (R2 = 0.514, p < 0.001). Men had a greater cold-induced percentage decline in SCV PDFF than women (mean difference 2.3%, p = 0.022), but this relationship was no longer significant (p = 0.479) when age and total body fat were included in the multivariate regression model.

Figure 2

Cold-induced percentage decline in supraclavicular proton density fat fraction, age, adiposity, and glucose homeostasis

(A–F) Relationship of the cold-induced percentage decline in supraclavicular (SCV) proton density fat fraction (PDFF) (%) with age (A), body fat percentage (B), waist circumference (C), hemoglobin A1c (HbA1c) (D), fasting glucose (E), and 2 h 75 g oral glucose tolerance test (OGTT) glucose (F) in men (blue circles) and women (red triangles).

∗p < 0.05, Pearson’s correlation; §p < 0.05, Spearman’s correlation. See Table S3.

Cold-induced percentage decline in SCV fat fraction and glucose homeostasis

A greater cold-induced percentage decline in SCV PDFF was associated with lower HbA1c (rho = −0.31, p = 0.017), fasting plasma glucose (rho = −0.29, p = 0.026) and 2 h 75 g oral glucose tolerance test (OGTT) plasma glucose (rho = −0.42, p = 0.018) (Figures 2D–2F). The latter remained related to cold-induced percentage decline in SCV PDFF while controlling for age, gender, and total body fat (rho = −0.44, p = 0.022). These data suggest that in addition to the relationship with hepatic fat content, higher cold-stimulated BAT activity is linked to better glucose homeostasis.

Relationship of gut microbiota to cold-induced percentage decline in SCV fat fraction

To determine if gut microbiota characteristics might link lower BAT activity to higher hepatic fat content in humans, as it does in rodents, we examined these relationships in our cohort. Alpha diversity of the fecal microbiota, as measured using the Shannon and inverse Simpson indices, was not related to cold-induced percentage decline in SCV PDFF and was not different between high- and low-BAT groups, and this was true whether or not age and total body fat percentage were included in the model (Figure 3A; Table S3). Furthermore, there was no distinct clustering of BAT groups on a principal coordinates analysis (PCoA) plot using the Bray-Curtis dissimilarity metric (Figure 3B), and there was no significant relationship of BAT groups with microbial community structure (permutational multivariate analysis of variance [PERMANOVA] R2 = 0.014, p = 0.74). The presence of three genera was related to cold-induced percentage decline in SCV PDFF after controlling for age, although these relationships were tempered when total body fat was included in the regression model (Table S4), and these bacteria were present in only one participant in our cohort. The only genus with a different abundance between high- and low-BAT groups was Succiniclasticum, which was more abundant in the low-BAT group (Table S4). It was, however, only present in two participants in the cohort and therefore does not appear in Figure 3C. Collectively, although the abundance of few bacteria was related to the cold-induced percentage decline in SCV PDFF and different between BAT groups (Table S4), these bacteria were present in only a few participants in our cohort.

Figure 3

Relationship of gut microbiota to cold-induced percent decline in supraclavicular fat fraction

(A) Alpha diversity between brown adipose tissue (BAT) groups (low and high BAT activity).

(B) Principal coordinate analysis (PCoA) plot on Bray-Curtis dissimilarity distances between BAT groups.

(C) Relative abundance of top 25 genera between BAT groups.

For (A), data are represented as the median, interquartile range, and 95% range of the data. The colored points and lines are the point estimates of the means from the regression models and the 95% confidence intervals of the model estimates of the means, respectively. For (B), axes 1 and 2 captured 16.4% and 7.3% in the variation between samples, respectively. See Tables S3 and S4 and Figure S3.

To examine whether the lack of differences between groups might be a function of our sample preparation and analysis, we compared microbiome signatures of people with and without NAFLD. Consistent with previous studies,, gut microbiota characteristics differed on the basis of NAFLD status in several distinct areas. Alpha diversity was lower in the fecal samples from individuals with NAFLD compared with those with normal hepatic fat content, independent of age and total body fat (p = 0.022 and p = 0.046 for the Shannon and inverse Simpson indices, respectively) (Figure S3A; Table S3). Similarly, and consistent with previous findings,, there was a distinct clustering between NAFLD groups on the PCoA plot (Figure S3B). Additionally, NAFLD status had a small (3%) but significant (p = 0.012) effect on the variation in gut microbial community structure as measured using the Bray-Curtis dissimilarity index. In summary, although we identified characteristics of the gut microbiota that differed between those with and without NAFLD, we identified no characteristics that were related to both BAT activity and hepatic steatosis.

Microbial transfer to mice

To assess whether BAT activity in humans is microbially transmissible, germ-free mice (three or four per donor) were colonized with fecal material from human donors who had high BAT (n = 3) or low BAT (n = 4) activity (Table S5). Our laboratory has previously shown that this method of microbial transfer and the 8 week host exposure time is sufficient to transmit a microbial community that can promote insulin resistance and poor glucose control when colonizing germ-free mice. After 8 weeks, fecal microbiota of stool samples from the colonized mice was found to cluster with the stool sample from their respective donor (PERMANOVA R2 = 0.39, p = 0.0001) (Figure S4A). Using an UpSet plot to evaluate how many amplicon sequence variants (ASVs) were shared between donor samples and recipient mice, we found 26 and 15 ASVs were uniquely shared between matched donors and recipients of high-BAT and low-BAT samples, respectively (Figure S4B). This was compared with 3 ASVs in discordant pairs of donors and recipients (i.e., high-BAT donors and low-BAT recipient mice), suggesting successful transfer of unique microbiota in each group. The colonized mice who received fecal material from high-BAT donors had similar body mass (Figure 4A), fat mass (Figure 4B), and oxygen consumption (Figure 4C) as the colonized mice who received fecal material from low-BAT donors. Also, recipient mice from high- and low-BAT donors had comparable oxygen consumption and interscapular surface temperature at baseline and when challenged with the β3-adrenergic receptor agonist CL-316,243, which specifically increases UCP1-mediated thermogenesis (Figures 4D–4F). Similar results were collected using C57Bl6J-colonized specific-pathogen-free (SPF) mice gavaged with fecal material from human high- or low-BAT donors (Figure S5). These data indicate that transmission of fecal microbes from humans discordant for BAT activity did not result in alterations in body mass, fat mass, or BAT activity in the recipient mice.

Figure 4

Human brown adipose tissue activity and NAFLD status are not transmissible via gut microbes in germ-free mice

(A and B) Weekly body mass (A) and percentage fat mass (B).

(C) Average oxygen consumption (Vo2) during light and dark cycles and over a 24 h period (average of light and dark) after 5 weeks of colonization.

(D and E) Oxygen consumption (D) and dorsal interscapular surface temperature of anesthetized mice following saline or CL-316,243 administration after 6–7 weeks of colonization (E).

(F) Representative infrared images from each group after saline or CL-316,243.

(G) Liver triglycerides after 7–8 weeks of colonization (n = 9 H-BAT/NAFLD− and n = 12 L-BAT/NAFLD+).

Data are expressed as mean ± SEM. Significance was determined using Student’s t test or two-way repeated-measure analysis of variance (ANOVA) and Sidak’s post hoc test; †p < 0.05 versus saline. See Figures S4 and S5 and Table S5.

Our findings contrast with those of studies done in rodents in which browning of white adipocytes and enhancement of BAT activity can be transmissible via fecal microbes. Cold-exposed mice have augmented BAT activity and distinctive gut microbial profiles.48, 49, 50 Microbiota transfer from cold-exposed mice to germ-free mice promotes white adipocyte browning through the upregulation of UCP1 expression. After cold stimulation, exposed mice have increased production of bacterially derived compounds such as butyrate that can augment BAT function. Although these results suggest that BAT metabolism can be altered via the gut microbiota in rodents, we have not identified similar findings in humans. No detectable changes in body composition, whole-body energy expenditure, or UCP1-mediated thermogenesis were identified in mice that received gut microbes from humans with differing levels of BAT activity. These findings are consistent with the lack of association of the characteristics of the fecal microbiota with BAT activity.

Limitations of the study

Although this study is one of the largest human studies to use MRI in the evaluation of BAT in cold-stimulated conditions, metabolic health, and gut microbiota, our study is not without limitations. It is uncertain how our findings might apply to groups that were not included in our study, including older adults and individuals with diabetes or other conditions managed with pharmacotherapy. We also did not have individuals with severe obesity or high levels of liver fat in our study, because of the weight limitations of the MRI scanner (weight must be less than 136 kg) and the exclusion of individuals using pharmacotherapy for metabolic conditions. Despite these limitations, we had findings related to the relationship between the gut microbiome and NAFLD consistent with other studies that included individuals with more severe NAFLD (Figure S3C).,,, Another limitation of our study was the relatively small sample size. The study was powered to evaluate the relationship of BAT activity and hepatic fat accumulation, independent of age, gender, and total body adiposity as outlined in STAR Methods. Further studies in a larger group of individuals with higher hepatic fat levels would contribute to our understanding of the generalizability of our findings to individuals who are older and to those with more severe NAFLD. Importantly, we examined the microbiota-NAFLD relationship independent of total body adiposity. Most previous studies have controlled for BMI instead of adiposity., The limitations of BMI as a surrogate for adiposity have been previously published. We did not collect ethnicity data. Ethnicity has been shown to affect the gut microbiome but not BAT activity., Although it is unlikely that ethnicity influenced our conclusions related to BAT activity, further investigation in an ethnically diverse group would be needed to demonstrate this definitively. Last, the human microbiota can be influenced by diet. Although the self-reported macronutrient and dietary fiber intake did not differ between those with and without NAFLD (Table S2), we did not standardize dietary intake, and thus differences may have contributed to the lack of relationship between the fecal microbiota and BAT activity.

Conclusions

The present study demonstrates that in human adults, lower BAT activity (cold-induced percentage decline in SCV PDFF) is associated with greater hepatic fat content and higher measures of glycemia, independent of age, gender, and total body fat. Contrary to our hypothesis, we did not identify characteristics of the gut microbiota related to BAT activity or evidence that BAT activity can be transmitted via fecal transplantation from humans to germ-free mice. Thus, our findings highlight a clear link between BAT activity and hepatic fat content, supporting the premise that stimulation of BAT activity may be a potential therapeutic target for the management of NAFLD in humans. We did not, however, find evidence that the fecal microbiota plays an intermediary role in the relationship between BAT activity and NAFLD.

STAR★Methods

Key resources table

Resource availability

Lead contact

Requests for further information regarding resources utilized or data should be directed to and will be fulfilled by the Lead Contact, Dr. Katherine M. Morrison (kmorrison@mcmaster.ca).

Materials availability

This study did not generate any new reagents or materials.

Experimental model and subject details

Study population

Participants in this two-visit, observational, cross-sectional study were enrolled from the City of Hamilton, Ontario, Canada, and the surrounding area between June 2016 and March 2018. Male and female participants aged 18 to 60 years were included in the study unless they met any of the following exclusion criteria: use of nicotine or tobacco products, self-reported alcohol intake of ≥ 7 drinks/week and ≥ 3 drinks/day for males or ≥ 2 drinks/day for females, contraindications for MRI imaging (claustrophobic, implanted metal, metallic injuries, recent tattoos obtained outside of Canada), inability to fit inside or greater than the weight limit (i.e., more than 136 kg) of the MRI scanner, pregnancy or a breastfeeding mother, use of one of the excluded medications including adrenergic and serotoninergic agents, metformin and other diabetes medications (Table S6), a history of bariatric surgery, liver transplantation or medical conditions i.e., liver diseases and untreated hormonal imbalance (Table S6).

Method details

Study approval

This study was approved by the Hamilton Integrated Research Ethics Board (HiREB) and was conducted utilizing Good Clinical Practise (GCP) principles. All participants provided informed consent. All the animal experiments were approved by the McMaster University Animal Ethics Committee and conducted under the Canadian guidelines for animal research.

Study visits

Participants were asked to fast for at least 8 hours, to refrain from vigorous physical activity for at least 48 hours, avoid intake of caffeine for at least 12 hours, and to abstain from serotonin-rich foods (e.g., banana, tomato, kiwi, walnut, avocado, pineapple, and plum) for at least 24 hours prior to both visits. Study participants completed a prospective, dietary diary for the 24-hour period prior to Visit 2. Visits 1 and 2 took place in the morning and were scheduled at least 7 days apart, but ideally within one month of each other. To standardize Visit 2 to the same point in the menstrual cycle for female participants who were pre-menopausal, the visit occurred within 7 days after the commencement of menstruation.

Visit 1 occurred at the McMaster University Medical Centre, Hamilton, ON, Canada. Measures included anthropometry (height, weight, and waist circumference), blood pressure, body composition assessment, and a 75 g OGTT, with blood samples collected at baseline and 2 hours. The one study participant with a known diagnosis of type 2 diabetes (n = 1) had a fasting blood sample only and no OGTT was done.

Visit 2 took place at the Imaging Research Centre at St. Josephs Healthcare, Hamilton. The mean outdoor temperature at 1 hour, 24 hours, and 7 days prior to this visit were obtained retrospectively from data collected at the McMaster Weather Station, operated by the Hydrometeorology and Climatology Lab in the School of Geography and Earth Sciences, McMaster University. Upon arrival, participants changed into standardized cotton tank tops and shorts. They were then acclimated to room temperature for 30 minutes. The temperature in the room where participants underwent acclimation and cold exposure was recorded (Wireless Forecast Station with Pressure History Model WS-9037U-IT; La Crosse Technology, La Crosse, WI, USA). The baseline (pre-cold) MRI scan was conducted of the SCV region, liver, and abdomen. Participants were then exposed to a standardized cold exposure of 18°C for 3 hours using a full-body water-perfused suit. Shivering intensity was measured using electromyography (EMG) as described below. The MRI scan of the neck was repeated immediately after the cold exposure.

Cold exposure protocol

The cold exposure protocol was based on a protocol previously shown to be adequate for BAT stimulation, measured using 18F-FDG PET/CT.,, Participants were fitted with a high-density liquid conditioned suit (LCS; Two Piece, Med-Eng, Ottawa, ON, Canada), and 18°C water was circulated for 3 hours using a temperature- and flow-controlled circulation bath (Isotemp 6200 R28, Fisher Scientific, Waltham, MA, USA). Inlet and outlet water temperature entering and exiting the suit were recorded during the cold exposure at 15 s intervals with a dynamic range of 16 bit using a data logger (PowerLab; ADInstruments, Sydney, Australia) connected to two quick-disconnect thermocouples (TMQSS-020G-2; OMEGA Engineering, Stamford, CT, USA) fixed to the inlet and outlet manifolds. The inlet to outlet temperature difference (Δ Outlet-inlet) of the suit was calculated for the last 30 minutes of the cold exposure.

Measurement of shivering intensity with electromyography (EMG)

The shivering intensity was measured at baseline and during cold exposure using surface EMG (Trigno Wireless System, Delsys, MA, USA with Trigno Snap-Lead Sensor connected to pre-jelled Norotrode 20 Bipolar SEMG electrodes, Myotronic, WA, USA). Electrodes were placed on the bellies of the trapezius, vastus lateralis, and vastus medialis muscles (right side only). A series of maximal voluntary contractions (MVC) were performed on each muscle (three trials of maximum contraction to be reached by verbal encouragement alternate with 30 s rest between trials) and the shivering intensity of each muscle was calculated from root-mean-square (RMS) values calculated from raw sEMG, normalized to MVC RMS and corrected for baseline RMS measured at room temperature while lying down, as previously described.

MRI acquisition protocol

MRI-derived proton density fat fraction (MRI-PDFF) measures the ratio of the total density of fat mobile protons to the total density of fat and water mobile protons and therefore mirrors the triglyceride concentration in tissues. Accordingly, MRI-PDFF differentiates BAT and surrounding tissues including white adipose tissues and muscles based on their physical properties. Thus, BAT-containing depots showed lower PDFF than white adipose tissue.,

MRI-PDFF can be used to measure BAT activity; the reduction in the supraclavicular (SCV) PDFF in response to cold exposure is correlated with 18F-FDG uptake measured using PET-CT., Further, cold-stimulated reductions in SCV-PDFF have been corroborated with histological findings, UCP1 immunohistochemistry,, and observations that triglyceride mobilization is central to BAT thermogenesis. MRI-PDFF is also a well-established, non-invasive, quantitative, and accurate means to measure hepatic fat content.

MRI scans were performed using a 3-Tesla (T) whole-body MRI scanner (Discovery 750; GE Healthcare, Waukesha, WI, USA). Quantitative terative Decomposition of water and fat with cho symmetry and east-squares estimation (IDEAL-IQ) was used since it is a confounder-corrected 3D gradient multi-echo MRI sequence that provides an accurate measure of tissue triglyceride content by using multiple spectral modeling of adipose tissues while accounting for T2∗ decay. This pulse sequence generates six distinct image contrasts: water-only, fat-only, in-phase, out-of-phase, corrected proton density fat fraction (PDFF), and R2∗ images. LAVA-FLEX is a 3D gradient dual-echo MRI sequence that generates four distinct image contrasts: water-only, fat-only, in-phase, and out-of-phase. Due to its nature, the LAVA-FLEX pulse sequence provides a higher resolution image with sharper tissue boundaries, compared to IDEAL-IQ.

BAT activity

Supraclavicular measurements of PDFF were acquired using the IDEAL-IQ pulse sequence. Using a Head/Neck/Spine (HNS) coil, with an additional attachment that provided signal from the anterior portion of the chest, axial images were taken from the C2/C3 disc to the T4/T5 disc (IDEAL-IQ, slice thickness 4mm, 50 slices, flip angle = 4°, TE/TR = 1.2/7.8ms, FOV = 380 mm, matrix size = 224x192 (i.e., image resolution = 1.70 × 1.98 × 4mm), acceleration factor 2, scan time = 2.7 min). Segmentation and analysis of all MRI images were conducted using Analyze Pro (Version 1; Mayo Clinic, Biomedical Imaging Resource, AnalyzeDirect, Overland Park, KS, USA) by a single reader. The SCV was segmented using semi-automated and manual segmentation tools and the accuracy of segmentation was confirmed by an experienced radiologist (NS). The vertebral column was used as the reference for the start (C5-C6 disc) and end (T1-T2 disc) of segmentation. A fat mask, generated from the fat-only images, was first applied to the PDFF map at the C7-T1 disc to isolate adipose tissues and exclude background noise from the MRI image. A FF threshold set at 30 to 100% was then applied to the PDFF image. A region of interest (ROI) was manually drawn over the SCV region defined by the trapezius muscle posteriorly, the sternocleidomastoid muscle medially, and the clavicle inferiorly. The ROIs were then post-processed using a one-time 2D erosion (1x3 voxels) to correct for any inherent partial volume effects in the images. T2∗ images were generated from the R2∗ data by use of the relationship: T2∗ = 1/R2∗. A T2∗ mask was applied to the FF map to help differentiate BAT from white adipose tissue (WAT) by isolating voxels with a T2∗ between 2 and 25ms, based on a published report that tissue with a T2∗ value of ≥ 26ms consists mostly of muscle, fluid or white adipocytes. The lower range of 2ms was selected due to the MRI sequence’s limitation in detecting very low T2∗ values. All voxels that satisfied the above criteria were averaged and classified as SCV PDFF.

Hepatic fat

Hepatic fat was measured from axial scans obtained in a single breath-hold using a 32-channel torso array coil (NeoCoil, Pewaukee, WI, USA) from 7cm above to 7cm below the L4-L5 disc (IDEAL-IQ, slice thickness 8mm, 32 slices, flip angle = 3°, TE/TR = 1.0/6.5ms, FOV = 340 mm, matrix size = 160x160 (i.e., image resolution = 2.1 × 1.7 × 8mm), acceleration factor 2 × 1.5, scan time = 18 s). A multi-slice segmentation that included the entire liver was undertaken using Analyze Pro software. ROIs were drawn over the entire liver using a tool that “snaps” to the edges of regions where changes in voxel intensities are high. The in-phase image was used as a reference when the boundaries of the liver were not clearly defined in the fat-fraction image. The ROIs were then post-processed using a 2D erosion (3x3 voxels) to correct for partial volume effects. All voxels that satisfied the above criteria were averaged and classified as hepatic fat.

Abdominal subcutaneous and visceral adipose tissue areas

For abdominal subcutaneous adipose tissue (SAT) and visceral adipose tissue (VAT) were measured from a single slice located 6 cm above the L4-L5 disc. Abdominal SAT PDFF was calculated from a single slice at L4-L5 disc obtained in a single breath-hold scan (IDEAL-IQ, Slice thickness 5mm, 32 slices, flip angle = 3°, TE/TR 6.1/0.9ms, FOV = 400mm, matrix size = 160x160 (image resolution 2.5 × 2.0 × 5mm), acceleration factor 2 × 1.5, scan time = 18 s). A semi-automated tool (i.e., 2D seeded region growing) was used to segment abdominal SAT. Voxels with low PDFF values (below 30%) were excluded to isolate adipose tissues from muscle. A 2D erosion (3x3 voxels) was applied to attenuate partial volume effects. All voxels that satisfied the above criteria were averaged and classified as SAT PDFF. To measure the area of abdominal SAT and VAT, the same slice was used from the higher resolution, single breath-hold LAVA-FLEX scan (LAVA-FLEX, slice thickness 5mm, 32 slices, flip angle = 3°, TE/TR = 1.3/4.1ms, FOV = 400mm, matrix = 320x224 (i.e., spatial resolution = 1.25 × 1.78 × 5mm), acceleration factor 1.5, scan time = 14 s). A fat mask was applied to exclude non-adipose tissues from analysis while SAT and VAT were separated semi-automatically based on seed points. Retroperitoneal adipose tissue was excluded manually from VAT. Finally, the ROIs were post-processed using 2D erosion (3x3 voxels) to attenuate partial volume effects.

Anthropometric measures

Weight, height, and waist circumference were each measured three times by trained research personnel and these values were averaged. Weight (kg) was measured using an electronic platform scale (BMI Scale Model 882; Seca, Hamburg, Germany) to the nearest 0.1 kg. Height (cm) was measured using a wall-mounted stadiometer (Height Measuring Rod Model 240; Seca, Hamburg, Germany) to the nearest 0.1 cm. BMI was calculated as weight in kilograms divided by height squared in meters. Waist circumference (cm) was measured at the midpoint between the costal margin and the iliac crest at the end of expiration using a weighted measuring tape (Pull Type Spring Scale; Ohaus, Parsippany, NJ, USA) set at 750 g.

Body composition

Body composition (total body fat (%), fat mass (kg), lean mass (kg)) was assessed using dual-energy X-ray absorptiometry (DEXA) (Lunar Prodigy Advance 8743; GE Healthcare, Waukesha, WI, USA). All scans were reviewed by one individual (KMM) to ensure the consistency of ROIs.

Blood samples

Blood samples were collected at baseline (fasting) and 2-hours after the consumption of a 75 g glucose drink and were analyzed in the Clinical Laboratory at Hamilton Health Sciences (Hamilton, ON, Canada). Fasting total cholesterol, HDL-cholesterol, triglyceride were measured using Architect C16 (Abbott) on serum samples and LDL-cholesterol was calculated using the Friedewald equation. Alanine aminotransferase (ALT), aspartate aminotransferase (AST), fasting, and 2h plasma glucose were measured using Architect C4000 or c4100 (Abbott, Chicago, IL, USA) on heparinized plasma samples. HbA1c was analyzed in whole-blood samples on the Capillary 2 Flex Piercing (Sebia, Lisses, France) instrument. Dysglycemia was based on HbA1c and was classified according to the American Diabetes Association guidelines, prediabetes (HbA1c = 5.7 – 6.4%) and diabetes (HbA1c ≥ 6.5%).∗

Dietary intake

All study participants were instructed to complete a diary of all foods and drinks consumed (measured and including brand names) in the 24 hour period prior to the study visit. Information from the diet diary was entered into the nutritional software Food Processor by ESHA. Total calories, macronutrients (protein, carbohydrate, total and saturated fat), and dietary fiber intakes were calculated utilizing the software.

Stool samples

Stool samples were collected by participants within 48 hours of Visit 2. Participants were instructed to immediately transfer 10-20 g of fecal sample into the 50ml sterile specimen container and storage packaging provided, store in a cooler pack in their freezer, and bring the frozen samples to Visit 2. Samples were delivered frozen to the laboratory and stored at −20◦C until thawed and aliquoted. One portion (approximately 200mg) was reserved for nucleic extraction and up to five 2 mL cryovials were stored for future analyses or experiments. All aliquots were kept frozen at - 80◦C.

Microbiota profiling

DNA extraction was carried out as described previously, with few modifications, from 0.2g of stool sample that was mechanically lysed with 2.8mm ceramic beads (Mo Bio Laboratories Carlsbad, CA) and 0.1mm glass beads (Mo Bio Laboratories, Carlsbad, CA) for 2 cycles of 3 minutes at 3000rpm with a 45 s delay between cycles in 800 μL of 200mM of monobasic NaPO4 (pH = 8) and 100 μL of guanidinium thiocyanate buffer. Samples were further centrifuged at 13,500 g for 5 minutes and processed using the MagMAX-96 DNA Multi-sample kit (Life Technologies, Carlsbad, CA) via the MagMAX Express 96-Deep Well Magnetic Particle Processor (Applied Biosystems, Foster City, CA). Purified DNA was used for amplification of the 16S rRNA gene using paired-end reads of the V3 region. Amplification of the bacterial 16S rRNA gene for the germ-free mice microbiota was done for the V3-V4 region. 50 ng of DNA was used as a template with 1U of Taq, 1x buffer, 1.5 mM MgCl2, 0.4 mg/mL BSA, 0.2 mM dNTPs, and 5 pmoles each of 341F (CCTACGGGAGGCAGCAG) and 518R (ATTACCGCGGCTGCTGG) Illumina adapted primers, as previously described. The reaction was carried out at 94◦C for 5 minutes, 25 cycles of 94◦C for 30 s, 50◦C for 30 s, and 72◦C for 30 s, with a final extension of 72◦C for 10 minutes. PCR products were then visualized on a 1.5% agarose gel. Positive amplicons were normalized using the SequalPrep Normalization Plate kit (Thermo Fisher, Waltham, Massachusetts, USA). Triplicate amplifications were pooled for 250bp paired-end MiSeq Illumina sequencing in the McMaster Genomics Facility. Sequence data were processed through the DADA2 pipeline (divisive amplicon denoising algorithm 2) to produce ASVs (amplicon sequence variants). After processing, we had 4,881,511 reads representing 2,362 ASVs from the human stool samples. The number of reads per sample ranged from 11,497 to 143,413 with a median read number of 87,133,5. We removed all ASVs that were suspected to be host DNA including any ASVs (i) assigned to Kingdom Eukaryota, (ii) that lacked Phylum assignment, or (iii) that were assigned to Family Mitochondria. We also removed all ASVs with mean abundance < 5. After filtering, we had 4,771,682 reads representing 574 ASVs. The number of reads per sample ranged from 10,160 to 141,152 with a median of 84,666.

Germ free mice samples were sequenced on a 2x300bp paired-end MiSeq Illumina run. To compare germ-free mice data with the human data, read 2 (R2) sequences from the mice microbiota sequencing data and donor stool sequencing data were trimmed to the V3 region using cutadapt. This data underwent the same filtering steps as described above. After filtering, we had 2,927,423 reads representing 1291 ASVs. The number of reads per sample ranged from 27,010 to 134,635.

Filtered ASVs were processed for analysis using phyloseq package in the R statistics software. Alpha diversity (within-sample diversity) indices of community evenness and richness (Shannon Wiener index and Inverse Simpson index) were calculated after rarefying the samples to minimum sample depth. Bray-Curtis dissimilarity (between-sample diversity) on the relative abundance of all ASVs was calculated and principal coordinate analysis (PCoA) plots were generated with ggplot and were used to visualize the clustering of microbiota communities between BAT groups and NAFLD groups. Non-metric Multidimensional Scaling (NMDS) plots were generated and used to visualize the clustering of communities between human donor stool and stool from fecal transplant recipient mice. Upset plots were generated by considering any given ASV that appeared in any donor or recipient sample, grouped by high or low BAT, as presented and visualized using UpsetR. Relative abundance of genera was calculated by taking the ASV counts assigned for each member and dividing these counts by the sum of all ASV counts within a given sample. Stacked bar charts of the relative abundance of the top 25 genera were plotted using ggplot to visualize the differences within BAT and NAFLD groups.

Murine fecal microbiota transfer experiment

Male germ-free C57BL/6N mice (14-20 weeks old) were exported from McMaster University’s Farncombe Gnotobiotic Unit and immediately and continually colonized twice weekly using fecal slurries from human donors over 8 weeks. Fecal material was diluted 1:10 (wt:vol.) in sterile saline and dissolved mechanically. A brief centrifugation separated undissolved fecal material from the supernatant, which was aliquoted and frozen at −80°C until needed. Mice were gavaged 200μL of the fecal slurry twice weekly from their assigned donor. Male donors were selected based on BAT and NAFLD status, and each donor colonized 3-4 germ free mice. Mice were housed individually at thermoneutrality (30°C), given access to irradiated standardized chow and autoclaved water ad libitum, and handled in a level II biosafety hood to prevent bacterial contamination. Weekly body weight and adiposity measurements using the Bruker Minispec LF90-II BCA-Analyzer were taken. At week 5, metabolic monitoring was conducted at thermoneutrality using a Comprehensive Lab Animal Monitoring System (CLAMS, Columbus Instruments, OH, USA). Between weeks 6-7, UCP-1 mediated thermogenesis using CL-315,243 was assessed as previously outlined. Briefly, anesthetized mice were injected i.p. with either saline or CL-316,243 and placed in an enclosed stationary treadmill dorsal side up. Oxygen uptake was monitored using CLAMs with air sampled every 5 s. After twenty minutes, mice were removed from the treadmill and a dorsal thermal image was acquired using an infrared camera. These experiments were also conducted in control, wild-type C57BL/6J mice gavaged with fecal slurries three times a week from human high or low BAT donors, or their fecal material (i.e., autologous controls). Non-moving VO2 was calculated from CLAMS data. Mice were sacrificed between weeks 7-8 and harvested tissues were flash-frozen for further analyses. Hepatic TGs were extracted from 50mg of tissue and quantified using an adapted Bligh and Dyer lipid extraction method and quantified using a glycerol assay.

Quantification and statistical analysis

To evaluate BAT activity accounting for the pre-cold SCV PDFF, the percent change in SCV PDFF was calculated as [(pre-cold SCV PDFF - post-cold SCV PDFF)/ (pre-cold SCV PDFF)] ∗100 and is introduced in the main text as cold-induced percent decline in SCV PDFF. For microbiota analysis, participants were categorized into high BAT (n = 27, male = 21) and low BAT (n = 33, male = 16) groups. Individuals were considered with low BAT activity if the cold-induced percent decline in SCV PDFF was within 1 SD of zero.

Statistics

The sample size was based on our primary objective to examine the relationship of BAT activity to hepatic fat accumulation in adults, considering inclusion of the potential covariates age, gender and total body fat content. Based on the work of Harris, the suggested number of participants should equal 50 + m (where m is the number of independent variables). Thus, a sample size of 54 participants was needed to have sufficient power to address the primary question. Additional analyses included examining the relationship of the characteristics of the gut microbiome to BAT activity which we undertook in two ways. Determination of sample size for human microbiota studies is challenging but referring to previous similar studies can be used to estimate sample size for a new study. To the best of our knowledge, there have been no studies that examined microbiota and BAT in humans. The sample size to examine the relationships of the fecal microbiota to NAFLD is similar to a previous study. For the microbial transfer, 3-4 germ free or SPF mice per donor were previously found to be sufficient to show the transfer of microbial phenotype.

SPSS Statistics (version 27; IBM, North Castle, NY, USA), GraphPad Prism (version 8; GraphPad Software, La Jolla, CA, USA), RStudio (version 1.2.1335; R Core Team, 2013), and R (Foundation for Statistical Computing, Vienna, Austria) were used for all statistical analyses and for graphing. Normality was assessed by following the procedures outlined by Tabachnick and Fidell where variables with Zskewness and/ or Zkurtosis ≥ |3.29|were classified as non-normal. Non-parametric tests were used for analyses of non-normally distributed data. Data were presented as n (%) for categorical variables and mean (SD) for normally distributed variables or median [IQ1, IQ3] for skewed variables. Paired Student’s t test was used to compare MRI measures before and after cold exposure. Pearson (r) correlation coefficients were used to assess the linear association between normally distributed variables. Spearman (rho) correlation coefficients were used to assess the association between non-normally distributed variables and for non-linear monotonic correlations. Independent Samples t test and Mann-Whitney U test were used to compare differences of normally distributed and skewed variables respectively between NAFLD and BAT groups. Hierarchical multivariate regression analysis was used to examine whether associations between cold-induced percent decline in SCV PDFF (independent variable) and liver fat (dependent variable) are independent of covariates (age, gender, body fat, or VAT). Multivariate regression analysis was also done to examine the predictors (independent variables; age, gender, and body fat percentage) of cold-induced percent decline in SCV PDFF (dependent variable). Non-parametric partial correlation was used to examine the association between BAT activity and glycemia measures controlling for the effects of the covariates (age, gender, and body fat percentage). Curve fitting was used to assess whether relationships between variables were linear or non-linear. Quadratic functions were used in the regression models for cold-induced percent decline in SCV PDFF in relation to hepatic fat and for body fat in relation to cold-induced percent decline in SCV PDFF since they fit better than the linear models. A two-tailed p value of less than 0.05 was considered significant.

Hierarchical regression was used to examine the difference in microbiota alpha diversity indices (the outcome) between BAT activity and NAFLD groups. Three models were used, model 1; without controlling for covariates, model 2; controlling for age, and model 3; controlling for age and total body fat. Permutational Analysis of Variance (PERMANOVA) test with 999 permutations using the adonis function in the vegan package of R software was performed to test the effect of BAT and NAFLD groups on beta diversity controlling for age and total body fat.

DESeq2, a negative binomial generalized linear regression model, was used to examine the differentially abundant genera in relation to cold-induced percent decline in SCV PDFF, between BAT and between NAFLD groups through three models (model 1; without controlling for covariates, model 2; controlling for age and model 3; controlling for age and total body fat). DESeq2 was reported to be one of the best tools for the analysis of count data for differentially expressed genes in terms of precision, accuracy, and sensitivity. DESeq2 provides an estimate log2 fold of change (LFC) of the read counts for each bacterial gene. As pre-cold SCV PDFF and cold-induced percent decline in SCV PDFF have different scales, Z-scores of these two variables were used to make comparable effect sizes across them. As such, for each one-unit change (1 SD) of the cold-induced percent decline in SCV PDFF, we reported the log2 (i.e., the change in each genus or the number of reads that are uniquely assigned for each gene). For BAT and NAFLD groups, we reported the log2 (the change in the ratio between the mean reads for high BAT activity or NAFLD positive to the mean of reads for low BAT activity or NAFLD negative). We used the library depth size factor, i.e., the number of bacterial reads in each sample in the presented analysis. DESeq2 uses a Wald test that is the shrunken estimate of LFC divided by its standard error to give a Z-statistic. To adjust for multiple testing, the false discovery rate was controlled using Benjamini and Hochberg’s method. An adjusted p value of < 0.05 was considered significant.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Human blood samples (serum and plasma)	This study	N/A
Human stool samples	This study	N/A
Critical commercial assays
ALT	Abbott	7D56-21
AST	Abbott	7D81-21
GGT	Abbott	7D65
HDL	Abbott	3K33-21
Cholesterol	Abbott	7D62
TG	Abbott	7D74-21
Glucose	Abbott	3L82
HbA1c	Sebia	https://www.sebia.com/tests/hba1c/
CL-316,243	Tocris	1499 (Batch: 7A/173953)
2,2,2-Tribromoethanol (Avertin)	Sigma-Aldrich	T48402
Tekland Irradiated Global 18% Protein Rodent Diet	Envigo	2918
Glycerol Standard Solution	Sigma-Aldrich	G7793
Free Glycerol Reagent	Sigma-Aldrich	F6428
Ceramic beads	Mo Bio Laboratories	13114- 50
Glass beads	Mo Bio Laboratories	13116-50
MagMax-96 DNA Multi-sample kit	Life Technologies	4413022
SequalPrep Normalization Plate kit	Thermo Fisher	A1051001
Software and algorithms
SPSS v27	IBM	RRID:SCR_002865
Graphpad Prism v8	GraphPad Software	RRID:SCR_002798
RStudio v1.2.1335	Foundation for Statistical Computing	RRID:SCR_000432
Analyze Pro v1	Analyze Direct	https://analyzedirect.com/analyzepro-help-videos/
EMGworks v4.3		https://delsys.com/emgworks
Amide	Amide	http://amide.sourceforge.net/index.html
Other
Discovery MR750 3.0T (MRI scanner)	GE Healthcare	https://www.gehealthcare.com/courses/discovery-mr750-30t
Lunar Prodigy Advance 8743 (DEXA scanner)	GE Healthcare	https://www.gehealthcare.com/products/bone-and-metabolic-health/prodigy
High-Density Liquid Conditioned Two-Piece Suit (cooling suit)	Med-Eng	https://www.med-eng.com/
Isotemp 6200 R28 (refrigerated / heated bath circulator)	Fisher Scientific	https://www.fishersci.com/shop/products/fisher-scientific-isotemp-refrigerated-heated-bath-circulators-6-8-8-6l-115v-6-refrigerated-bath-circulator-6200-r28-28-200-c/13874676
Trigno Wireless System (EMG system)	Delsys	https://delsys.com/trigno
Norotrode 20 Bipolar (EMG electrodes)	Myotronic	https://www.myotronics.com/supplies/all-electrodes/norotrode-20-bipolar-semg-electrode
PowerLab (temperature data logger)	ADInstruments	https://www.adinstruments.com/products/powerlab-daq-hardware
TMQSS-020G-2 (thermocouples)	OMEGA Engineering	https://www.omega.ca/en/temperature-measurement/temperature-probes/probes-with-integral-connectors/p/JMQSS
MATLAB (custom-designed algorithm for EMG analysis)	Mathworks	RRID:SCR_001622
BMI Scale Model 882 (electronic platform scale)	Seca	https://www.scalesgalore.com/seca882.htm
Height Measuring Rod Model 240 (wall-mounted stadiometer)	Seca	https://www.seca.com/en_ee/products/all-products.html#/?category=6
Pull Type Spring Scale (weighted measuring tape)	Ohaus	https://us.ohaus.com/en-US/SpringScales
BPM-200 (blood pressure machine)	BpTRU Medical Devices	https://www.medsource-sw.com/blood-pressure/item-bpm-200/
Head/Neck/Spine (MRI coil)	ScanMed	https://www.scanmed.com/adultheadneckspinemricoil
IDEAL-IQ (MRI sequence)	GE Healthcare	https://www.gehealthcare.com/products/magnetic-resonance-imaging/mr-applications/ideal-iq---body
LAVA-FLEX (MRI sequence)	GE Healthcare	https://www.gehealthcare.com/products/magnetic-resonance-imaging/mr-applications/lava-flex---body
32-channel torso array (MRI coil)	NeoCoil	http://www.neocoil.com/products.php
WS-9037U-IT (weather station)	La Crosse Technology	https://www.lacrossetechnology.com/products/ws-9037u-it
Architect C4000 or c4100	Abbott	https://www.corelaboratory.abbott/int/en/offerings/brands/architect/architect-c4000
Architect C16	Abbott	https://www.corelaboratory.abbott/int/en/offerings/brands/architect
Capillarys 2 Flex Piercing	Sebia	https://www.sebia.com/tests/hba1c/
Comprehensive Lab Animal Monitoring System (CLAMS)	Columbus Instruments	https://www.colinst.com/products/clams-comprehensive-lab-animal-monitoring-system
Bruker Minispec LF90-II BCA-Analyzer	Bruker	https://www.bruker.com/en/products-and-solutions/magnetic-resonance.html
Infrared camera, T650sc, emissivity of 0.98, FLIR Systems	Flir	https://www.flir.ca/support/products/t650sc/#Overview
Experimental models: Organisms/strains
Mouse: Germ Free C57BL/6N	Farncombe Gnotobiotic Unit, McMaster University	Denou et al.59