Primary ASCs were isolated from subcutaneous fat obtained by liposuction from six unrelated female donors ( Supplementary Table 1 ) after informed consent (approval by the Regional Committee for Research Ethics for Southern Norway, REK 2013/2102 and 2018-660). ASCs were cultured in DMEM/F12 (17.5 mM glucose) with 10% fetal calf serum and 20 ng/ml basic fibroblast growth factor (proliferation medium). Upon confluency, fibroblast growth factor was removed, and cells cultured for 72 h in DMEM/F12 (17.5 mM glucose) with 10% fetal calf serum (basal medium) before induction of differentiation. For white adipose differentiation, ASCs were induced with a cocktail of 0.5 µM 1-methyl-3 isobutyl xanthine, 1 µM dexamethasone, 10 µg/ml insulin and 200 µM indomethacin in basal medium. Differentiation media was renewed every 3 days until day 9, after which cells were maintained in DMEM/F12 (17.5 mM glucose) with 10% fetal calf serum and 10 µg/ml insulin. For beige adipose differentiation, media were supplemented with 1 µM Rosiglitazone (Sigma, R2408) until day 15. Samples were harvested 15 days after induction. All differentiation experiments were done in at least three biological replicates between passage 3 and 8.

Cells were differentiated for 15 days on 12-mm diameter coverslips in 24-well plates. Cells were washed 3 times with PBS before fixation in 4% paraformaldehyde for 10 min. Cells were then incubated in Bodipy (1µg/ml, Invitrogen D3922) for 15 min and washed 3 times in PBS before mounting in DAKO Fluorescence Mounting Medium (S3023, Agilent). Images were acquired on an IX81 microscope (Olympus) fitted with epifluorescence, an 100× 1.4 NA objective mounted on a piezo drive, and a DeltaVision personalDV (Applied Precision, Ltd.) imaging station. The Bodipy surface per field was quantified on threshold-adjusted images using ImageJ software.

Proteins were resolved by gradient 4-20% SDS–PAGE, transferred onto nitrocellulose membranes (BioRad) and blocked with 5% BSA. Membranes were incubated for 1h at room temperature or overnight at 4°C using the following antibodies: Perilipin1 (Progen, GP29), CITED1 (Novus, H00004435-M03), UCP1 (Abcam, 23841), PPARγ (Thermofisher, MA5-14889), PEMT (Novus, NBP1-59580), and γTubulin (Sigma, T5326). Proteins were visualized using IRDye-800-, IRDye-680-, or HRP-coupled secondary antibodies. Bands were quantified by densitometry (Image Lab; BioRad) using γTubulin for normalization. Uncropped membranes are presented as Supplementary Data 1 - 4 .

RNA sequencing (RNA-seq) was done in biological triplicates for 6 independent donors. Total RNA was isolated from samples harvested at differentiation endpoint (day 15) using the RNeasy kit (QIAGEN). PolyA-selected RNA was sequenced from paired-end libraries (TruSeq Stranded mRNA kit; Illumina) using Novaseq platform (Illumina). Reads were filtered with fastp, aligned to the hg38 genome (GENCODE 32 annotation) with hisat2, and counted using featureCounts ( Supplementary Table 2 ). Low abundance genes were filtered using filterByExpr from edgeR and then normalized using the trimmed mean of M values method (29). Beige and white gene expression was compared pairwise for each donor using the robust eBayes method with limma-voom adjustment (30). Less than 100 genes changed according to the fraction of reads assigned to genes (adjusted p-value < 0.01), when this was included as a continuous variable in the model.

RNA was isolated using RNeasy kit (QIAGEN) and 1 µg was used for cDNA synthesis using High-Capacity cDNA Reverse Transcription Kit (ThermoFisher). RT-PCR was done using IQ SYBR green (Bio-Rad Laboratories) with SF3A1 as a reference gene. PCR conditions were 95°C for 3 min and 40 cycles of 95°C for 30 s, 60°C for 30 s, and 72°C for 20 s. PCR primers are listed in Supplementary Table 3 .

