Cold-associated mammokines preserve adipocyte identity

Almost four decades of research suggest a dynamic role of ductal epithelial cells in adipocyte adaptation in mammary gland white adipose tissue (mgWAT), but factors that mediate such communication are not known. Here, we identify a complex intercellular crosstalk in mgWAT revealed by single-cell RNA-seq (scRNA-seq) and comprehensive data analysis suggest that epithelial luminal cells during cold exposure undergo major transcriptomic changes that lead to the expression of an array of genes that encode for secreted factors involved in adipose metabolism such as Adropin (Enho), neuregulin 4 (Nrg4), angiopoietin-like 4 (Angptl4), lipocalin 2 (Lcn2), milk fat globule-EGF factor 8 (Mfge8), Insulin-like growth factor-binding protein 1 (Igfbp1), and haptoglobin (Hp). To define the mammary epithelial secretome, we coin the phrase “mammokines”. We validated our cluster annotations and cluster-specific transcriptomics using eight different adipose scRNA-seq datasets including Tabula Muris and Tabula Muris Senis. In situ mRNA hybridization and ex vivo isolated mgWAT luminal cells show highly localized expression of mammokines in mammary ducts. Trajectory inference demonstrates that cold-exposed luminal cells have similar transcriptional profiles to lactation post-involution (PI), a phase defined by reappearance and maintenance of adipocytes in the mammary gland. Concomitantly, we found that under cold exposure female mgWAT maintains more adipogenic and less thermogenic potential than male scWAT. Conditioned media from isolated mammary epithelial cells treated with isoproterenol suppressed thermogenesis in differentiated beige/brown adipocytes and treatment of beige/brown differentiated adipocyte with LCN2 suppresses thermogenesis and increases adipogenesis. Finally, we find that mice lacking LCN2 show markedly higher colddependent thermogenesis in mgWAT than controls, and reconstitution of LCN2 in the mgWAT of LCN2 knockout mice promotes inhibition of thermogenesis. These results show a previously unknown role of mammary epithelium in adipocyte metabolism and suggest a potentially redundant evolutionary role of mammokines in maintaining mgWAT adiposity during cold exposure. Our data highlight mammary gland epithelium as a highly active metabolic cell type and mammokines could have broader implications in mammary gland physiology and lipid metabolism.

differences, and iv) mammary gland development (Fig. 1B). This integrated data set allowed us to precisely annotate various cell types present in mgWAT of female mice. Further subclustering analysis based on known cell type marker genes from our integrated approach identified clusters luminal-HS-AV, and myoepithelial cells ( Fig. 1C and 1D).
To gain insight into the remodeling of stromal cells under adrenergic stress, we segregated the cumulative tSNE-plot into RT and COLD treatment by animal replicate (Fig. 1E). The tSNE and dot plots reveal global changes in relative proportions of SVF clusters between RT and cold ( Fig. 1E and 1F). Among all the clusters APCs, macrophages, and luminal cells showed marked differences in cell type percentages (Fig. 1F) and appeared to have large differences in their global transcriptomic profiles in the t-SNE two-dimensional projection where cells from RT and cold were segregated (Fig. 1E). To quantitatively determine the transcriptional impact of cold treatment on individual cell types, we plotted fold change and adjusted p-values as a function of cluster types and found a high degree of transcriptional variation in APCs, macrophages, and luminal HS and AV under cold conditions (Fig. 1G). For the most varied cell types (APC, macrophage, and luminal cells), we subclustered these cell types and found two subclusters of APCs (APC1 and APC2) 16 Fig. 1B, 1G, and 1H). The tSNE plots highlighted global transcriptional differences between RT and cold in subclusters (Extended Data Fig. 1B and 1H) and highly expressed genes in those subclusters (Extended Data Fig. 1C and 1I). We then analyzed differentially expressed genes (DEGs) in APC, luminal, and macrophage clusters which were represented as a volcano plot with the ratio of fold change between cold and RT as a function of adjusted p-value. APCs showed an increase in expression of adipogenic genes such as Cebpd and Klf9 17-20 and a decrease in genes encoding for collagens suggesting and supporting the previous observation of adipogenesis under cold exposure (Extended Data Fig. 1D and 1E) 16 . Macrophages on the other hand namely CD14 macrophages and FN1 macrophages showed marked reduction in relative cell fraction and volcano plots demonstrated a decrease in Tlr2, Rel and Nfkb singling in Fn1 macrophages and increase in the expression of Apoe2, Plin2 and Ddit4 in Cd14 and Fn1 macrophage (Extended Data Fig. 1J and 1K).

