A non-coding variant linked to metabolic obesity with normal weight affects actin remodelling in subcutaneous adipocytes

Abstract

Recent large-scale genomic association studies found evidence for a genetic link between increased risk of type 2 diabetes and decreased risk for adiposity-related traits, reminiscent of metabolically obese normal weight (MONW) association signatures. However, the target genes and cellular mechanisms driving such MONW associations remain to be identified. Here, we systematically identify the cellular programmes of one of the top-scoring MONW risk loci, the 2q24.3 risk locus, in subcutaneous adipocytes. We identify a causal genetic variant, rs6712203, an intronic single-nucleotide polymorphism in the COBLL1 gene, which changes the conserved transcription factor motif of POU domain, class 2, transcription factor 2, and leads to differential COBLL1 gene expression by altering the enhancer activity at the locus in subcutaneous adipocytes. We then establish the cellular programme under the genetic control of the 2q24.3 MONW risk locus and the effector gene COBLL1, which is characterized by impaired actin cytoskeleton remodelling in differentiating subcutaneous adipocytes and subsequent failure of these cells to accumulate lipids and develop into metabolically active and insulin-sensitive adipocytes. Finally, we show that perturbations of the effector gene Cobll1 in a mouse model result in organismal phenotypes matching the MONW association signature, including decreased subcutaneous body fat mass and body weight along with impaired glucose tolerance. Taken together, our results provide a mechanistic link between the genetic risk for insulin resistance and low adiposity, providing a potential therapeutic hypothesis and a framework for future identification of causal relationships between genome associations and cellular programmes in other disorders.

Extended data Fig. 1 ∣. The pleiotropic 2q24.3 MONW locus is associated with increased risk for type 2 diabetes and decreased adiposity related traits and maps to sparse enhancer signatures in adipocytes.

(a) Schematic overview for the 2q24.3 metabolic risk locus dissection. Aim of step (top, bold); methods/experiments used (middle); key finding/result of each step (bottom). (b) PheWAS of trait associations at the rs3923113-tagged haplotype of a meta-analysis https://t2d.hugeamp.org/. Colors represent trait classes while individual rs3923113 variant association p-values are shown on the Y axis. Direction of effect is indicated by orientation of triangles, upward: increase, downward: decrease (c) The 2q24.3 MONW locus spans 19 non-coding SNPs in high linkage disequilibrium with rs3923113 (LD r2 > 0.8). The region of association localizes to a >55 kb interval in an intergenic region between COBLL1 and GRB14. (d) Annotation panel and color key for the twenty-five state chromatin model70. Rows represent chromatin states abbreviations, columns are emission parameters, corresponding to the frequency with which each mark is expected in each state (left table) and genome coverage and median enrichments of relevant genomic annotations (right panel). TssA: Active TSS, TssAFlnk: Flanking Active TSS, TxFlnk: Transcription at gene 5’ and 3’, Tx: Strong Transcription, TxWk: Weak Transcription, EnhG: Genic enhancers: Enh: Enhancers, ZNF/Rpts: ZNF genes & repeats, Het: Heterochromatin, TssBiv: Bivalent/Poised TSS, BivFlnk: Flanking Bivalent TSS/Enhancer, EnhBiv: Bivalent Enhancer, ReprPC: Repressed Polycomb, ReprPCWk: Weak Repressed Polycomb, Quies: Quiescent/Low. (e) Stranded allele-specific chromatin accessibility measures at the haplotype using ATAC-seq data in differentiating adipocytes from a heterozygous individual. For each day of differentiation of an individual heterozygous, the number of reads overlapping with 20 non-coding SNPs in the haplotype, ordered by their start position and strand relative to the position of the variant, are shown. More reads indicate higher activity in haplotype 1 (non-risk, blue) compared to haplotype 2 (risk yellow). x-axis: offset from SNP position (bp), y-axis: stranded read count. (f) Replication of the effect at time 0 (mesenchymal stem cells) with ATAC-seq. (g) BMI-dependent variant association analysis. Bar plots represent the beta of the rs6712203 association with type 2 diabetes following BMI stratification. The cohort analysed is the UK Biobank self-identified white British individuals (total N = 327,960; N = 109198 with BMI < 25, N = 140539 with BMI between 25 and 30, and N = 78223 with BMI > = 30), and overlay of data points is not practical. Betas and 95% confidence intervals are shown, derived from a two-sided generalized linear model on outcome adjusted for demographic covariates (age, sex, genotyping array, 40 PCs).