Reads were aligned to the hg38 genome with STAR using the GENCODE 32 reference annotation. ENCODE options from the STAR manual were used as well as –outSAMstrandField intronMotif –outSAMtype BAM Unsorted ( Supplementary Table 2 ). Split reads with at least 6 bp of anchor sequence and an intron length between 20 bp and 1 Mb were extracted using regtools (31). Introns sharing splice sites were clustered using LeafCutter’s leafcutter_cluster_regtools.py script, requiring at least 30 reads per library supporting each cluster for introns of up to 1 Mb in size.

Differences in intron excision were tested between white and beige triplicates with donor as a confounding variable. Strict prefiltering required that an intron was found in at least 14 (of 36) libraries and that clusters had at least 12 libraries per condition with 20 or more spliced reads. The differential_splicing function from LeafCutter R package version 0.2.9 was used to implement the Dirichlet-multinomial generalized linear model (15). LeafCutter assigns a ΔPSI to each intron and a p-value to each cluster, which contains introns that share splice sites and thus belong to alternative isoforms of that gene or genes. Effect size and cluster significance output files were merged into a single table containing both ΔPSI and cluster p-value ( Supplementary Table 4 ). LeafCutter p-values are adjusted by FDR.

Introns were assigned to transcripts, with the priority for annotations: GENCODE 32 > RefSeq GCF_000001405.40-RS_2023_10 > FANTOM CAT robust (32). In other words, only if an intron was missing from the GENCODE annotation was it compared to the RefSeq annotation and so on. Normalized TRIFID prediction scores were downloaded for the GENCODE 37 and RefSeq110 annotations. Where multiple transcripts were possible for a single intron excision event, TRIFID scores were averaged, but if no score was found -0.1 was used. To calculate a TRIFID difference, the top two intron-excision events were compared for clusters with at least one significant junction (| ΔPSI | > 0.1 & p.adjust < 0.05) ( Supplementary Table 5 ).

For ChIP-Seq profiles at TSSs, for each significant exon-exon junction, the most upstream annotated TSS was selected to avoid oversampling of junctions. For Unibind differential enrichment (beige vs white) at TSSs, overlapping -2 kb/+0.5 kb windows were merged.

Day 15 differentiated white and beige samples were trypsinized, resuspended in HBSS, 0.5% BSA, and centrifuged 200g for 5 min to isolate floating mature adipocytes. Purified nuclei were fixed with 1% formaldehyde; lysed for 10 min in ChIP lysis buffer (1% SDS, 10 mM EDTA, 50 mM Tris-HCl, pH 8.0, proteinase inhibitors, 1 mM PMSF, 20 mM Na Butyrate) and sonicated for 30 sec ON/OFF for 10 min in a Bioruptor^®Pico (Diagenode) into ~200 base-pair fragments. After sedimentation, the supernatant was diluted 10 times and chromatin incubated with anti-H3K27ac (Diagenode c15410174), anti-H3K4me3 (Diagenode c15410003), anti-H3K27me3 (Diagenode c15410069) or anti-H3K4me1 (Diagenode c15410037) antibodies, each at 2.5 μg/106 cells, for 2 h at 4°C. ChIP samples were washed, cross-links reversed, and DNA eluted for 2 h at 68°C. DNA was purified using phenol-chloroform isoamylalcohol and dissolved in H₂O. Libraries were prepared using a Microplex kit (Diagenode) and sequenced on a NextSeq 500 or NovaSeq (Illumina).

ChIP-Seq reads for histone modifications were filtered with fastp, aligned to the genome with bowtie2, and deduplicated ( Supplementary Table 2 ). Log ratio ChIP/Input tracks were generated using deepTools bigwigCompare from reads per genome coverage (RPGC) normalized tracks, created with bamCoverage. Profile plots were generated using deeptools computeMatrix and plotProfile.