Identification of mammokines in mgWAT
Luminal cells showed a remarkable transcriptional difference in cell clusters between RT and cold, implicating a potential remodeling of the luminal epithelium upon cold exposure ( Fig.   1G-1I). To gain insight into the crosstalk between cell types within the mgWAT SVFs at RT and cold, we applied a bioinformatics analysis, CellPhoneDB 21 , an established method for inferring cell-cell communication from single cell RNA-seq. CellPhoneDB leverages a vast repository of curated receptors, ligands and their interactions to infer cell-cell communication. We focused our analysis on the communication between APCs and previously reported 24h cold-exposed mgWAT mature adipocytes single nuclei RNA-seq data (SNAP-seq) 22 with luminal epithelial cells within mgWAT (Extended Data Fig. 2). We noticed a prominent crosstalk between APCs and luminal cells via Wnt, Tgf, Tnf, BMPs, Ephrin signaling. Even though overall transcriptomic dynamics of luminal cells were increased, however, comparing with RT, cold stressed mice showed an overall decrease in the communication between the APCs and luminal cells (Extended Data Fig. 2). To probe for cold-induced ligand-receptor interaction between APCs and luminal cells, we then focused our CellPhoneDB analysis on these cell types and found increases in NRGs, Wnt7B, BDNF and FGF1 signaling. These data indicate that under cold temperature luminal cells could elicit adipogenic signals to APCs by factors such as FGF1 and Wnts (Extended Data Fig. 2B) that could lead to an increase in adipogenic transcription factors such as CEBPD and KLF9 23 . To test crosstalk between mature adipocytes and luminal cells, we analyzed the data reported in this study as a function of our published SNAP-seq 22 data, which is the first reported single nuclei RNA-seq data for adipocytes in mgWAT in females. Based on our analysis, we find that luminal cells have a distinct ligand-receptor interaction with mature adipocytes which are distinct to those observed in APCs (Extended Data Fig. 2C and 2D). We had previously reported 14 adipocyte subclusters in mgWAT 22 and among them we identified Cluster 12/14 as hormonally active adipocyte cluster, cluster 4 and 6 as traditional lipid storing adipocyte clusters, and cluster 9 as thermogenic adipocyte cluster. Our integrated CellPhoneDB data now show cold-induced signaling between luminal cells and i) Cluster 12/14 cells via receptor KIT ligand (KITLG) 24,25 , ii) Cluster 4 and 6 via Wnt4 26 , and iii) cluster 9 by inhibitin beta B (Inhbb) 27 (Extended Data Fig. 2C and 2D). Adipocyte are important for ductal morphogenesis and development 7 and our Cellphone DB data point toward a coordinated efforts of different kinds of mgWAT adipocytes in inducing luminal epithelium development. Besides these ligands, we also find an increase in VEGF, TGF, NOTCH, AREG, and, NECTIN4 signaling between adipocytes and luminal cells.
Some of the mammokines are known to be expressed by other cell types; however, we find that cold-induced mammokines and cold-induced luminal epithelial genes (ci-LEG) have distinct expression pattern and are highly enriched in EPCAM+ve epithelial cluster of mgWAT.
Furthermore, Adrb2 and Adrb1 expression in luminal cells, suggest that these cells may directly respond to cold-induced SNS activation [38][39][40] (Fig. 1K). We next assessed the expression of mammokines and ci-LEG across different adipose depots in young and aged male and female mice.
We queried the expression profiles of mammokines and ci-LEG in the Tabula Muris and Tabula Senis data set that include scRNA-seq from the SVFs of gonadal adipose tissue-GAT, mesenteric adipose tissue-MAT, and subcutaneous adipose tissue-SCAT, and mgWAT. The Tabula Muris scRNA-seq annotations from mgWAT showed five distinct cell clusters that include basal myoepithelial cells, two epithelial cells (HS and AV), endothelial cells , and stromal cells.
(Extended Data Fig. 3A). Mammokines and ci-LEG showed mostly luminal epithelial cells enrichment and their expression did not vastly differ across various ages (Extended Data Fig. 3A-3C). SVF scRNA-seq from GAT, MAT, SCAT from males and female mice showed specificity of mammokines and ci-LEG expression in female SCAT. However, genes such as Angptl4, Glul, and F3 were also expressed in APCs, endothelial cells, and myeloid cell populations (Extended Data Fig. 3D-3I). t-SNE plots of normalized gene expression levels for cold-induced mammokines in mgWAT (our study), male scWAT SVFs and mature adipocytes 22 show that Angptl4 and PRLR are also expressed by most mature adipocytes; however, other cold-induced mammokine genes showed relatively localized expression in luminal epithelial cells (Extended Data Fig. 3J) Finally, to infer cold-induced transcription factor gene regulatory networks across all our scRNA-seq data sets we used SCENIC 41 (single-cell regulatory network inference and clustering) analysis to derive top 5 transcription factors enriched for the transcriptional state of each cell type.
SCENIC revealed highly dynamic cold-induced transcription factor activity across all cell types and we found increased activity of COX14, IRX2, ELF5 ZFP708B, and TfFAP2A in luminal AV cells and FOXA1, BHLHE41, TEAD2, TFAP2B, and GRHL1 in Luminal HS cells ( Fig. 1L and Extended data 4A-4C). Regulon activity implicated BHLHE41 as a master regulator of genes in luminal HS cells, and ELF5 showed a cold-specific gene transcriptional activity and could be directly involved in the expression of mammokines and ci-LEG such as Lcn2 and Igfbp5 (Extended Data Fig. 4D). ELF5 42 is an important transcriptional regulator of mammary epithelial development and our data indicate ELF5 as an adrenergic-activated master regulator in luminal cells that could enhance expression of mammokines under cold stress.