Extended data Fig. 2 ∣. Conditional ana- lyses implicating rs6712203 in the genetic control of anthropometric traits and type 2 diabetes.

a, Conditional analyses implicating rs6712203 in the genetic control of anthropometric traits and type 2 diabetes. Each panel represents a different trait/sex/conditional analysis window, and all panels have an X axis corresponding to 100 kb on either side of the rs6712203 variant. The Y axis shows, for each variant in the window, the association strength for the given trait conditioned on the variants noted in White British participants in UK Biobank with the sex shown, and red lines indicate the significance threshold 5 × 10-8). −log10 p-values are shown, derived from a two-sided generalized linear model on outcome adjusted for demographic covariates (age, sex, genotyping array, 40 PCs).

Extended data Fig. 3 ∣. Chromatin inter- actions and CRISPRi of 2q24.3 locus identify COBLL1 as target genes.

(a) Cross-cell type conserved genome-wide higher order chromatin interactions for the 2q24.3 locus analyzed by Hi-C assays in human fibroblasts (left) and NHEK primary normal human epidermal keratinocytes (right), chr2: 163,556,000 - 167,558,000 (hg19), binned at 2 kb resolution. (b) Cas9 protein expression in dCas9 hWAT compared to the parental hWAT cell line. (c) mRNA expression of COBLL1 and GRB14 in response to increasing amounts of lentiviral sgRNA vectors (2 sgRNAs, virus volume 50 μl and 500 μl) targeting TSS regions of each gene compared to non-targeting controls (NT, 2 sgRNAs). Columns are means of individual sgRNAs indicated by different symbols. (d) COBLL1 protein expression normalized to b-actin in dCas9 hWATs transduced with sgRNAs targeting COBLL1 or GRB14 compared to controls. Top panel: Image of gel of representative sgRNA targeting NT, COBLL1 or GRB14. Bottom panel: plot of protein expression; 2 sgRNA for each target in 2 replicates. (e) Representation of 1,181 bp region flanking the COBLL1 intronic variant rs6712203 at the 2q24.3 MONW locus showing individual sgRNAs (n = 6) targeting the rs6712203 flanking regulatory region used in the CRISPRi experiments. (f, g) mRNA expression of (f) COBLL1 and (g) GRB14 in undifferentiated dCas9-hWAT preadipocytes at 6 days post lentiviral transduction with sgRNAs targeting TSS regions (red: COBLL1 TSS; blue GRB14 TSS) and the rs6712203-flanking regulatory element at position 1 to 6 as depicted in (e). Data are mean +/− SEM of 3 independent experiments. **** P < 0.0001, *** P = 0.0004, ** P = 0.006, * P = 0.013 – 0.036, two-tailed Student’s t test. (h) Predicted binding of POU2F2 between the two alleles using the Intragenomic Replicate Method (Cowper-Sal lari et al. 2012). As in Fig. 2d with different kmer counts.

Extended data Fig. 4 ∣. COBLL1 regulates actin cytoskeleton remodeling.

(a) COBLL1 expression in subcutaneous and visceral AMSCs throughout adipogenic differentiation, N = 4 biologically independent experiments, t-test two-sided, data represent median + 95% CI. (b) COBLL1 gene expression enrichment across 142 tissues (A-D) from enrichment profiler36. COBLL1 probes 203641_s_at and 203642_s_at were used for coregulation analysis (E-F). (c) Correlation with COBLL1 probe ILMN_1761260 using microarray data from lean and individuals with obesity. (d) Enrichment of pathways in the HCI (upper panel) and WikiPathways (lower panel) gene set lists from Enrichr, plotted as in Fig. 3A (KEGG), with p-value thresholds corresponding to the FDR cutoffs in those data. p-values are derived from a hypergeometric test. (e) COBLL1 expression in subcutaneous adipose tissue before and after a very low caloric diet (VLCD, upper panel, n = 18), corresponding body weight (lower panel), Wilcoxon signed-rank test.

