Mammary ductal epithelium controls cold-induced adipocyte thermogenesis


 Sympathetic activation during cold exposure increases adipocyte thermogenesis via the expression of mitochondrial protein uncoupling protein 1 (UCP1)1. The propensity of adipocytes to express UCP1 is under a critical influence of the adipose microenvironment and varies among various fat depots 2-7. Here we report that cold-induced adipocyte UCP1 expression in a female mouse subcutaneous white adipose tissue (scWAT) is regulated by mammary gland ductal epithelial cells in the adipose niche. Single cell RNA-sequencing (scRNA-seq) show that under cold condition glandular alveolar and hormone-sensing luminal epithelium subtypes express transcripts that encode secretory factors involved in regulating adipocyte UCP1 expression. We term mammary duct secretory factors as “mammokines”. Using whole-tissue immunofluorescence 3D visualization, we reveal previously undescribed sympathetic nerve-ductal points of contact and show that sympathetic nerve-activated mammary ducts limit adipocyte UCP1 expression via cold-induced mammokine production. Both in vivo and ex vivo ablation of mammary ductal epithelium enhances cold-induced scWAT adipocyte thermogenic gene program. The mammary duct network extends throughout most scWATs in female mice, which under cold exposure show markedly less UCP1 expression, fat oxidation, energy expenditure, and subcutaneous fat mass loss compared to male mice. These results show a previously uncharacterized role of sympathetic nerve-activated glandular epithelium in adipocyte thermogenesis. Overall, our findings suggest an evolutionary role of mammary duct luminal cells in defending glandular adiposity during cold exposure, highlight mammary gland epithelium as a highly active metabolic cell type, and implicate a broader role of mammokines in mammary gland physiology and systemic metabolism.


