Single nucleus RNA-sequencing integrated into risk variant colocalization discovers 17 cell-type-specific abdominal obesity genes for metabolic dysfunction-associated steatotic liver disease

Abstract

Background: Abdominal obesity increases the risk for non-alcoholic fatty liver disease (NAFLD), now known as metabolic dysfunction-associated steatotic liver disease (MASLD).

Methods: To elucidate the directional cell-type level biological mechanisms underlying the association between abdominal obesity and MASLD, we integrated adipose and liver single nucleus RNA-sequencing and bulk cis-expression quantitative trait locus (eQTL) data with the UK Biobank genome-wide association study (GWAS) data using colocalization. Then we used colocalized cis-eQTL variants as instrumental variables in Mendelian randomization (MR) analyses, followed by functional validation experiments on the target genes of the cis-eQTL variants.

Findings: We identified 17 colocalized abdominal obesity GWAS variants, regulating 17 adipose cell-type marker genes. Incorporating these 17 variants into MR discovers a putative tissue-of-origin, cell-type-aware causal effect of abdominal obesity on MASLD consistently with multiple MR methods without significant evidence for pleiotropy or heterogeneity. Single cell data confirm the adipocyte-enriched mean expression of the 17 genes. Our cellular experiments across human adipogenesis identify risk variant -specific epigenetic and transcriptional mechanisms. Knocking down two of the 17 genes, PPP2R5A and SH3PXD2B, shows a marked decrease in adipocyte lipidation and significantly alters adipocyte function and adipogenesis regulator genes, including DGAT2, LPL, ADIPOQ, PPARG, and SREBF1. Furthermore, the 17 genes capture a characteristic MASLD expression signature in subcutaneous adipose tissue.

Interpretation: Overall, we discover a significant cell-type level effect of abdominal obesity on MASLD and trace its biological effect to adipogenesis.