Cold-induced mgWAT luminal cells resemble lactation post-involution transcriptomic state
In agreement with the data in Extended Fig. 3, our RNAscope fluorescent in situ hybridization (FISH) analysis showed a highly localized expression of mammokines Enho, Lrg1, Lcn2, Hp, Nrg4, and Wnt4 in Epcam+ve mammary ductal cells ( Fig. 2A and Extended Data Fig. 5). To probe for a possible role of cold induced mammokines and ci-LEG in mgWAT physiology, we next performed trajectory inference using Slingshot 43 from mammary gland EPCAM+ve cells at nulliparous (NP), gestation (G), lactation (L), and lactation post-involution (PI) 12 (Fig. 2B) . Our analysis showed distinct nodal points for different mgWAT states (NP, G, L, and PI) and coclustering of RT and cold luminal HS, AV, and HS-AV cells ( Data 6B and 6C and 7B and 7C). PI is a stage in the mammary gland where adipocytes reappear and fill the fat pad 8 . Taken together, our CellPhoneDB analysis, SCENIC, DEG analysis, and trajectory inference suggest that mgWAT luminal cells of female mice during cold-exposure could produce mammokines that may act to promote adiposity in mgWAT by facilitating adipogenesis and adipocyte expansion, and blockade of thermogenesis. Our integrated analyses further supports the data in Extended Data Fig. 1A, where we report a markedly less propensity of females mgWAT to beige and higher capacity of mgWAT to maintain adiposity following cold exposure.

Cold-induced mammokines inhibit adrenergic-dependent adipose thermogenesis
To experimentally validate our analysis and increase in EPCAM+ve luminal cell population (Fig.   3A), we first purified EPCAM+ve and EPCAM-ve cells from mgWAT of RT and cold-exposed mice and then stained for EPCAM and CD49F for FACS analysis. As shown in Fig. 3B, we saw three distinct populations of cells that were EPCAM lo CD49F lo (stromal), EPCAM hi CD49F hi (luminal), and EPCAM +ve CD49F hi (basal) cells. Luminal cells were only enriched in EPCAM+ve purified cells and showed a significant cold-dependent increase in cell population compared to RT condition, in agreement with our scRNA-seq data. We also noticed a marked reduction in the basal cell population upon cold treatment in EPCAM+ve selected cells, indicating a differential response to cold stress by luminal and basal cell in ducal epithelium. To test whether luminal cells directly responded to cold-induced SNS activation, we tested mammokine expression on isolated primary mgWAT EPCAM+ve and EPCAM-ve cells from RT or cold mice or cultured primary mgWAT EPCAM+ve and EPCAM-ve cells that were either treated or not treated with isoproterenol. Cold-induced mammokines were highly enriched in EPCAM+ve cells and showed increased expression upon isoproterenol or cold treatment (Fig 3C-3E). To test if adrenergicdependent secretion from mammary epithelial cells blocked adipocyte thermogenesis, we treated normal mammary epithelial cell NMuMG with varying concentration of isoproterenol and then treated brown/beige differentiated 10T1/2 cells with the conditioned media from NMuMG cells.
Consistent with the data in Fig. 3D, we find a dose-dependent increase in cold-induced mammokines and reciprocal decrease in thermogenic genes in conditioned media treated beige/brown differentiated 10T1/2 cells ( Fig. 3F and 3G). We further directly tested the roles of mammokines LRG1 and LCN2 by treating beige differentiated preadipocytes with recombinant LRG1 or LCN2. As shown in Fig 3H and 3I, LRG1 and LCN2 caused a significant reduction in thermogenic gene expression.