PPARγ and MED1 ChIP-Seq were downloaded from Gene Expression Omnibus (GEO) accession GSE59703 (33). Reads from replicates were merged and aligned to hg38 using Bowtie2 ( Supplementary Table 2 ). Peaks were detected using MACS3. ChIP/Input ratio tracks were generated using bamCompare in deepTools. PPARγ and MED1 ChIP/Input mean ratios were calculated for +/- 250 bp windows around the TSS using multiBigwigSummary from deepTools.

Energy minimized structures of PEMT isoforms were predicted using Alphafold via ColabFold (34). Molecular graphics and analyses performed with UCSF ChimeraX (35).

Fisher’s exact test, implemented in ClusterProfiler (36), was used for overrepresentation analysis of genes in ontologies and pathways from MSigDb v2023.1 (37). For bioinformatic analyses, two-way ANOVA were conducted in R, followed by pairwise t-tests or Wilcoxon tests with Holmberg adjustment for multiple comparisons.

To characterize the alternative splicing events that may specify the beige thermogenic phenotype of human adipocytes, we differentiated primary ASCs from 6 unrelated, normal-weight female subjects (S1-S6) into beige or white adipocytes, in the presence or absence of rosiglitazone, respectively. After 15 days, differentiation efficiency is similar in all 6 ASC line-derived adipocytes, based on lipid droplet accumulation ( Figures 1A, B ) and Perilipin1 expression ( Figures 1C, D , Supplementary Figure 1 ). As expected, rosiglitazone treatment elicits the expression of the beige adipocyte markers CITED1 and UCP1 in all donors ( Figure 1C ), albeit to a different extent, confirming that efficient beige adipogenesis is achieved with all 6 ASC lines. Accordingly, transcriptome analysis shows that differentially expressed genes (DEGs) in beige vs white adipocytes from all subjects are significantly enriched for Gene Ontology (GO) terms “Mitochondrial matrix” (GO:0005759), “Monocarboxylic acid metabolic process” (GO:0032787) and “Fatty acid metabolic process” (GO:0006631) ( Figure 1E , Supplementary Table 4 ). Further analysis of the thermogenic gene expression signature using BATLAS beige-curated gene list (38) clearly segregates white from beige differentiated adipocytes ( Figure 1F ), with the median expression of BATLAS beige markers being significantly increased in beige adipocytes ( Figure 1G ). Hence, PPARγ agonist treatment efficiently elicits a beige adipocyte transcriptional program in all 6 ASC lines.

We next assessed differential splicing events in the white and beige adipocyte transcriptomes using LeafCutter, an annotation-free method that allows the identification and quantification of alternative splicing events by focusing on intron excision events, which are inferred from reads spanning exon–exon junctions (15). LeafCutter analysis identifies 777 exon junctions with robust differential usage between white and beige adipocytes in all 6 subjects ( Figure 2A , Supplementary Table 5 ). Strikingly, genes with highly significant differences in intron excision, quantified by percent spliced-in (ΔPSI), include beiging markers (CITED1, PANK1), genes involved fatty acid and phospholipid remodeling (FAR2, PEMT, PLD1, LPIN1) or glucose metabolism (PC) as well as PPARG itself. Amongst the cell type specific exon junctions detected by LeafCutter, 723 could be annotated to known transcripts from either Gencode, RefSeq, or Fantom databases ( Figure 2B ). However, the remaining 54 cryptic exon junctions display canonical 5’ and 3’ motifs, and mostly result from differential 5´ splice site usage ( Supplementary Figures 2A-C ). Differential exon junctions detected in the beige vs white adipocyte transcriptomes map to 562 genes. Strikingly, only a minor fraction of these differentially spliced genes (DSGs) are also differentially expressed ( Figure 2C ), indicating that differential isoform expression constitutes an additional mechanism for gene expression regulation during beige adipocyte differentiation.

