Single-cell transcriptomic analysis reveals dynamic activation of cellular signaling pathways regulating beige adipogenesis

AnimalsAll animal protocols were approved by the Institutional Animal Care and Use Committees at Seoul National University (SNU-200302-3, SNU-191113-2, and SNU-191022-2). Mice were fed a standard chow diet, housed at 22 °C, and maintained on a 12-h light/12-h dark cycle with free access to food and water at all times. Male mice were used for the experiments. C57BL/6 mice (6 weeks old, male) were purchased from the Central Lab. Animal, Inc (Seoul, Korea). PDGFRA-Cre [stock #013148, C57BL/6-Tg(PDGFRA-Cre)1Clc/J], PDGFRA-CreER [stock #018280, B6N.Cg-Tg(PDGFRA-Cre/ERT) 467Dbe/J], ROSA26-LSL-ZsGreen [stock #007906, B6.Cg-Gt(ROSA)26Sortm6(CAG-ZsGreen1)Hze/J], and ROSA26-LSL-tdTomato [stock#007909, B6.Cg-Gt(ROSA)26Sortm9(CAG-tdTomato)Hze/J] mice were purchased from the Jackson Laboratory (Bar Harbor, ME, USA). DPP4-CreER [stock#RMRC13230] mice were purchased from NARLabs (National Laboratory Animal Center, Taipei, Taiwan). DPP4-CreER, PDGFRA-Cre, and PDGFRA-CreER mice were crossed with ROSA26-LSL-tdTomato or ZsGreen mice to trace DPP4+ cells and PDGFRA+ cells during CL treatment or for transplantation assays. For Cre recombination, double transgenic mice and wild-type controls were treated with tamoxifen dissolved in sunflower oil (75 mg/kg; Sigma‒Aldrich, St. Louis, MO, USA) by oral gavage every 5 consecutive days3. Additionally, we included vehicle (oil)-treated control groups to confirm the specific induction of Cre recombinase activity in adipose tissues upon tamoxifen treatment. Experiments were started 10 days after the last dose of tamoxifen. For beta3 adrenergic receptor stimulation, mice were treated with CL316243 (1 mg/kg/day; Sigma‒Aldrich) by intraperitoneal injection for up to 7 days. For EdU labeling of proliferating cells, mice were injected with EdU (0.4 μmol/mouse; Thermo Fisher, Waltham, MA, USA) at the indicated times. For the transplantation assay, PDGFRA+ DPP4+ and PDGFRA+ DPP4− cells were isolated from PDGFRA-tdTomato mice and subcutaneously injected with Matrigel (Corning Inc., Corning, NY, USA), as described previously3. Matrigel transplants were collected and analyzed 3 weeks after injection.Cell culturesC3H10T1/2 cells (ATCC, Manassas, VA, USA) were cultured as previously described11. For genetic modulation using siRNA, cells were transfected with 20 nM predesigned siRNA targeting Gli1 (Bioneer, Daejeon, Korea, 14632) or Rbpj (Bioneer, 19664) or negative control siRNA (Bioneer, SN-1013) using INTERFERin® transfection reagent (Polyplus, Illkirch-Graffenstaden, France, #409-10) according to the manufacturer’s instructions. For overexpression of GLI1 and the GFP-tagged GLI2 protein, pLUT7-HA-GLI1 (Addgene, Watertown, MA, USA, plasmid #62967) and the pCEFLmGFP-Gli2 vector (Addgene, plasmid #37672) were obtained and transfected into the cells by using a JetPrime Transfection Kit (Polyplus, #114-15) according to the manufacturer’s instructions. For inhibition of DPP4 activity, sitagliptin (Sigma‒Aldrich) was used.WAT stromovascular cell fractionation and FACS analysisFor flow cytometry analysis, stromovascular cell (SVC) fractions from mouse iWAT were isolated as described previously3 and used for scRNA-seq analysis. For EdU detection, fixed SVCs were first processed for the Click-it reaction (ClickTech EdU Cell Proliferation Kit; Thermo Fisher) according to the manufacturer’s instructions, followed by cell surface marker staining. For cell surface marker staining, the following antibodies were used: anti-PDGFRA-APC (rat, 1:50; Biolegend, San Diego, CA, USA), anti-DPP4-FITC (rat, 1:50; Biolegend), anti-CD45-Pacific Blue (rat, 1:200; Biolegend), anti-NOTCH1-APC (hamster, 1:100; Biolegend), anti-PTCH1 (goat, 1:100; R&D Systems, Minneapolis, MN, USA), and donkey anti-goat IgG-Alexa Fluor 488 (1:200; Invitrogen, Waltham, MA, USA). Analytic cytometry was performed using a BD FACS LSRFortessa X-20 or FACS Lyric system (BD Biosciences, Franklin Lakes, NJ, USA). A FACSAria III system (BD Biosciences) was used to sort the cell surface marker-labeled cells. The raw data were processed using FlowJo software v.10 (TreeStar Inc., Ashland, OR, USA). The sorted cells were plated and cultured at an initial concentration of 1\(\times\)105 cells/ml in growth medium and then switched to adipogenic induction medium, as described previously3.HistologyAdipose tissue was processed to obtain histological sections, and 5 μm-thick paraffin sections were subjected to immunohistochemical analysis as previously described3. EdU analysis was performed according to the manufacturer’s instructions (Click-iT™ Plus EdU Alexa Fluor™ 488 Imaging Kit; Thermo Fisher). Species-matched IgGs were used or the primary antibodies were omitted as nonspecific controls for immunohistochemistry. All quantitative analyses of histologic samples was carried out in a blinded manner. The antibodies used for immunohistochemistry included tdTomato (rabbit, 1:100; mouse, 1:100; Clontech, Tokyo, Japan), UCP1 (0.5 mg/ml rabbit; Abcam, Cambridge, United Kingdom), PDGFRA (goat, 1:200; R&D Systems), PTCH1 (goat, 1:100; R&D Systems), and NOTCH1 (hamster, 1:50; Biolegend). For lipid staining, BODIPY or LipidTOX (Invitrogen) was used. For nuclear staining, DAPI (Sigma‒Aldrich) was used.Analysis of mRNA expression levelsRNA was extracted using TRIzol® reagent (Invitrogen), and 1 μg of RNA was reverse transcribed using a cDNA synthesis kit (High-capacity cDNA Reverse Transcription Kit; Applied Biosystems, Waltham, MA, USA). One hundred nanograms of cDNA was subjected to qRT‒PCR in 20-μl reaction volumes (iQ SYBR Green Supermix; Bio-Rad, Hercules, CA, USA) with 100 nM primers. qRT‒PCR was performed using SYBR Green dye and a CFX Connect Real-time system (Bio-Rad) for 45 cycles, and the fold change in expression for all the samples was calculated by using the 2−ΔΔCt method. Peptidylprolyl isomerase A (PPIA) was used as a housekeeping gene for mRNA expression analysis. The primers used for qRT‒PCR are listed in Supplementary Table 2.Oil red O stainingFor lipid droplet staining, Oil Red O staining was performed according to the manufacturer’s protocol. Briefly, cells were fixed with 4% paraformaldehyde for 30 min, washed with 60% isopropanol, and stained with Oil Red O (0.5% Oil Red O; Sigma‒Aldrich) in 60% isopropanol. The stain was extracted in isopropanol and measured at 492 nm using a microplate spectrophotometer (MultiSkan GO, Thermo Fisher).Dpp4 promoter reporter analysisHEK293T cells (ATCC) were transfected with DPP4 reporter vectors (MPRM34034; Genecopoeia, Rockville, MD, USA) with a JetPrime transfection kit. Additionally, pBABEPuro-HA-GLI1 (Addgene plasmid #62967), pLPCX-NICD1 (Addgene plasmid #44471), or pUC19 control plasmids (TaKaRa Bio, Tokyo, Japan) were cotransfected as indicated. At 24 h post-transfection, the cells were treated with either the Hedgehog inhibitor GANT61 (Cayman, Ann Arbor, MI, USA) or the gamma-secretase inhibitor DAPT (Cayman) to suppress Hedgehog and Notch signaling, respectively. After 24 h, the cell culture media were collected, and Gaussia luciferase (GLuc) and secreted alkaline phosphatase (SeAP) activities were measured via a Secrete-Pair Dual Luminescence Assay Kit (Genecopoeia) according to the manufacturer’s instructions. Ten microliters of media was used for each luminescence reaction. Luminescence signals were measured using a microplate luminometer (Centro LB 960, Berthold Technologies USA LLC, Oak Ridge, TN, USA). GLuc luciferase activity was normalized to SeAP activity.Hedgehog and notch signaling reporter analysisFor stable cell line establishment, HEK293T cells were transfected with pMuLE_EXPR_CMV-eGFP_TOP-NLuc1.1_12GLI-FLuc_CBF-GLuc (Addgene, plasmid #113862), pMD2.G (Addgene, plasmid #12259), and psPAX2 (Addgene, plasmid #12260) at a ratio of 4:1:3 using a JetPrime transfection kit. Lentiviral particles were collected over 2 days, and C3H10T1/2 cells were transduced with the particles and polybrene (Sigma‒Aldrich). Following transduction, the cells were selected with G-418 (AG Scientific, San Diego, CA, USA) for one week. Alternatively, transient transfection of pMuLE_EXPR_CMV-eGFP_TOP-NLuc1.1_12GLI-FLuc_CBF-GLuc into HEK293T cells was performed. For the GLuc assay, 10 μl of medium (without lysis) was added to 96-well white-bottom microplates, and 100 μl of coelenterazine buffer (Genecopoeia) was added. For the firefly luciferase (FLuc) assay, cells were lysed in Reporter Lysis Buffer (Promega) and centrifuged at 10,000 × g for 10 min at 4 °C, after which the supernatant was collected for analysis. Ten microliters of sample was loaded on 96-well white-bottom microplates, and 50 μl of luciferin buffer (Promega) was subsequently added. The luciferase activity of each construct was analyzed using a microplate luminometer (Centro LB 960). Luciferase activity values were normalized to the protein concentration of the lysate.scRNA-seq experimentsSVC fractions were isolated from mouse iWAT, and PDGFRA+ cells were subsequently enriched from the fractions using FACS, as described previously3. The sorted cells were stained with Trypan blue and resuspended at a concentration of 1 × 105 to 2 × 106 cells/ml. Cell viability was estimated to be greater than 80%. Using these cells, we performed single-cell library preparation according to the manufacturer’s protocol using the Chromium Single-cell 3′ Reagent Kit v2 (10× Genomics, Pleasanton, CA, USA). Briefly, cell suspensions were loaded onto each channel of a Chromium Single Cell A Chip along with reverse transcription master mix and single-cell 3′ gel beads to capture 3000 cells per channel. Gel beads-in-emulsions (GEMs) were generated by running the Chromium Controller, followed by reverse transcription using a Mastercycler Nexus system (Eppendorf AG, Hamburg, Germany). cDNA was amplified (12 PCR cycles) and purified using SPRIselect beads (Beckman Coulter, Brea, CA, USA). Purified cDNA was processed sequentially for enzymatic fragmentation, end repair, and A-tailing. Sequencing libraries with unique sample indices were pooled for multiplexing and subsequently sequenced using the Illumina HiSeq 2500 platform in high-output mode.scRNA-seq data processingThe Cell Ranger software suite (v2.0.0) provided by 10X Genomics was used to demultiplex sample barcodes, process cell barcodes (CBs) and UMIs, and construct a gene-by-cell count matrix. The raw BCL files were demultiplexed into FASTQ files using the cellranger mkfastq in Cell Ranger software. For each sample, paired-end reads in the FASTQ files were reformatted for processing CBs and UMIs, mapped to the mouse reference genome (GRCm38) by STAR (v2.5.2b) with an Ensembl GTF file (release 89), and quantified using CellRanger count in Cell Ranger software with the option of the expected number of recovered cells set to 3000 (–expect-cells = 3000). The output files of each sample were aggregated into a single gene-by-cell count matrix without normalization using cellranger aggr. We then filtered out CBs that corresponded to empty droplets after determining a threshold (called “a knee point”) in the barcode rank plot. All CBs with a total UMI less than the knee point were discarded. Next, we removed low-quality cells that had greater than 10% UMIs associated with mitochondrial-encoded genes and greater than 99.5% of the genes not expressed. The thresholds were chosen by visually inspecting outliers in the 2-dimensional principal component analysis (PCA) score plots of cell-level quality control (QC) metrics calculated by the scater R/Bioconductor package (v1.4.0), as suggested by Ilicic et al.12. Finally, for removal of cell-specific biases in the combined gene-by-cell count matrix, the raw counts were normalized by pooled size factors that were estimated from the scran R/Bioconductor package (v1.5.8). The normalized counts were then log2-transformed with a pseudocount of 1.scRNA-seq data analysisHighly variable genes (HVGs) were identified using the trendVar and decomposeVar functions of the scran package with a false discovery rate (FDR) ≤ 0.05 and biological variability >0.1. Because the combined dataset consists of two technical batches, the batch information was incorporated into a design matrix, which was subsequently used as an input into the trendVar function. To visualize each cell in a 2-dimensional space, we reduced the dimensionality using the t-SNE algorithm implemented in the RunTSNE function of the Seurat package (v2.1.0) with the first 20 principal components (PCs) obtained for the HVGs and perplexity = 30. To filter out immune cells based on hematopoietic markers, we grouped cells into 4 clusters using the k-means clustering algorithm. The cells were further grouped into 24 clusters using the FindClusters function of the Seurat package with the first 50 PCs for the HVGs. Of the 24 clusters, immune-related clusters identified via previous k-means clustering were filtered out for downstream analysis.After analysis of HVGs in nonimmune cells as described above, cells were visualized in a 2-dimensional t-SNE plot using the RunTSNE function of the Seurat package with the first 50 PCs and perplexity = 30. For clustering, we used the FindClusters function of the Seurat package with the first 50 PCs on the HVGs. The initial marker genes of each cluster were identified using the Wilcoxon rank sum test in the FindAllMarkers function of the Seurat package with an adjusted P value < 0.05 and a fold change >1.3. We then applied partial least square-discriminant analysis (PLS-DA) to the initial marker genes and computed the variable importance in projection (VIP), which represents the relative contribution of each gene to the separation of a cluster from the other clusters13. For each cluster, to determine a cutoff of significant VIPs, we estimated an empirical distribution of the null hypothesis for VIPs (i.e., a gene has no contribution to the separation) by applying the Gaussian kernel density estimation method14 to VIPs resulting from 1000 random permutations of cells. The adjusted P values of the VIPs of the individual genes were computed by the one-tailed test using the empirical null distribution. Finally, among the marker genes selected for each cluster by the rank sum test, we further selected the final marker genes as genes with an adjusted P value < 0.05 according to the VIP test.For trajectory inference analysis, we performed feature selection by identifying HVGs in cells belonging to the Progen, Prolif, and Diff groups and clustering these cells with the HVGs as described above. We selected the union of the marker genes across all clusters as a set of feature genes. We then normalized the raw counts by the scran size factor and reduced the dimension of the normalized matrix of the feature genes to 2 using the DDRTree method in the Monocle (v2.6.1)15 package. For each trajectory, we ordered cells along the trajectory using the orderCells function of the Monocle package. Finally, the union of the marker genes was clustered based on their smoothed expression profiles along the trajectory using the plot_pseudotime_heatmap function of the Monocle package, and the genes showing early, middle, or late upregulation along each trajectory were identified based on the clustering results. For each gene, the smoothed relative gene expression profile along each of the trajectories was plotted using the RollApply function of the zoo R package (v1.8-3) with a window size of 150 (control) or 300 (CL-treated).Enrichment analysis of GO biological processes and pathwaysFunctional enrichment analysis of the marker genes was performed using DAVID software16. The enrichment of Gene Ontology biological process (GOBP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways associated with the marker genes was considered to be associated with P < 0.05.Network analysisThe network model for the selected marker genes was visualized using Cytoscape17. The nodes in the network model were arranged based on their locations in their associated pathways in the KEGG pathway database, and the edges in the network model were defined to display their activation/inhibition information obtained from the KEGG pathways.

Hot Topics

Related Articles