LCN2 is a mammokine involved in maintaining mgWAT adiposity
Overall, our broad scRNA-seq data and targeted experimental data consistently identify LCN2 as a mammokine that is i) highly enriched in mammary luminal epithelial cells, ii) induced by both cold and isoproterenol treatment, iii) expressed mostly in the PI phase of mammary gland, iv) and potentially regulated by the activity of ELF5 in luminal cells under cold stress. The secretion of LCN2 by luminal AV and HS-AV cells could potentially be a mechanism of luminal cells to block excess lipid mobilization, thermogenesis, and preserve adiposity. Consistent with data in Fig.1 and showed that LCN2 exogenous expression causes a significant decrease in thermogenic genes such as Ucp1, Cidea, Ppara and increase in adipogenic genes such as Lep, Mmp12 45,46 , Zfp423 47 , and Lbp 48 . LCN2 expression also led to an increase in Aldh1a1, which was recently shown to inhibit adipose thermogenesis by downregulating UCP1 levels 49,50 . We validated our RNA-seq data by directed qPCR in both LCN2 reconstituted mgWAT of LCN2KO and WT mice (Fig. 4C). Finally, to test the role of LCN2 under cold stress, we exposed female LCN2KO and WT mice to 24 h cold exposure. LCN2KO mice mgWAT showed more beiging/browning and gene expression analysis showed a significant increase in thermogenic genes such as Ucp1 and decreases in adipogenic genes such as Lep indicating that LCN2 is potentially one of the limiting factors involved in regulating the propensity of mgWAT to beige ( Fig. 4D and 4E). However, paradoxically, we noticed that cold exposed LCN2KO mice mgWAT showed increase in de novo lipogenic genes, which could suggest a potential negative feedback loop to compensate excess for lipid mobilization 51 . Compared to males, LCN2KO females had significantly decreased in body and mgWAT weight (Extended Fig. 10A and 10B). We noticed both sexes had an increase in plasma free glycerol and decrease in free fatty acid (FFA), implicating a potential local effect of mammokines on the local interaction between epithelium and adipocyte metabolism. However, plasma triglyceride levels were significantly lower in male LCN2KO compared to controls whereas TG levels were comparable between LCN2KO and WT mice (Extended Data Fig. 10C-10E). Overall, our scRNAseq analysis and both our tissue-specific gain-of-function and loss-of-function experimental data show LCN2 as a factor expressed in luminal cells that function to maintain adiposity in mgWAT during cold-stress.          and surrounding SVF at room temperature were identified using SCENIC.

(B)
The top five transcription factor regulons for all cell types in the mammary gland (this study) and surrounding SVF at 4 degrees were identified using SCENIC. Genes in red were derived from the room temperature data, genes in blue were derived from the 4-degree data, and genes with both were present in the regulons in both conditions. Red arrow show cold-induced mammokines.

Figure 5
(A-E) RNAscope FISH for indicated probes from RT or cold exposed female mice mgWAT.     C-E) Plasma free glycerol, free fatty acid (FFA), and triglyceride (TG) levels from male and female at RT or 24 cold-exposed.