To assign functionality to differential exon junction usage, we mapped each junction to its matching transcripts and leveraged the TRIFID database to score protein functionality (18). TRIFID predicts protein functionality based on a number of transcript features including the presence of functional domains, cross-species conservation, transcript length and annotation database agreement. The TRIFID score (between 0 and 1) represents the likely functional relevance of protein isoforms. Both beige- and white-specific transcripts using spliced-in exon junctions have significantly higher average TRIFID functionality scores compared to the total transcriptome and to non-enriched exon-exon junctions from the same genes ( Figure 2D ), supporting a functional outcome for alternative isoform expression. Indeed, differential TRIFID score (ΔTRIFID) for DSGs predicts higher functionality for PC, CA5B, CITED1, RTN4 and PPARG isoforms in beige adipocytes ( Figure 2E , Supplementary Table 6 ), and genes related to “Adaptive thermogenesis”, “Fatty acid metabolism” and “Mitochondrial matrix” show either an increase or decrease in predicted isoform functionality in beige adipocytes ( Figure 2F ). Altogether, these results argue for a functionally relevant impact of alternative isoform expression on the regulation of beige vs white adipocyte functions.

To understand the mechanism driving beige-specific intron excision events, we first mapped these events along transcripts. We find a large majority of differential events occurring in the first intron of transcripts ( Figure 3A ), suggesting that DSGs between white and beige adipocytes may result from alternative TSS usage. This is further supported by the absence of consistent changes in expression of splicing regulators between white and beige adipocytes, and the minor overlap with splicing-driven alternative transcript expression described for adrenergic-induced adipocyte beiging ( Supplementary Figures 3 , 4 ) (28). We thus examined the epigenetic state of white and beige specific TSSs in mature adipocytes by ChIP-seq. Enrichment profiles for the active histone marks H3K4me3 and H3K27ac at beige TSSs are similar in beige and white adipocytes ( Figures 3B, C , Supplementary Figures 5A, B ), indicating that beige TSSs are already in an active state in white adipocytes. In line, we find a similar bimodal H3K4me1 enrichment around white and beige TSSs, and no enrichment for the repressive histone mark H3K27me3 is detected at beige TSSs in white adipocytes ( Supplementary Figures 5C-G ). Hence, differential TSS usage in beige vs white adipocytes is not driven by a transition in chromatin state.

We reasoned that beige TSS activation might result from differential binding of transcription factors. Using Unibind’s binding site predictions (39), we find PPARγ and CEBPα transcription factor binding sites as the most differentially enriched in beige compared to white TSSs ( Figure 3D ). Analysis of PPARγ ChIP-Seq data in human white and beige adipocytes (33) confirms that a higher proportion of beige TSSs overlap PPARγ peaks compared to white TSSs and non-significantly enriched TSSs ( Figure 3E ). Accordingly, average PPARγ enrichment levels are higher at beige TSSs in both beige and white adipocytes ( Figures 3F, G ). Since cell type-specific enhancers are enriched for MED1 binding in beige vs white adipocytes (33), we examined MED1 binding in each TSS subclass ( Figure 3H ). Similar to PPARγ, MED1 signals are significantly enriched at beige vs white TSSs ( Figure 3I ). However, rosiglitazone treatment does not result in increased MED1 binding, indicating that promoter-enhancer contacts at beige TSSs are also established in white adipocytes ( Figures 3H, I ). Altogether, our results show that beige-specific TSSs are in an active chromatin state and bound by PPARγ in both white and beige adipocytes, implying that rosiglitazone activation of PPARγ triggers cell type-specific isoform expression in beige adipocytes.