Extended data Fig. 5 ∣. Knockdown of COBLL1 affects actin remodeling processes in differentiating adipocytes along with adipocyte differentiation insulin sensitivity and lipolysis rate.

(a) COBLL1 expression in siCOBLL1 and siNT at day 0, 3 and 14 of differentiation, N = 3 biologically independent experiments, t-test two-sided. knock-down efficiency 80%, mean values + SEM. (b–d) Morphological profiles of siCOBLL compared to siNT AMSCs at day 0 (b) day 3 (c) and day 9 (d) of differentiation, t-test two-sided, significance level < 5%FDR. (e) Actin and COBLL1 staining in siCOBLL1 compared to siNT subcutaneous adipocytes at day 9 using phalloidin and COBLL1 antibody staining (HPA053344, Alexa-Fluor 488), magnification x63/oil. Scale bar = 52.8 um. Representative results from N = 3 independent experiments. (f) Cells_Children_ LargeBODIPY_objects_count in siCOBLL1- and siNT AMSCs at day 3, 9, 14, N = 3 biologically independent experiments, t-test two-sided, significance level < 5% FDR. (G) qPCR-based gene expression of COBLL1 and adipocyte marker genes GLUT4, FASN, LIPE, PPARG, PLIN1, FABP4, CEBPA, ADIPOQ in siCOBLL1 and siNT AMSCs at day 14 of differentiation, t-test two-sided, N = 4 biologically independent experiments, mean values +/− SEM. (h) qPCR-based leptin gene expression in shCOBLL1 compared to shEV adipocytes. Data are represented as median + 95% CI, one-way ANOVA with Tukey’s HSD test, N = 4 biologically independent experiments (i) Correlation of COBLL1 mRNA with LEP mRNA in subcutaneous adipose tissue from 24 lean individuals measured by Illumina microarrays. The pearson’s correlation coefficient r and p-value are depicted (j) Schematic of siCOBLL1 KD and AMSCs differentiation. (k) UMAP-based dimensionality reduction of LipocyteProfiler features in siCOBLL1 and siNT AMSCs. (l) Actin and COBLL1 staining in siCOBLL1 and siNT visceral adipocytes at day 14 using phalloidin and COBLL1 antibody staining (HPA053344, Alexa-Fluor 488), magnification x63/oil. Representative result from N = 2 independent experiments, scale bar = 52,8um (m) Representative Oil-Red-O lipid staining in SGBS adipocytes following lentiviral COBLL1 knock-down (shCOBLL1, knock-down efficiency 69%) and GRB14 (shGRB14, knock-down efficiency 61%) compared to empty vector control (shEV), scale bar = 15 mm. (n) GPDH metabolic activity in shCOBLL1, shGRB14 and shEV SGBS adipocytes, one-way ANOVA with Tukey’s HSD test, mean + 95% CI, N = 4 biologically independent experiments (o) Basal and insulin-stimulated 3H-2-deoxyglucose uptake in shCOBLL1, shGRB14 and shEV SGBS adipocytes, one-way ANOVA with Tukey’s HSD test, mean + 95% CI, N = 4 biologically independent experiments, 1st and 3rd quartiles (box) and median (middle line) are indicated, p = 4.3 × 10-8. (p) qPCR-based GLUT4 gene expression in shCOBLL1, shGRB14 and shEV adipocytes, one-way ANOVA with Tukey’s HSD test, mean + 95% CI, N = 4 biologically independent experiments.

Extended data Fig. 6 ∣. Lipocyte profiles of risk versus non-risk haplotype carriers.

(a–c) Differences in morphological profiles between TT (n = 7) and CC (n = 6) allele carriers at day 0 (a), day 3 (b) and day 8 (c) in subcutaneous AMSCs (multi-way ANOVA, significance level < 5% FDR). (d–f) Differences in morphological profiles between TT (n = 7) and CC (n = 6) allele carriers at (d) day 0, (e) day 3 and (f) day 8 in visceral AMSCs (multi-way ANOVA, significance level < 5% FDR).