Animal Studies
C57BL/6 WT males and female mice (#000664) and LCN2KO (#24630) were acquired from Jackson Laboratory and maintained in a pathogen-free barrier-protected environment (12:12 h light/dark cycle, 22°C-24°C) at the UCLA and Mount Sinai animal facility. For the time course cold exposure experiment, WT mice at 8-10 weeks of age were singly housed at 4C room in a non-bedded cage without food and water for first 6 h; thereafter food, water, and one cotton square were added. For the 24 h harvest, 3 h before harvest, food, water, and cotton square were removed and then mice were harvested. At the end of the experiment, mgWATs were resected for analysis. For overexpression studies, recombinant adeno-associated virus serotype 8 (AAV8) expressing LCN2 or GFP was generated and injected as described previously 31 . Animal experiments were conducted in accordance with the Mount Sinai and UCLA Institutional Animal Care and Research Advisory Committee.

RNA-Seq
RNA isolation, library preparation, and analysis were conducted as previously described 31  Diego, CA). Reads were aligned to the mouse genome mm10 using STAR 52 or HISAT2 53 aligner and quantified using the Bioconductor R packages as described in the RNA-Seq workflow 54 . P values were adjusted using the Benjamini-Hochberg procedure of multiple hypothesis testing 54 .

Single cell isolation from SVF
Single cell SVF populations from adipose tissue were isolated as described previously 22,55 . The fourth inguinal white adipose tissue (iWAT) depot mgWAT from mice exposed to cold stress (4°C) or room temperature for 24 hr was dissected and placed on a sterile 6-well tissue culture plate with ice-cold 1X DPBS. Excess liquid was removed from fat pads by blotting. Each tissue was cut and minced with scissors and then placed in 15 ml conical tubes containing digestion buffer (2 ml DPBS and Collagenase II at 3 mg/ml; Worthington Biochemical, Lakewood, NJ, USA) for 40 min of incubation at 37°C with gentle shaking at 100 rpm. Following tissue digestion, enzyme activity was stopped with 8 ml of resuspension media (DMEM/F12 with glutamax supplemented with 15%FBS and 1% pen/strep; Thermo Scientific, CA). The digestion mixture was passed through 100 μm cell strainer and centrifuged at 150 x g for 8 min at room temperature.
To remove red blood cells, the pellet was resuspended and incubated in RBC lysis buffer (Thermo Scientific, CA) for 3 min at room temperature, followed by centrifugation at 150 x g for 8 min.
The pellet was then resuspended in resuspension media and again spun down at 150 x g for 8 min.
The cell pellet was resuspended in 1 ml of 0.01% BSA (in DPBS) and passed through a 40 μm cell strainer (Fisher Scientific, Hampton, NH, USA) to discard debris. Cell number was counted for 10X Genomics single cell application.

SVF single cell barcoding and library preparation
To yield an expected recovery of 4000-7000 single cells, an estimated 10,000 single cells per channel were loaded onto Single Cell 3' Chip (10X Genomics, CA). The Single Cell 3' Chip was placed on a 10X Genomics instrument to generate single cell gel beads in emulsion (GEMs). Chromium Single Cell 3' v3 Library and Cell Bead Kits were used according to the manufacturer's instructions to prepare single cell RNA-Seq libraries.

Illumina high-throughput sequencing libraries
Qubit Fluorometric Quantitation (ThermoFisher, Canoga Park, CA, USA) was used to quantify the 10X Genomics library molar and a TapeStation (Aligent, Santa Clara, CA, USA) was used to estimated library fragment length. Libraires were pooled and sequenced on an Illumina HiSeq 4000 (Illumina, San Diego, CA, USA) with PE100 reads and an 8 bp index read for multiplexing.
Read 1 contained the cell barcode and UMI and read 2 contained the single cell transcripts.

Single cell data pre-processing and quality control
To obtain digital gene expression matrices (DGEs) in sparse matrix representation, paired end reads from the Illumina HiSeq 4000 were processed and mapped to the mm10 mouse genome using 10X Genomics' Cell Ranger v3.0.2 software suite. Briefly, .bcl files from the UCLA Broad Stem Cell Research Center sequencing core were demultiplexed and converted to fastq format using the 'mkfastq' function from Cell Ranger. Next, the Cell Ranger 'counts' function mapped reads from fastq files to the mm10 reference genome and tagged mapped reads as either exonic, intronic, or intergenic. Only reads which aligned to exonic regions were used in the resulting