LeafCutter analysis reveals that white and beige adipocytes express distinct proportions of PPARG isoforms (see Figure 2A ). Indeed, exon junction reads spanning PPARG1 exon 1 and 2 are overrepresented in white adipocytes, while PPARG2 first exons are most expressed in beige adipocytes ( Figures 4A-C ). In agreement with whole TSS analyses, PPARG1 and PPARG2 TSSs are similarly enriched with broad domains of H3K4me3 and H3K27ac, marking active promoters ( Figure 4B ). However, only the PPARG2 TSS region shows an enrichment for PPARγ transcription factor binding and several MED1 peaks, suggestive of promoter-enhancer interaction ( Figure 4B ). Interestingly, PPARG2 TSS activation does not result in a global increase in PPARG expression in beige adipocytes ( Figure 4D ), but rather in a switch in the proportion of PPARG isoforms as confirmed by isoform-specific RT-qPCR ( Figure 4E ). Importantly, we show that both isoforms of PPARγ are readily detected at the protein level ( Figure 4F ), and rosiglitazone treatment results in a shift in PPARγ isoform expression with PPARγ2 becoming the main isoform in all ASC lines ( Figure 4G ), while total PPARγ protein expression remains unchanged ( Figure 4H ). Thus, altered balance of PPARγ isoform expression in beige adipocytes may underlie the activation of beige-specific TSSs upon rosiglitazone treatment.

Amongst highly significant DSGs is PEMT, encoding the phosphatidylethanolamine methyltransferase ( Figures 5A–C ; see Figure 2A ), an enzyme found at mitochondria-associated membranes and involved in the regulation of thermogenesis in vivo (40). PEMT encodes two major protein isoforms: short PEMT-S is liver-specific and highly active, while the long variant PEMT-L has an additional 37 N-terminal amino acids, reduced activity, and is expressed at low levels across a broad range of tissues (41) ( Figure 5D ; Supplementary Figures 6 and 7 ). Unlike most DSGs, PEMT expression is strongly induced in beige adipocytes ( Figure 5E ; see Figure 2C ) and RNA-seq read alignment suggests this is driven by upregulation of a novel isoform ( Figures 5A, B ).

Indeed, LeafCutter identifies a switch from PEMT-L to a cryptic isoform of PEMT (hereafter referred to as PEMT-C) in beige adipocytes, with a first exon-exon junction matching a predicted transcript from RefSeq (XM_006721418.5) ( Figures 5A-C ), with a TRIFID functionality score similar to that of the PEMT-L isoform (see Figure 2F ). However, no RNA-seq reads map over the first 250 bp of XM_006721418.5’s first intron, suggesting PEMT-C rather arises from a second downstream TSS, which is readily detected in Fantom robust CAGE-seq clusters ( Figure 5B , right panel), and is a prominent PEMT TSS in differentiated adipocytes based on Fantom5 CAGE-Seq signal ( Supplementary Figure 8 ). While PEMT-C promoter is similarly bound by PPARγ and MED1 in white and beige adipocytes, its upregulation in beige adipocytes correlates with an enrichment for the active promoter marks H3K4me3 and H3K27ac around its TSS ( Figure 5B ); this implies that the recruitment of chromatin modifiers underlies PEMT-C isoform expression in beige adipocytes.

Analysis of PEMT isoform expression using primers spanning isoform-specific exon junctions reveals that PEMT upregulation in beige adipocytes results from the induction of PEMT-C isoform only, while PEMT-L expression remains constant ( Figure 5F ). Importantly, only one PEMT isoform is detected at the protein level, with a molecular weight similar to that of the canonical PEMT-S isoform ( Figure 5F ). In agreement with RNA-seq and qPCR analysis, we find PEMT protein expression is significantly increased in beige adipocytes derived from all 6 ASC lines ( Figures 5G, H ). Thus, beige adipocytes express high levels of a novel PEMT isoform variant which is structurally close to the highly active PEMT-S variant ( Figure 5D ). Altogether, our results indicate that differential alternative TSS usage significantly impacts the beige adipocyte proteome.