Extended data Fig. 7 ∣. Generation of COBLL1 mutant mice using CRISPR/Cas9 editing.

(a) Overview of the CRISPR/− Cas9 strategy to delete ~20 kb of the Cobll1 gene. The gRNA-targeting sequences (gRNAs) are underlined, and the PAM sequences are indicated in bold. Exons are represented as thick black boxes introns are indicated as black lines with arrows, and the yellow boxes indicate the DNA-targeting region. Red hexagon indicates a stop codon gen- erating a Cobll1 truncated protein. Agarose gel showing the PCR products generated from DNA containing success- fully targeted Cobll1 from F0 mouse tail genomic DNA. The 308 bp band corresponds to the genomic deletion. (b) A real-time quantitative PCR of levels of Cobll1 mRNA in white adipose tissue (WAT), liver and kidney of Cobll1 WT, Cobll1 heterozygous (+/−) and null knockout Cobll1 (−/−) animals to confirm the Cobll1 ablation in knockout animals. Each group was analyzed using 5 different mice and the values were expressed as the mean ± s.e.m and P values by Student’s t-test the experiment was repeated independently two times with similar results. (c) Pie chart illustrating non-redundant differential features per channel and class of measurement at day 8 of subcutaneous adipocyte differentiation in rs6712003 homozygous risk compared to non-risk carriers. (d, e) Differences in morphological profiles between AMSCs from Cobll1−/− mice (n = 3) and WT (n = 4) at (h) day 0 (i) day 2 (t-test two-sided, significance level < 5%FDR).

Fig. 1 ∣. The pleiotropic 2q24.3 MONW locus is associated with increased risk for T2D and decreased adiposity-related traits, and maps to sparse enhancer signatures in adipocytes.