Identification of cell clusters
To achieve high resolution cell type identification and increased confidence in our cell type clustering we brought in external publicly available single cell data from SVF and mammary tissues. Specifically, we included single cell data from 9 datasets comprising 91,577 single cells from the mammary gland and multiple adipose depots, across 4 different single cell platforms ( Table 1). These external datasets and the SVF data from this study were all independently normalized using sctransform 56 and integrated using Seurat 57,58 v3.1.5. The single cell expression profiles were projected into two dimensions using UMAP 59 or tSNE 60 and the Louvain 61 method for community detection was used to assign clusters. This integrated data was only used to identify and define the cell types. All plots which are not explicitly designated as integrated with at least one external dataset and all downstream analyses (e.g. differential expression analyses) were conducted on non-integrated data to retain the biological effect of the cold treatment. Visualization of the non-integrated data was conducted on a subsampled dataset where all samples had the same number of cells to give an equal weight to each sample, however, all downstream analyses (e.g. differential expression analyses) leveraged the full dataset.

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

Quantitative assessment of global transcriptome shifts
Quantitative assessment of global transcriptome shifts in each cell type due to cold treatment was assessed using a previously established Euclidean distance-based approach 11 . Briefly, average gene expression profiles for each condition (room temperature and cold treated) were generated for each cell type. After representing the gene expression as a z-score and only considering the top 1000 most highly expressed genes in each cell type, the Euclidean distance between room temperature and cold treated conditions was calculated for each cell type. To assign significance to the distances, a null distribution was estimated by calculating the Euclidean distance between two groups of randomly sampled cells of the given cell type. This permutation approach was repeated 1000 times and compared to the Euclidean distance from the true group labels to determine an empirical p-value. This p-value was adjusted for multiple testing with a Bonferroni correction.

Differential gene expression analysis
Within each identified cell type, cold treated and room temperature single cells were compared for differential gene expression using Seurat's FindMarkers function (Wilcoxon rank sum test) in a manner similar to Li et al. 15 . Differentially expressed genes were identified using two criteria: (i) an expression difference of >= 1.5 fold and adjusted p-value < 0.05 in a grouped analysis between room temperature mice (n = 2) and cold treated mice (n = 2); (ii) an expression difference of >= were considered statistically significant.

CellPhoneDB analysis
To predict cell-cell interactions with the single cell data, we used the CellPhoneDB method 21 .
Breifly, CellPhoneDB has a curated database of >2400 interactions in the categories of proteinprotein interactions, secreted and membrane proteins, and protein complexes. This database was used to identify enriched receptor-ligand interactions between two cell types based on the expression of a receptor by one cell type and the corresponding ligand by another cell type.
Significance was obtained using a permutation approach. We applied CellPhoneDB to each sample separately and only considered interactions between cell types significant if it was observed in both samples a particular group (e.g. room temperature or cold treated). Since the curated database of interactions is based on human genes, we first converted our mouse genes to their human orthologues.

Lineage trajectory inference
Prior to lineage trajectory inference, the SVF single cell data from this study was integrated with a dataset demonstrating the differentiation dynamics of mammary epithelial cells 7 using Seurat as described in the "Identification of cell clusters" section. Lineage trajectory inference was performed on this integrated dataset by applying slingshot 43 which was wrapped in a docker container as part of the dynverse 72 .

Gene regulatory network inference
Gene regulatory network inference was performed with pySCENIC 41 following the workflow described by Van de Sande et al 73 . Briefly, starting with counts data, gene modules which are coexpressed with transcription factors were identified with GRNBoost2 74 . Next, candidate regulons were created from transcription factortarget gene interactions and indirect targets were pruned based on motif discovery with cisTarget 41 . Finally, regulon activity was quantified at cellular resolution with AUCell 41 which allowed for the prioritization of regulons for each cell type based on the quantified activity.

Real time qPCR
Total RNA was isolated using TRIzol reagent (Invitrogen) and reverse transcribed with the iScript cDNA synthesis kit (Biorad). cDNA was quantified by real-time PCR using SYBR Green Master Mix (Diagenode) on a QuantStudio 6 instrument (Themo Scientific, CA). Gene expression levels were determined by using a standard curve. Each gene was normalized to the housekeeping gene 36B4 and was analyzed in duplicate. Primers used for real-time PCR are available upon request.
Sections were baked at 60°C for 1 hour, deparaffinized, and baked again at 60°C for another hour prior to pre-treatment. The standard pre-treatment protocol was followed for all sectioned tissues. In situ hybridization was performed according to manufacturer's instructions using the

CONFLICTS OF INTEREST
The authors declare no conflicts of interest.