Funding: NIH grants R01HG010505, R01DK132775, and R01HL170604; the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant No. 802825), Academy of Finland (Grants Nos. 333021), the Finnish Foundation for Cardiovascular Research the Sigrid Jusélius Foundation and the Jane and Aatos Erkko Foundation; American Association for the Study of Liver Diseases (AASLD) Advanced Transplant Hepatology award and NIH/NIDDK (P30DK41301) Pilot and Feasibility award; NIH/NIEHS F32 award (F32ES034668); Finnish Diabetes Research Foundation, Kuopio University Hospital Project grant (EVO/VTR grants 2005-2021), the Academy of Finland grant (Contract no. 138006); Academy of Finland (Grant Nos 335443, 314383, 272376 and 266286), Sigrid Jusélius Foundation, Finnish Medical Foundation, Finnish Diabetes Research Foundation, Novo Nordisk Foundation (#NNF20OC0060547, NNF17OC0027232, NNF10OC1013354) and Government Research Funds to Helsinki University Hospital; Orion Research Foundation, Maud Kuistila Foundation, Finish Medical Foundation, and University of Helsinki.

Keywords: Abdominal obesity; Colocalization; Expression quantitative trait loci (eQTL); Genome-wide association study (GWAS); Metabolic dysfunction-associated steatotic liver disease (MASLD); Single-nucleus RNA-sequencing (snRNA-seq); Waist-hip ratio adjusted for body mass index (WHRadjBMI).

Fig. 1

Single-nucleus RNA-sequencing (snRNA-seq) and colocalization analyses identify adipose cell-type marker genes underlying abdominal obesity. a, UMAP illustration of 11 cell-type clusters in snRNA-seq data of subcutaneous adipose tissue from 8 individuals with obesity in the KOBS cohort. ASPC indicates adipose stem and progenitor cells; HPC, hematopoietic stem cell; LymphEndo, lymphatic endothelial cells; and UMAP, Uniform Manifold Approximation and Projection. b, Main adipose cell-type proportions estimated in the 262 bulk RNA-seq data using Bisque. The box shows the 25th and 75th percentiles, the centre line shows the medians, and the whiskers extend to the 5th and 95th percentiles. c, Cis-eQTL effect for PPP2R5A in the subcutaneous adipose tissue. Boxplots show association between the genotypes of rs10779574 and the normalised expression of PPP2R5A from the adipose bulk RNA-seq data of the 262 individuals in the KOBS cohort. The box shows the 25th and 75th percentiles, the centre line shows the medians, and the whiskers extend to the 5th and 95th percentiles. d, Comparison of WHRadjBMI GWAS and adipose cis-eQTL SNPs (left panel) and regional overview of WHRadjBMI GWAS (top right panel) and adipose cis-eQTL loci (bottom right panel) demonstrates a significant colocalization of the WHRadjBMI GWAS and adipose cis-eQTL SNP rs10779574, targeting the gene PPP2R5A. The axes show the −log₁₀ of p-values from the GIANT and UK Biobank WHRadjBMI GWAS meta-analysis and −log₁₀ of p-values from the subcutaneous adipose cis-eQTL analysis in the KOBS cohort (n = 262). The colocalized cis-eQTL SNP, rs10779574, is represented by a purple diamond. Colors represent LD (r²) with colocalized cis-eQTL SNP. Chr indicates chromosome; eQTL, expression quantitative trait locus; GWAS, genome-wide association study; LD, linkage disequilibrium; Mb, mega base; P, p-value; SNP, single nucleotide polymorphism; and UMAP, Uniform Manifold Approximation and Projection.

Fig. 2

Mendelian randomization (MR) demonstrates a putative causal effect of WHRadjBMI on MASLD using multiple MR methods. Significant putative causal effects of WHRadjBMI on MASLD are demonstrated by multiple MR methods using adipose cell-type-aware WHRadjBMI GWAS cis-eQTL variants as instrumental variables (IVs). Forest plot of MR analyses shows causal estimates (betas) with 95% confidence intervals and p-values computed using cML-MA-BIC, inverse variance weighting, MR-PRESSO and weighted median. No outlier SNPs were detected by the MR-PRESSO outlier test in the MR analysis, and no significant evidence (p-value>0.05) for horizontal pleiotropy was observed by the MR-PRESSO global test. CI indicates confidence interval; IVW, inverse variance weighting; MASLD, metabolic dysfunction-associated steatotic liver disease; Padj, p-value adjusted for multiple testing using Bonferroni; and WHRadjBMI, waist-to-hip ratio adjusted for body mass index.

Fig. 3

Nuclei expression of the 17 adipose cell-type-aware abdominal obesity genes shows preferential average expression in adipocytes, and longitudinal expression of the 10 adipocyte and ASPC marker genes changes significantly across human adipogenesis. a,b, The module scores of the 17 adipose cell-type-aware GWAS cis-eQTL target genes in the subcutaneous adipose snRNA-seq data (n = 8) show preferential average expression of the genes in the adipocyte cell-type. Module scores are calculated using the AddModuleScore function in Seurat as an average expression level of the 17 adipose cell-type-aware GWAS eQTL target genes while subtracting the aggregate expression level of the control feature sets for each nucleus. a, UMAP illustration of the adipose cell-type clusters where each dot represents a nucleus coloured by the module score of the 17 abdominal obesity genes. The p-value was computed using the Wilcoxon rank sum test to evaluate the difference in module scores of adipocyte vs non-adipocyte nuclei. b, Dot plot shows a higher average module score in the adipocytes. The size of each dot represents the percent of cells with a module score >0 in each cell-type and the colors represent an average module score for each cell-type. c, Human adipogenesis experiment shows longitudinal expression changes of the 10 adipocyte and ASPC marker genes during differentiation of human primary preadipocytes to adipocytes. Human primary preadipocytes were differentiated for 14 days and RNAs were collected at 6 time points for bulk RNA-sequencing. The genes were grouped into 4 distinct clusters based on the high probability of cluster assignment in their longitudinal expression trajectories during adipogenesis, detected using DPGP. Colors represent expression of the 10 adipocyte and ASPC marker genes from the 17 abdominal obesity genes quantified by bulk RNA-sequencing, and counts were normalised and scaled using ImpulseDE2. Gene-wise expression trajectory fits were obtained by implementing the impulse model, and the longitudinal differential expressions were evaluated using ImpulseDE2. ASPC indicates adipose stem and progenitor cells; HPC, hematopoietic stem cell; LymphEndo, lymphatic endothelial cells; and UMAP, Uniform Manifold Approximation and Projection.

Fig. 4

Changes in adipose marker gene expression during human adipogenesis are reflected in regional chromatin accessibility at the cis-eQTL SNP sites. a,d, Comparison of WHRadjBMI GWAS and adipose cis-eQTL SNPs demonstrates a significant colocalization of the WHRadjBMI GWAS and adipose cis-eQTL SNP rs2509963 targeting the gene AHNAK (a) and rs6866204 targeting the gene SH3PXD2B (d). The axes show the −log₁₀ of p-values from the GIANT and UK Biobank WHRadjBMI GWAS meta-analysis and −log₁₀ of p-values from the subcutaneous adipose cis-eQTL analysis in the KOBS cohort (n = 262). The colocalized cis-eQTL SNPs are represented by a purple diamond. Colours represent LD (r²) with colocalized cis-eQTL SNP. b,e, Bulk RNA-seq data of differentiating preadipocytes collected at 6 time points show significant differential expression (adj. p-value<0.05) of the adipose cell-type marker genes AHNAK (b) and SH3PXD2B (e) longitudinally during adipogenesis. c,f, Longitudinal bulk ATAC-seq data of differentiating preadipocytes collected at 6 time points show significant differentially accessible (adj. p-value < 0.05) chromatin regions that include the colocalized WHRadjBMI GWAS cis-eQTL SNPs targeting the adipose cell-type marker genes AHNAK (c) and SH3PXD2B (f). Gene-wise expression and chromatin accessibility trajectory fits were obtained by implementing the impulse model, and the longitudinal differential expressions were evaluated using ImpulseDE2. Significance of the longitudinal differential expressions and differential chromatin accessibility were measured using ImpulseDE2 and corrected for multiple testing using Bonferroni corrected p-value < 0.05. Adj. p indicates adjusted p-value; chr, chromosome; eQTL, expression quantitative trait locus; GWAS, genome-wide association study; LD, linkage disequilibrium; Mb, mega base; P, p-value; and SNP, single nucleotide polymorphism.

Fig. 5

Oil Red O (ORO) lipid staining reveals altered lipid accumulation with the knockdown of PPP2R5A and SH3PXD2B in human SGBS preadipocytes during adipogenesis. siRNA-mediated knockdown (KD) of two abdominal obesity genes PPP2R5A and SH3PXD2B in the differentiating human SGBS preadipocytes shows disruption in lipid accumulation. The cells were stained with ORO for each condition at three time points, and the intensity of the ORO staining was quantified by measuring the absorbance of 492 nm wavelength light and normalizing to the cell number. a, Cell images of differentiating SGBS preadipocytes stained with ORO (red color) taken at 3 time points for each condition using the EVOS Core XL microscope at 20x zoom. Rows indicate experimental conditions as follows: non-transfected controls, scrambled siRNA control for PPP2R5A (60 nM), scrambled siRNA control for SH3PXD2B (150 nM), PPP2R5A knockdown (60 nM), and SH3PXD2B knockdown (150 nM). Columns indicate the number of days from initiating the differentiation of preadipocytes: 0-day (baseline), 2-days (2D), and 7-days (7D). b, Relative ORO intensity in each knockdown condition at each time point compared to their respective scrambled controls. The average ORO stain intensity from 2 to 4 biological replicates each with 3 technical replicates for each condition (colors) at each time point (dots) with ±standard deviation (error bars) was compared (fold change) to the average ORO stain intensity of the non-transfected controls at the respective time points. Significant differences in ORO intensities between the knockdown and scrambled control samples are shown (∗, p-value < 0.05) by the t-test.

Fig. 6

Knockdown of PPP2R5A and SH3PXD2B in human SGBS preadipocytes identify altered expression of key adipogenesis genes. a,c, Results of the differential expression (DE) analysis using the bulk RNA-sequencing data (see Methods) from the siRNA-mediated PPP2R5A (a) and SH3PXD2B (c) knockdown samples compared to the scrambled control samples. The top 30 significantly upregulated (blue) and 30 significantly downregulated (orange) genes by log fold-change in knockdowns compared to the respective scrambled controls (logFC) at the 7-day (7D) time point and the respective knockdown genes (gold) are shown for each time point sorted by the logFC at the 7D time point. All significant DE genes for each KD (passing multiple testing correction using FDR <0.05) are shown in the Supplementary Tables S10 and S11. b,d, Average expression, in counts per million (CPM), of the knockdown genes and selected key adipogenesis, fat storage regulator, and adipose tissue function genes are shown for each knockdown condition and scrambled controls at each time point. Dots represent average CPMs from 2 to 4 technical replicates with error bars indicating ± standard deviation. Time point is represented by the x-axis and the average expression (CPM) by the y-axis. The colours represent average expression in the knockdown vs respective scrambled controls. Significant differences in the average expression of each gene between the knockdown and scrambled controls samples at each time point are shown (∗, p-value<0.05; ∗∗, p-value<0.01; ∗∗∗, p-value < 0.001) by the t-test. FDR indicates false discovery rate; and NS, non-significant (FDR ≥ 0.05).

Fig. 7

Expression profiles of the 17 abdominal obesity genes and DGAT2 in the adipose tissue demonstrate differences based on the liver phenotype status. a,b,c, DGAT2 expression quantified in the adipose tissue and liver from the bulk RNA-sequencing data of 262 individuals with obesity in the KOBS cohort shows higher expression in the adipose tissue when compared to the liver. Within the adipose tissue, DGAT2 expression is significantly higher in the individuals with obesity and diagnosed with steatosis (n = 154) (a), fibrosis (n = 115) (b), or metabolic associated steatohepatitis (MASH) (n = 81) (c) using the liver histology-based assessment when compared to the controls with obesity and healthy livers (n = 86 for all comparisons). In contrast, the liver bulk RNA-sequencing data show no significant difference in the liver DGAT2 expression based on the liver metabolic dysfunction-associated steatotic liver disease (MASLD) diagnosis groups. The x-axis represents the tissue where DGAT2 is expressed and the y-axis represents DGAT2 expression measured in Trimmed Mean of M-values (TMM) normalised counts per million (CPM). The colours represent the liver conditions (0 for the controls and 1 for steatosis, fibrosis, and MASH). Significant differences in the DGAT2 expression between individuals with a liver disease vs controls were evaluated in each tissue by the Wilcoxon rank sum test using TMM normalised, log-transformed CPMs adjusted for technical factors and cell-type proportion estimates of their respective tissues (see Methods). d, The first principal component of the subcutaneous adipose expression of the 17 abdominal obesity genes is lower in the individuals with obesity and liver steatosis (n = 154) when compared to the controls with obesity and healthy livers (n = 86). The x-axis represents the hepatic liver diagnosis status (0 = controls and 1 = cases) and the y-axis represents the first principal component (PC1) of the principal component analysis (PCA) performed using the TMM normalised, log-transformed CPMs of the 17 abdominal obesity genes after adjusted for technical factors and adipose cell-type proportion estimates (see Methods). The box shows the 25th and 75th percentiles, the centre line shows the medians, and the whiskers extend to the 5th and 95th percentiles. The significance of the difference in PC1 of the adipose expression of the 17 abdominal obesity genes between the 2 groups was evaluated using the Wilcoxon rank sum test. a-d, The box shows the 25th and 75th percentiles, the centre line shows the medians, and the whiskers extend to the 5th and 95th percentiles, and the significant differences between the two groups are annotated (∗, p-value < 0.05; ∗∗, p-value<0.01; ∗∗∗, p-value < 0.001). NS indicates non-significant (p-value ≥ 0.05).