a, Phenome-wide association studies of trait associations at the rs3923113-tagged haplotype in the UKB. The colours represent trait classes while individual rs3923113 variant association P values are shown on the y axis. The direction of effect is indicated by the orientation of the triangles: upward, increase; downward, decrease. b, Chromatin state annotations for the 55-kb-long MONW risk locus. Genomic intervals are shown across 127 human cell types and tissues reference epigenomes profiled by the Roadmap Epigenomics project, based on a 25-state chromatin state model (see Extended Data Fig. 1d for the color code of the chromatin states) learned from 12 epigenomic marks using imputed signal tracks at 25-nucleotide resolution. Chromatin states considered included Polycomb repressed states (grey, H3K27me3), weak enhancers (yellow, H3K4me1 only), strong enhancers (orange, also H3K27ac) and transcribed enhancers (green, also H3K36me3) (https://www.freepatentsonline.com/y2022/0243178.html). ES, embryonic stem; ESC, embryonic stem cell; GI, gastrointestinal; iPSC, induced pluripotent stem cell; MSC, mesenchymal stem cell.

Fig. 2 ∣. rs6712203 is a functional variant at the 2q24.3 MONW locus.

a, Phylogenetic conservation analysis and CNN-based prediction of chromatin accessibility for 20 highly linked (LD = r² > 0.8) variants at the 2q24.3 locus (tag SNP rs3923113). Phylogenetic conservation scores of jointly conserved motifs using PMCA are shown on the x axis. PMCA was used to identify orthologous regions in 20 vertebrates and to scan the 120-bp sequence context for groups of transcription factor binding site motifs whose sequence, order and distance range is conserved across species. Scores indicate the count of non-overlappingjointly conserved transcription factor binding site motifs whose relative positions within the window are conserved. Predicted relative change in chromatin accessibility (SNP accessibility difference SAD scores) in preadipocytes (day 0 of differentiation) for each SNP comparing alleles on haplotype 1 and haplotype 2 is shown on the y axis. A deep CNN Basset was trained on genome-wide ATAC–seq data assayed in preadipocytes. Alleles were assigned to each SNP in the haplotype and evaluated for predicted accessibility. b, Both PMCA and Basset predicted rs6712203 as a functional variant. For the C allele, there are no nucleotide variants that reduce binding in the rs6712203 region. Each position on the x axis represents a nucleotide and the four values in the heatmap correspond to all possible substitutions. c, For the T allele, in silico saturation mutagenesis suggests that binding loss of the POU2F motif that overlaps rs6712203, including the C allele itself, results in reduced predicted chromatin accessibility. d, Intragenomic replicates predict a higher binding affinity of POU2 family transcription factors for the T than C allele to both strands. Offsets from instances of the given k-mer sequence are shown on the x axis. Estimated affinity of binding (https://www.freepatentsonline.com/y2022/0243178.html) is shown on the y axis. A model with 8-mers is shown; alternatives with 6-mers through to 9-mers are shown in Extended Data Fig. 3b. e, Generation of rs6712203 CRISPR–Cas9 engineered lines starting from SGBS preadipocytes (heterozygous for rs6712203) edited to the homozygous risk (CC, yellow) and non-risk (TT, blue) alleles and qPCR-based gene expression measurement of COBLL1 conditional on the regulator POU2F2. n = 3 biologically independent experiments. f, Schematic model of the regulatory circuitry under the genetic control of rs6712203.

Fig. 3 ∣. The 2q24.3 effector gene COBLL1 affects actin remodelling processes in differentiating adipocytes.

a, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of genes correlated with COBLL1 in adipocytes using Enrichr^,. n = 30 donors, P values were derived from a hypergeometric test. b, Schematic of siCOBLL1 experiments in AMSCs across differentiation. AMSCs from a normal-weight female donor were silenced 3 d before induction and LipocyteProfiling (days 0, 3, 8 and 14 of differentiation), n = 3. c, Image-based profiles of siCOBLL1-treated compared to non-targeting siRNA-treated AMSCs on day 14; two-sided t-test, significance level FDR < 5%, n = 3. d, Pie chart illustrating differential features per channel and measurement class comparing siCOBLL1 and the non-targeting siRNA control on day 14. e,f, Spatial intensity distribution of AGP informative for the actin cytoskeleton in the centre of the cytoplasm: Cytoplasm_RadialDistribution_FracAtD_AGP_1of4 (e) and juxtaposed to the plasma membrane (Cytoplasm_RadialDistribution_RadialCV_AGP_4of4) (f). Two-sided t-test, n = 3 biologically independent experiments, median ± 95% confidence interval (CI). g, Representative images of COBLL1 knockdown and non-targeting siRNA control at days 0 and 14 of differentiation. COBLL1 staining: anti-COBLL1 primary and donkey anti-rabbit IgG H&L secondary antibodies. Actin staining: phalloidin–Atto 565. Olympus FLUOVIEW FV1000 CLSM Inverse microscope (40× magnification). Images were processed with Image J. Dotted square: image zoom-in. h,i, Texture of BODIPY stain (Cells_Texture_Correlation_Lipid_10_01) (h) and granularity measures of the BODIPY stain (Cells_Granularity_3_BODIPY) (i) in siCOBLL1 knockdown and non-targeting siRNA; two-sided t-test, n = 3 biologically independent experiments, median ± 95% CI.j, Oil Red O staining in SGBS adipocytes after stable lentiviral COBLL1 knockdown (shCOBLL1) versus empty vector control (shEV); scale bar, 15 mm. k, GPDH metabolic activity in differentiated shCOBLL1 compared to non-targeting siRNA adipocytes. Paired Student’s t-test, median ± 95% CI, n = 3. l, Basal and insulin-stimulated ³H-2-DG uptake in differentiated shCOBLL1 compared to shEV adipocytes. One-way analysis of variance (ANOVA) with Tukey’s honestly significant difference (HSD) test, median ± 95% CI, n = 3 biologically independent experiments. m, Basal and isoproterenol-stimulated lipolysis rate measured using glycerol release in differentiated shCOBLL1 compared to shEV adipocytes. One-way ANOVA with Tukey’s HSD test, median ± 95% CI, n = 3 biologically independent experiments. n, Western blots for lipolysis-relevant proteins assayed in basal or isoproterenol/IBMX-stimulated shCOBLL1 versus shEV adipocytes (n = 2).

Fig. 4 ∣. The rs6712203 MONW risk haplotype affects actin remodelling in adipocytes and adipocyte lipid storage capacity.

a, Schematic of adipocyte differentiation and LipocyteProfiling of subcutaneous AMSCs derived from TT (n = 7) and CC (n = 6) allele carriers of rs6712203 using LipocyteProfiler. b,d, Differences in morphological profiles between TT (n = 7) and CC (n = 6) allele carriers at day 14 in subcutaneous (b) and visceral (d) AMSCs (multi-way ANOVA, significance level FDR < 5%). c, Pie chart illustrating non-redundant differential features per channel and class of measurement at day 14 of subcutaneous adipocyte differentiation in rs6712003 homozygous risk carriers compared to non-risk carriers. e,f, Spatial intensity distribution of AGP in the centre of the cytoplasm near the nucleus in subcutaneous adipocytes derived from TT (n = 7), TC (n = 14) and CC (n = 6) carriers of rs6712203 (Cytoplasm_RadialDistribution_ FracAtD_AGP_1of4) (e) andjuxtaposed to the plasma membrane (Cytoplasm_ RadialDistribution_RadialCV_AGP_4of4) (f) throughout differentiation. Multi-way ANOVA; data represent the median ± 95% CI. g,h, Lipid droplet count (Cells_Children_LargeLipidObjects_Count) (g) and intensity of BODIPY stain (Cells_Intensity_IntegratedIntensity_Lipid) (h) throughout differentiation. Multi-way ANOVA; data represent the median ± 95% CI.

Fig. 5 ∣. Cobll1-deficient mice are leaner and display metabolically dysfunctional phenotypes.

a, Schematic of differentiation and LipocyteProfiling at three time points (days 0, 2 and 10) of AMSCs derived from Cobll1^−/− mice (n = 3) and WT mice (n = 4). b, Morphological profiles of the AMSCs of Cobll1^−/− mice compared to the AMSCs of WT mice at day 10; two-sided t-test, significance level FDR ≤ 5%. c, Pie chart illustrating non-redundant differential features per channel and class of measurement comparing the AMSCs of Cobll1^−/− and WT mice at day 10 of differentiation. d–g, Lipid droplet count (Cells_Children_LargeLipidObjects_Count) (d), granularity of BODIPY stain (Cells_Granularity_3_BODIPY) (e), intensity of BODIPY stain (Cells_Intensity_IntegratedIntensity_Lipid) (f) and texture of actin cytoskeleton (Cytoplasm_Texture_Entropy_AGP) (g) at day 10 of differentiation. Two-sided t-test; data represent the median ± 95% CI. h, Oil Red O staining of differentiated murine AMSCs. i, GPDH activity of differentiated murine AMSCs was assessed by measuring the decrease in NADH at 340 nm. Data represent the mean ± s.e.m. *P < 0.05 compared to the WT group. j, Representative photograph of 14-week-old WT and Cobll1^−/− mice fed a normal chow. The yellow dotted lines delineate perigonadal WAT. k, Mouse body weight across time. Data are expressed as the mean ± s.e.m. *P < 0.05, **P < 0.01 and ***P < 0.001 compared to the WT group. l, Body composition (fat mass/body weight). Data are expressed as the mean ± s.e.m. *P < 0.05 compared to the WT group. m, Body length measurements of WT and Cobll1^−/− mice (n = 6). Data are expressed as the mean ± s.e.m. *P < 0.05 compared to the WT group. NS, not significant. n, BMD analysed by dual energy X-ray absorptiometry (DEXA). Data are expressed as the mean ± s.e.m. *P < 0.05 compared to the WT group. o, IPGTT in WT and Cobll1^−/− and Cobll1^+/− mice. The inset graph shows the area under the curve (AUC) of the blood glucose concentration levels measured during IPGTT. Data represent the mean ± s.e.m. **P < 0.01 compared to the WT group. In i,k–o statistical significance was determined by Student’s t-test. In l–o the experiment was repeated independently three times with similar results.