adipogenic 10T1/2 cells with nontransformed mouse mammary gland (NMuMG cells) (derived 174
from "normal" mammary epithelium). Seeding density of as low as 2.5% NMuMG cells resulted 175 in a significant reduction of Ucp1 and beiging potential of 10T1/2 cells compared to the pure 176 culture ( Fig. 3B and Extended Data Fig. 3D-3F). The higher the fraction of NMuMG cells in the 177 co-culture, the lower the relative expression of Ucp1 and other thermogenic genes such as Cox8B 178 and Ppargc1a measured by RT-qPCR (Fig. 3B). To test the ex vivo and in vivo role of ductal 179 epithelial cells in adipose thermogenesis, we compared aged-matched female mgWAT with three 180 complimentary ductal ablated models; i) estrogen receptor-alpha (ERa) knockout mice (Esr1 KO), 181 ii) male mice which lack or possess only rudimentary glandular ducts, iii) inguinal (ducts) or 182 dorsolumbar (no ducts) portion of mgWAT from 5-week-old female mice. First, to test the role of 183 ductal cells in adipose thermogenesis, we isolated SVFs from male iWATs and female mgWATs 184 and differentiated them into beige adipocytes in the presence or absence of isoproterenol. We 185 found that SVFs from male iWATs show markedly higher beiging and isoproterenol-mediated 186 Ucp1 expression compared to EPCAM+ SVFs from female mgWAT ( Fig. 3C and Extended Data 187 Fig. 3G). In agreement with the ex vivo data, female mgWATs showed significantly less 188 expression of cold-induced thermogenic genes such as Ucp1, Cox8b, and Ppargc1a compared to 189 male iWATs (Fig. 3D). Since mgWATs make up almost all of the subcutaneous fat mass in female 190 mice, we reasoned that highly reduced adipocyte thermogenic gene expression could potentially 191 influence whole body energy metabolism. To test this, we performed indirect calorimetry on age-192 matched male and female mice at RT and 24 hr cold exposure using a metabolic chamber. We 193 found that female mice showed highly reduced energy expenditure (EE), oxygen consumption 194 (VO2), and carbon dioxide production (VCO2) during cold exposure compared to male mice ( test the role of ductal epithelium in adipose thermogenesis 37 . 5-week-old female mice were 213 exposed to cold and inguinal and dorsolumbar regions were dissected to assess thermogenic 214 transcripts. Epcam transcripts were present only in the inguinal region, however, thermogenic 215 genes were mostly similar between inguinal and dorsolumbar region except for Ucp1 where we 216 saw an increase in Ucp1 levels in inguinal part (Fig. 3H). Chi et. al reported regional differences 217 between inguinal and dorsolumbar region and there could be a regional control of Ucp1 expression 218 in 5-week-old mice independent of ductal cells 38 . However, consistent with our data, we noticed significantly less expression of genes involved in lipid mobilization such as Pnpla2 and Hsl, in 220 inguinal regions compared to dorsolumbar, indicating that the ductal epithelium is potentially 221 inhibiting cold-induced lipid mobilization. All these data point, for the first time, toward a possible 222 unique mechanism of SNS-activated mgWAT under cold stress to limit thermogenesis and 223 preserve adiposity. 224 Our results show a unique possible SNS-mediated crosstalk between mammary ductal cells 225 and adipocytes to control adipocyte thermogenesis. Our differentially expressed genes in 226 mammary ductal cells under adrenergic stimulation showed upregulation of genes that encode 227 factors such as Angptl4, Enho, Lrg1, and Lcn2, which are known to play inhibitory roles in 228 adipocyte thermogenesis. These factors also showed high enrichment in mammary Epcam+ cells 229 by scRNA-seq and RNAscope in situ hybridization. Among these factors, Lcn2 was previously 230 shown to inversely correlate with Ucp1 expression in female gonadal WAT (gWAT) 24 . Female 231 gWATs, which lack mammary ducts, express exceedingly high levels of Ucp1 compared to male 232 gWATs 24 . Reciprocally, we find that cold-stressed female mgWATs express high levels of Lcn2 233 and low levels of Ucp1 compared to males ( Fig. 3D and Extended Data Fig. 4A). The secretion of 234 Lcn2 by luminal AV and HS-AV cells could potentially be a mechanism of luminal cells to block 235 excess thermogenesis and preserve adiposity. Consistent with data in Fig.1  Lcn2 under cold stress, female Lcn2 KO and WT mice were exposed to cold (4°C) for 24 hr. 260 Compared to controls, cold-exposed Lcn2 KO mice mgWAT showed more beiging/browning and 261 gene expression analysis showed a significant increase in thermogenic genes such as Ucp1 262 indicating that Lcn2 is potentially one of the limiting factors involved in regulating the propensity 263 of mgWAT to beige (Extended Data Fig. 4J). Overall, our scRNA-seq analysis and both our tissue-264 specific gain-of-function and loss-of-function experimental data show Lcn2 as a factor expressed 265 in luminal cells that could function to inhibit thermogenesis and maintain adiposity in mgWAT 266 during cold-stress. 267 In summary, our studies have uncovered a direct inhibitory role of mammary duct luminal 268 cells in adipocyte Ucp1 expression and thermogenic gene program. SNS fibers directly innervate 269 The Krt8rtTA-TetO-DTA mouse model was described previously 39 . For the time course cold 285 exposure experiment, WT mice at 8-10 weeks of age were singly housed at 4°C room in a non-286 bedded cage without food and water for first 6 h; thereafter food, water, and one cotton square 287 were added. For the 24 h harvest, 3 h before harvest, food, water, and cotton square were removed 288 and then mice were harvested. At the end of the experiment, mgWATs were resected for analysis. 289 For overexpression studies, recombinant adeno-associated virus serotype 8 (AAV8) expressing 290 LCN2 or GFP was generated and injected as described previously 24 . Indirect calorimetry was 291 performed using Promethion Systems (Sable Systems, Las Vegas, NV). Animals were placed 292 individually in chambers at ambient temperature (22.0 °C) for 21 hr followed by 24 hr cold (4.0 293 °C) with 12 hr light/dark cycles. Animals had free access to food and water. Respiratory 294 measurements were made in 5 min intervals after initial 7-9 hr acclimation period. Single cell SVF populations from adipose tissue were isolated as described previously 3,31 . The 316 fourth inguinal white adipose tissue (iWAT) depot mgWAT from mice exposed to cold stress 317 (4°C) or room temperature for 24 hr was dissected and placed on a sterile 6-well tissue culture 318 plate with ice-cold 1X DPBS. Excess liquid was removed from fat pads by blotting. Each tissue 319 was cut and minced with scissors and then placed in 15 ml conical tubes containing digestion 320 buffer (2 ml DPBS and Collagenase II at 3 mg/ml; Worthington Biochemical, Lakewood, NJ, 321 USA) for 40 min of incubation at 37°C with gentle shaking at 100 rpm. Following tissue digestion, 322 enzyme activity was stopped with 8 ml of resuspension media (DMEM/F12 with glutamax 323 supplemented with 15%FBS and 1% pen/strep; Thermo Scientific, CA). The digestion mixture 324 was passed through 100 μm cell strainer and centrifuged at 150 x g for 8 min at room temperature. 325 To remove red blood cells, the pellet was resuspended and incubated in RBC lysis buffer (Thermo 326 Scientific, CA) for 3 min at room temperature, followed by centrifugation at 150 x g for 8 min. 327 The pellet was then resuspended in resuspension media and again spun down at 150 x g for 8 min. Stem Cell Research Center sequencing core were demultiplexed and converted to fastq format 352 using the 'mkfastq' function from Cell Ranger. Next, the Cell Ranger 'counts' function mapped 353 reads from fastq files to the mm10 reference genome and tagged mapped reads as either exonic, 354 intronic, or intergenic. Only reads which aligned to exonic regions were used in the resulting 355 DGEs. After combining all four sample DGEs into a single study DGE, we filtered out cells with 356 (1) UMI counts < 700 or > 30,000, (2) gene counts < 200 or > 8,000, and (3) mitochondrial gene 357 ratio > 10%. This filtering resulted in a dataset consisting of 42,052 genes across 12,222 cells, with 358 approximately 2,300 -4,650 cells from each sample. A median of 2,411 genes and 7,252 359 transcripts were detected per cell. 360 361

Identification of cell clusters 362
To achieve high resolution cell type identification and increased confidence in our cell type 363 clustering we brought in external publicly available single cell data from SVF and mammary 364 tissues. Specifically, we included single cell data from 9 datasets comprising 91,577 single cells 365 from the mammary gland and multiple adipose depots, across 4 different single cell platforms 366 (Table 1). These external datasets and the SVF data from this study were all independently 367 normalized using sctransform 50 and integrated using Seurat 51,52 v3.1.5. The single cell expression 368 profiles were projected into two dimensions using UMAP 53 or tSNE 54 and the Louvain 55 method 369 for community detection was used to assign clusters. This integrated data was only used to identify 370 and define the cell types. All plots which are not explicitly designated as integrated with at least 371 one external dataset and all downstream analyses (e.g. differential expression analyses) were 372 conducted on non-integrated data to retain the biological effect of the cold treatment. Visualization 373 of the non-integrated data was conducted on a subsampled dataset where all samples had the same 374 number of cells to give an equal weight to each sample, however, all downstream analyses (e.g. 375 differential expression analyses) leveraged the full dataset.   Cell type-specific gene expression signatures were generated by identifying genes with expression 384 levels two-fold greater (adjusted p-values < 0.05) than all other cell types. To ensure consistency 385 across samples, Seurat's FindConservedMarkers function (Wilcoxon rank sum test with a meta p-386 value) was applied across each sample. 387

Resolving cell identities of the cell clusters 389
To identify the cell type identity of each cluster, we used a curated set of canonical marker genes 390 derived from the literature (Supplementary Table 1) to find distinct expression patterns in the 391 cell clusters. Clusters which uniquely expressed known marker genes were used as evidence to 392 identify that cell type. Cell subtypes which did not express previously established markers were 393 identified by general cell type markers and novel markers obtained with Seurat's 394 FindConservedMarkers function were used to define the cell subtype. 395 396

Differential gene expression analysis 397
Within each identified cell type, cold treated and room temperature single cells were compared for 398 differential gene expression using Seurat's FindMarkers function (Wilcoxon rank sum test) in a 399 manner similar to Li et al. 19 . Differentially expressed genes were identified using two criteria: (i) 400 an expression difference of >= 1.5-fold and adjusted p-value < 0.05 in a grouped analysis between 401 room temperature mice (n = 2) and cold treated mice (n = 2); (ii) an expression difference of >= were determined by using a standard curve. Each gene was normalized to the housekeeping gene 418 36B4 and was analyzed in duplicate. Primers used for real-time PCR are previously described 3,31 419 and presented in Table 2. 420 421

RNAScope Fluorescence in situ hybridization (FISH) 422
mgWAT from RT or cold exposed mice (Jackson Laboratory, #000664) was fixed in 10% formalin 423 overnight, embedded with paraffin, and sectioned into unstained, 5 μm-thick sections. Sections 424 were baked at 60°C for 1 hour, deparaffinized, and baked again at 60°C for another hour prior to 425 pre-treatment. The standard pre-treatment protocol was followed for all sectioned tissues. In situ 426 hybridization was performed according to manufacturer's instructions using the RNAscope

Sample collection 496
Immediately after cold exposure mice were anaesthetized with isoflurane (3%) and perfused with 497 heparinized saline followed by 4% paraformaldefyde (PFA) (Electron Microscopy Sciences, 498 Hatfield, PA, USA). Fat pads were carefully dissected and postfixed overnight in 4% PFA at 4 ᵒ C. 499 On the following day the tissue was washed 3 times in PBS before proceeding with the optical 500 clearing protocol. 501

Optical clearing 503
Whole-mount staining and clearing was performed using the Adipo-Clear protocol as previously 504 with 0.5% Trixon X-100, 0.1% Tween-20, 2 µg/ml heparin, as previously described 66 and applied 512 for 6 days at 37 ᵒ C. Following 5 washes with modified PTxwH over 1 day with the last wash 513 performed overnight, secondary antibodies were diluted in modified PTxwH (1:500) and samples 514 incubated for 6 days at 37 ᵒ C. Samples were washed 5 times over 1 day in modified PTxwH at 37 ᵒ C, 515 5 times over 1 day in PBS at RT, embedded in 1% agarose, dehydrated with a methanol gradient 516 in H2O (12%, 50%, 75%, 100%), washed 3 times for 1 hr in 100% methanol followed by 3 times 517 for 1 hr in DCM, before being transferred to dibenzylether (DBE) (Sigma-Aldrich) to clear. 518 Imaging 519 Z-stacked optical sections of whole fat pads were acquired with an Ultramicroscope II (LaVision 520 BioTec, Bielefeld, Germany) at a 1.3x magnification with a 4 µm step size and dynamic focus 521 with a maximum projection filter. Samples were then imaged in glass-bottom μ-dishes (81158, 522 Ibidi, Gräfelfing, Germany) using an inverted Zeiss LSM 710 confocal microscope with a 10x 523 (NA: 0.3) objective and a step size of 5 µm. 524

Image analysis 526
Imaris versions 9.6.0-9.7.2 (Bitplane AG, Zürich, Switzerland) were used to create digital surfaces 527 covering ducts, TH+ innervation and total sample volume (1.3x light sheet images and 10x 528 confocal images) to automatically determine volumes and intensity data. Volume reconstructions 529 were performed using the surface function with local contrast background subtraction. For 530 detection of EPCAM+ ducts in 1.3x light sheet images, a smoothing factor of 5 µm was used and 531 the threshold factor was set to correspond to the largest duct diameter in each sample. For detection 532 of TH+ nerves in 1.3x light sheet images, a smoothing factor of 3 µm and a threshold factor of 80 533 µm were used. For detection of EPCAM+ ducts in 10x confocal images, a smoothing factor of 534 3.35 µm was used and a threshold factor corresponding to the diameter of the thickest duct wall in 535 each sample was used. For detection of TH+ nerves in 10x confocal images, a smoothing factor of 536 2 µm and a threshold factor of 5 µm were used. In 10x confocal images, nerve/duct interactions 537 were defined by masking the TH channel using the TH+ nerve surface to remove any background, 538 and then masking it again using the EPCAM+ duct surface. This process revealed TH+ innervation 539 overlapping with EPCAM+ staining. A new surface was created to cover this overlapping TH+ 540 innervation using a smoothing factor of 2 µm and a threshold factor of 5 µm. 541

Statistical analyses 542
Data are shown as mean±S.E.M. Distribution was assessed by Shapiro-Wilk test. Significance was 543 determined by a two-tailed unpaired t test (parametric distribution) or by a Mann-Whitney test 544 (non-parametric distribution). Significance was set at an alpha level of 0.05. 545

RNA extraction and quantitative real-time PCR in organoid samples 566
To perform RNA extraction, isolated organoids were collected into kit lysis buffer. RNA was 567 extracted with Qiagen RNeasy Micro Kit. After nanodrop RNA quantification and analysis of 568 RNA integrity, purified RNA was used to synthesize the first-strand cDNA in a 30 μl final volume, 569 using Superscript II (Invitrogen) and random hexamers (Roche). Genomic contamination was 570 detected by performing the same procedure without reverse transcriptase. Quantitative PCR analyses were performed with 1 ng of cDNA as template, using FastStart Essential DNA green 572 master (Roche) and a Light Cycler 96 (Roche) for real-time PCR system. Relative quantitative 573 RNA was normalized using the housekeeping gene Gapdh. Analysis of the results was performed 574 using Light Cycler 96 software (Roche) and relative quantification was performed using the ddCt 575 method using Gapdh as a reference.