Integration of single cell omics with biobank data discovers trans effects of SREBF1 abdominal obesity risk variants on adipocyte expression of more than 100 genes

Study cohorts

Genotype quality control and imputation in the KOBS and RYSA cohorts

We genotyped DNAs of KOBS participants using the Infinium Global Screening Array-24 v1 (Illumina). We performed quality control (QC) on the genotype data using PLINK v1.9²¹ by excluding individuals with missingness >2% and removing unmapped, strand ambiguous, and monomorphic SNPs in addition to variants with missingness >2% and Hardy-Weinberg Equilibrium (HWE) p-value<10^-6. We further imputed biological sex using the ‘--sex-check’ function in PLINK v1.9²¹ and cross-checked with the reported sex of each individual.

We performed genotype imputation against the HRC reference panel version r1.1 2016²² on the Michigan imputation server after removing duplicate variants and variants with allele mismatch with the reference panel. Strand flips or allele switches needed to match the reference panel were performed on the server before haplotype phasing using Eagle v2.4²³ and genotype imputation using minimac4²⁴. We performed QC on the imputed genotype data by removing SNPs with imputation score R²<0.3 and HWE p-value<10^-6 for the downstream analyses.

We genotyped the RYSA participants’ DNAs using the Infinium Global Screening Array-24 v1 (Illumina). Genotype data QC, imputation, and post-imputation QC were performed for the RYSA cohort as described above.

Nuclei isolation in VAT biopsies from the KOBS cohort

To isolate nuclei from the snap-frozen VAT biopsies from the KOBS participants, we first combined and minced the VAT samples in a petri dish over dry ice and immediately transferred the minced tissue into a 500 μl chilled 0.1X Lysis Buffer consisting of 10 mM Tris-HCl, 10 mM NaCl, 3 mM MgCl₂, 0.1% Tween-20, 0.1% IGEPAL CA-630, 0.01% Digitonin, 1% BSA, 1 mM DTT, and 1 U/μL RNase inhibitor. After a 15-minute incubation period, the lysate mixed with a 500 μl chilled Wash Buffer containing 10 mM Tris-HCl, 10 mM NaCl, 3 mM MgCl₂, 1% BSA, 0.1% Tween-20, 1 mM DTT, and 1 U/μL RNase inhibitor was passed through a 70 μm Flowimi Cell Strainer into a 2 ml tube. Nuclei were centrifuged at 500 rcf for 5 minutes at 4°C and the supernatant was removed without disrupting the nuclei pellet, followed by a resuspension in a 1 ml chilled Wash Buffer. Nuclei were passed through a 40 μm Flowimi Cell Strainer into a 2 ml tube and centrifuged at 500 rcf for 5 minutes at 4°C. We removed the supernatant without disrupting the nuclei pellet and resuspended in a 30 μl chilled Diluted Nuclei Buffer containing 1X Nuclei Buffer (10x Genomics), 1 mM DTT, and 1 U/μl RNase inhibitor. We measured the concentration and overall quality of the nuclei using Countess II FL Automated Cell Counter after staining with trypan and DAPI and used the Single Cell Multiome ATAC + Gene Expression Reagent Kit (10x Genomics) for the joint snRNA- and snATAC-seq library construction. We analyzed the quality of cDNA and libraries using Agilent Bioanalyzer and sequenced the libraries on an Illumina NovaSeq X Plus (snRNA) and Next seq 500 (snATAC) with a target sequencing depth of 400 million read pairs for both snRNA- and snATAC-seq.

Nuclei isolation in SAT biopsies of the KOBS and RYSA cohorts

We performed SAT snRNA-seq experiments on the snap-frozen SAT biopsies from the KOBS participants, as previously described¹¹. Briefly, we first pooled approximately 100 mg of each biopsy and isolated nuclei from the pooled biopsies as described earlier¹¹. Next, we measured the concentration and overall quality of the nuclei using Countess II FL Automated Cell Counter after staining with trypan and DAPI and used the Single Cell 3’ Reagent Kit v3.1 (10x Genomics) for the library construction. We analyzed the quality of cDNA and gene expression library using Agilent Bioanalyzer and sequenced the library on an Illumina NovaSeq SP with a target sequencing depth of 600 million read pairs.

In this study, we used the operation time-point SAT snRNA-seq data from 68 participants with obesity of the RYSA cohort, generated as part of the full RYSA cohort, by randomly pooling 8 SAT samples/batch for nuclei isolation. We isolated nuclei and constructed libraries for SAT snRNA-seq from each batch as described above and sequenced libraries from all batches together on an Illumina NovaSeq S4 with a target sequencing depth of 400 million read pairs per batch.

Processing of the SAT snRNA-seq data from the RYSA cohort

We aligned the raw snRNA-seq data in a FASTQ format file for each batch against the GRCh38 human genome reference and GENCODE v42 annotations²⁵ using STAR v2.7.10b²⁶ with the ‘-- soloFeatures GeneFull’ option to account for full pre-mRNA transcripts. We evaluated the quality of the raw and mapped snRNA-seq data using FastQC. Next, we used DIEM v2.4.0²⁷ to remove empty droplets and droplets with high amounts of ambient RNA. We set droplets with UMI<500 as debris, used k=50 for the initialization step with k-means clustering, and otherwise used the default parameters. After the clustering step, we removed clusters of highly contaminated droplets characterized by low average UMIs, low average number of unique genes detected (nFeatures), high percentage of reads mapped to the mitochondrial genome (%mito), and high number of mitochondrial and ribosomal genes as top expressed features. Next, we used Seurat v4.3.0.1²⁸ to remove low-quality droplets with UMI<500, nFeatures<200, %mito>10%, and spliced read fraction≥75%, log-normalize gene counts using the ‘NormalizeData’ function, identify top 2,000 variable genes excluding mitochondrial and ribosomal genes using the ‘FindVariableFeatures’ function, scale the normalized gene counts to mean 0 and variance 1 using the ‘ScaleData’ function, perform principal component analysis (PCA) using the ‘RunPCA’ function, and cluster the nuclei with a standard Louvain algorithm, first 30 PCs, and a resolution of 0.5. For the remaining nuclei, we removed contaminated counts using DecontX²⁹ with the previously removed low-quality nuclei as the background and the Seurat cluster assignment as the ‘z’ and removed additional low-quality nuclei with UMI<500, UMI>30,000, nFeatures<200, and %mito>10% based on the remaining clean counts. To identify the originating participant of each nucleus, we used demuxlet from popscle software tool³⁰ with ‘-- min-MQ 30’ and high-quality imputed genotype data. We only included nuclei classified as a singlet and assigned the best matching participant to identify the originating participant for each nucleus. Next, we used DoubletFinder³¹ to identify and remove any remaining doublets. As DoubletFinder requires predicted number of doublets in each dataset, we performed pN-pK parameter sweeps on a subset of 10,000 nuclei to select pN of 0.25 and the most optimal pK value maximizing mean-variance normalized bimodality coefficient for each batch, as previously recommended²⁷.

After performing the QC, we used Seurat v4.3.0.1²⁸ to merge all remaining high-quality droplets from the batches and subset for nuclei originating from operation time-point samples of the 68 participants from the RYSA cohort included in this study. In the subset, we kept only the genes with at least 3 raw counts in at least 3 nuclei²⁹ and performed gene count normalization, variable gene identification, data scaling, and PCA, as described above. To account for variation in gene expression driven by batch effect, we used Harmony v1.0.3³² to integrate on batch and clustered the nuclei with a standard Louvain algorithm, first 30 reductions from Harmony, and a resolution of 0.5. Cell-type annotations were performed using SingleR v1.8.1³³, as described above, with the SAT single-cell and snRNA-seq data from the previously published adipose tissue atlas by Emont et al. as a reference³⁴.

Processing of the SAT VAT snRNA-seq data from the KOBS cohort

We aligned the VAT and SAT snRNA-seq data from the seven KOBS biopsies to the GRCh38 human genome reference with Ensembl annotations using CellRanger-arc -count v2.0.0³⁵ and CellRanger -count v6.1.1³⁵ for VAT and SAT, respectively. As the VAT sequencing from the seven KOBS VAT biopsies was generated using the 10x Genomics Multiome platform, we simultaneously profiled RNA and ATAC in each cell. We subset the data to include the nuclei passing the QC for both RNA and ATAC, as recommended by CellRanger-arc, while only including the VAT snRNA-seq data in the downstream analyses of this study. The KOBS SAT snRNA-seq data were generated as a single modality and processed using the same reference (GRCh38 human genome) but with the base software for CellRanger. Similarly to VAT, we only included the nuclei from SAT that passed the recommended filtering from CellRanger. For the remaining SAT and VAT nuclei, we separately further filtered out those containing ambient reads using DIEM²⁷ through the removal of droplets from a debris cluster that reflected high amounts of background RNA or low levels of nuclear RNA. To remove genetic doublets and demultiplex the KOBS VAT Multiome samples and SAT snRNA-seq samples back to their individuals of origin, we ran demuxlet³⁰ against the imputed genotype data of the individuals. We removed all SAT and VAT cells that were not classified as singlets. We then ran DecontX²⁹ to remove contaminated reads within each droplet and filtered out nuclei with high levels of ambient RNA, keeping those with UMIs over 200 and mitochondrial read percentage<10%. These QC processes resulted in 3,216 nuclei for VAT and 3,516 nuclei for SAT from the same individuals.

Integration and clustering of SAT and VAT snRNA-seq data from the KOBS cohort

To find differences between SAT and VAT cell-type level gene expression from the KOBS biopsy samples, we performed integration between the tissues using Seurat v4.3.0²⁸. The count data for nuclei from each tissue were first log-normalized using the NormalizeData function of Seurat²⁸, with the default scaling factor of 10,000. Using these normalized count data, the top 2000 variable genes per tissue were calculated using the FindVariableFeatures function. Afterwards, we integrated the two tissues using canonical correlation analysis (CCA) with the IntegrateData function of Seurat²⁸. Following integration, the normalized read counts were scaled with a mean of 0 and variance of 1. In the integrated samples with a total of 6,732 cells from across the SAT and VAT snRNA-seq data, we performed PCA and identified clusters using 30 PCs and a resolution of 0.8.

SnRNA-seq cell-type assignment and unique marker gene identification in the KOBS SAT and VAT biopsy samples

We used SingleR v2.0.0³³ with a previously published human single-cell adipose atlas³⁴ as a reference to annotate the data at the individual cell level³⁶. We ran the FindAllMarkers function in Seurat²⁸ using only.pos = TRUE and min.cells.group=50 to identify marker genes of each cell-type. We only included cell-types that consisted of at least 50 nuclei from the integrated tissues. P-values were adjusted for multiple testing using the Bonferroni approach. To identify marker genes unique per cell-type, we removed genes detected as marker genes (padj<0.05) for multiple cell-types. These marker genes are shown in Supplementary Table 1.

Tissue-specific marker gene identification

To find differences in the expression profiles of the three most prevalent SAT and VAT cell-types (adipocytes, macrophages, and adipose stem and progenitor cells (ASPC)) in the KOBS cohort, we identified three separate sets of cell-type marker genes. These three sets comprised cell-type marker genes specific to SAT (MGSS), VAT (MGSV), and shared between the tissues (MGSBT). We identified tissue-specific markers using the FindMarkers function of Seurat²⁸, using a log₂ fold change threshold of ≥0, to compare the SAT and VAT data of each cell-type. In this tissue-specific marker gene identification, we only included the genes identified as unique markers for each cell-type above. We adjusted p-values for multiple testing using FDR<0.05 for the total number of cell-type-specific markers from above. Using this adjustment method, we identified genes that are not only cell-type-specific but also show differential expression (DE) between VAT (MGSV) and SAT (MGSS). We also identified a third group of marker genes for each cell-type, i.e., those that showed no significant DE (FDR>0.05) between SAT and VAT. We classified these as marker genes shared between the SAT and VAT tissues (MGSBT). Overall, this resulted in 9 cell-type marker gene sets, i.e., three per cell-type (macrophages, adipocytes, and ASPC) with no overlap between them. These 9 sets are listed in Supplementary Table 2.

Identification of enriched functional pathways

After identifying the 9 sets of cell-type marker genes, as described above, we ran WebGestalt³⁷ to functionally characterize each set, using all genes with non-zero expression in at least two cells in the cell-type³⁴ irrespective of the tissue origin, as the background. For each set, we tested for significant (FDR<0.05) overrepresentation of genes from Gene Ontology (GO) biological processes, cellular components, and molecular functions. We then identified the top pathway genes (i.e. the genes that appeared in >1 significantly enriched pathways) by ranking each gene by the number of times it appeared in the identified pathways.

GWAS enrichment analysis of the cell-type marker genes

We assessed the 9 sets of tissue-specific and tissue-shared cell-type marker genes for significant enrichments of genetic associations with BMI, WHRadjBMI, T2D, and MASLD using MAGENTA v2.4 (Meta-Analysis Gene-set Enrichment of variaNT Associations)³⁸.

We first downloaded publicly available European-only GWAS summary statistics for T2D³⁹, while for WHRadjBMI and BMI, we used the UKB GWAS data we generated, as described above. For every outcome, we then used the summary statistics to generate association scores for each gene within the gene set based on the p-values of all variants within a 500kb upstream and downstream window. Gene score-cutoff tests using the 75th percentiles of the scores of all genes as cutoffs were performed to evaluate enrichment. While MAGENTA³⁸ filters out genes near one another, we retained all genes per set in our analyses due to our small gene set sizes.

Construction of marker gene -based regional and genome-wide PRSs for abdominal obesity, BMI, and type 2 diabetes in UKB

We built both marker gene -based regional and genome-wide PRSs for waist-hip ratio adjusted for BMI (WHRadjBMI), body mass index (BMI), and type 2 diabetes (T2D) for all unrelated Europeans from UKB¹⁸. In our analysis, WHRadjBMI was used as a well-established proxy for abdominal obesity¹¹. The genome-wide PRSs for these traits were constructed using all SNPs (see below), whereas the regional PRSs for each trait were created using the cis regional SNPs (gene+/-500kb) of each of the autosomal marker gene sets separately (MGSBT, MGSS, and MGSV) for the three main cell-types (adipocytes, macrophages, and ASPCs). We excluded the set that contained shared SAT/VAT macrophage marker genes as it comprised only 9 genes, and, therefore, it remained too small for reliable regional PRS building given the technical issues with overfitting in PRS⁴⁰. As WHRadjBMI is a highly sex-specific trait¹¹, we also built the regional WHRadjBMI PRSs for males (n= 88,988) and females (n=104,614) separately. Thus, in addition to the genome-wide PRSs, we constructed 24 regional PRSs for WHRadjBMI and 8 regional PRSs for T2D and BMI, respectively.

Construction of WHRadjBMI and BMI PRSs

To construct the WHRadjBMI and BMI PRSs, we randomly selected 50% of the unrelated Europeans from UKB (n=195,863) to generate GWAS summary statistics as the base data for the PRS, and applied the PRS model in the remaining 50% of the cohort (i.e., the target group). For the GWASs, the data were filtered to only retain variants with MAF >1% and INFO>0.8. We also removed individuals with >1% of their genotypes missing or extreme heterozygosity, as well as variants missing in >1% of the subjects or violating Hardy-Weinberg equilibrium from the genotype data target group⁴¹. We conducted the GWASs of WHRadjBMI and BMI using linear-mixed model approach of BOLT-LMM v2.3.6²³, with age, age², sex, the top 20 genetic PCs, testing center, and genotyping array as covariates, where we applied a rank-based inverse normal transformation for each outcome. For WHRadjBMI, given the sex-specificity of the outcome^11,42,43, we normalized twice, first in all individuals, and then for males and females separately. Next, we fit the PRS models and computed the PRSs for WHRadjBMI and BMI for the individuals in the target group using the split-validation mode from lassosum⁴⁴ (n=193,602) and the filtered GWAS summary statistics as the base data for the PRS models.

Construction of T2D PRSs

To construct the T2D PRSs, we used the publicly available T2D GWAS summary statistics of Europeans without UKB from the DIAGRAM Consortium⁴⁵ as the base data. We built the PRS models using the target group genotype data employing PRScs⁴⁶, with the provided LD matrix from the Europeans in the 1000 Genomes Project as the LD reference, and applied the models to the target group using PLINK.

Calculating the explained variance in the trait by the PRSs

We calculated the incremental variance explained (R²) in the trait against a null model containing the covariates of age, age², the top 20 genetic PCs, testing center, genotyping array, and sex using lassosum⁴⁴. For T2D, we calculated a delta AUC, similarly to the incremental variance from above.

We evaluated the significance of each regional PRS by ranking the observed incremental R² or delta AUC against the incremental R²/delta AUC of 10,000 PRSs similarly built from the cis regional variants of random gene sets of the same size as in the regional PRS set, using all expressed genes (counts>=1 in at least 2 nuclei)³⁴ in the cell-type of interest as a background. The P-value was defined as the number of permutations in which the R²/delta AUC is larger than the calculated PRS divided by 10,000.

Partitioned heritability assessment of PRS-enriched gene sets

We performed a partitioned heritability analysis with LD Score regression (LDSC)^47,48 to assess the genes with WHRadjBMI R² enriched PRSs for enrichment in the WHRadjBMI heritability relative to the genome, similarly as in Finucane et al.⁴⁸. Briefly, we used LDSC to estimate the LD scores from all SNPs in the genome, as well as from all SNPs residing with the cis regions of each PRS-enriched gene set. Scores were constructed in 76,758 randomly selected, unrelated Brits from UKB (35,257 males, 41,501 females) for computational efficiency, and as in Finucane et al.⁴⁸, we only included SNPs with MAF>5% for the analysis. We then used the same WHRadjBMI GWAS summary statistics as in the PRS analysis to compute the WHRadjBMI heritability of the genome and marker gene annotations and test each marker annotation for significant (p<0.05) heritability enrichment, defined as the proportion of heritability divided by the proportion of the genome from the annotation⁴⁸.

Longitudinal RNA- and ATAC-sequencing during human SAT primary preadipocyte differentiation

We previously performed a longitudinal adipogenesis experiment using human SAT primary preadipocytes and generated longitudinal bulk RNA- and ATAC-seq data on samples collected at the 0 day (0d), 1d, 2d, 4d, 7d, and 14d timepoints, with 4 isogenic replicates per timepoint, as described earlier^7,49. Peaks from the ATAC-seq were filtered to remove blacklisted regions and identify consensus peaks, as described in detail previously⁷.

Differential expression (DE) analysis across six adipogenesis time points

We examined the longitudinal expression patterns of the top pathway-enriched genes from the adipocyte MGSS set (n=43 genes) using ImpulseDE2 v0.99.10⁵⁰. We used the runImpulseDE2 function with default parameters and boolCaseCtrl=FALSE, boolIdentifyTransients=TRUE, and scaNProc=1. P-values were corrected for multiple testing using FDR<0.05. As the adipogenesis experiment was conducted in SAT preadipocytes, we only ran this DE analysis for genes in the adipocyte MGSS set.

Identification of longitudinally co-expressed clusters of adipocyte pathway genes and their regulators during human adipogenesis

To search for longitudinal co-expression patterns among the adipocyte pathway genes across human adipogenesis, we ran DPGP v0.1⁵¹ to cluster genes by their expression trajectories. We ran DPGP with the timepoint data of the adipogenesis experiment described above and only included the genes that we identified as significantly longitudinally DE (FDR<0.05) during human adipogenesis from ImpulseDE2⁵⁰ and expressed across adipogenesis (n=42).

Identification of co-accessible peaks during human adipogenesis within the cis region of SREBF1

To identify co-accessible peaks in the SREBF1 cis region during human adipogenesis, we used the longitudinal ATAC-seq data from the human SAT primary preadipocyte differentiation experiment described above. We first identified SAT adipocyte peaks within ±500 kb of SREBF1, a key SAT adipocyte marker gene. To determine temporal differential accessibility (DA) of these SAT peaks and identify clusters of longitudinally co-accessible peaks, we ran ImpulseDE2⁵⁰ as described above on the accessibility counts for the peaks, subset the data to the peaks passing an FDR threshold<0.05, and then clustered the DA peaks (n=120) using DPGP⁵¹, as described above. We only included the clustered peaks that passed a cluster assignment threshold of probability>0.9.

Identification of WHRadjBMI GWAS variants within the longitudinally co-accessible SAT adipocyte peaks in the SREBF1 cis region

We next investigated the 5 longitudinally co-accessible SAT adipocyte peaks discovered in the SAT adipogenesis experiment within the cis regions of SREBF1 for WHRadjBMI GWAS variants. Accordingly, we overlapped the positions of the genome-wide significant (p<5× 10^-8) WHRadjBMI GWAS variants using the summary statistics from a previously published extensive WHRadjBMI GWAS of the GIANT–UKB meta-analysis⁵² with the positions of the co-accessible SAT adipocyte peaks.

Identification of adipocyte MGSS genes by risk allele status of the WHRadjBMI GWAS variants in the SREBF1 region

We tested whether the risk allele carrier status of the seven WHRadjBMI GWAS variants residing in the longitudinally co-accessible SAT peaks in the SREBF1 region affects the SAT adipocyte expression of SREBF1 and other SAT adipocyte markers genes in trans using RYSA SAT adipocyte snRNA-seq data. Briefly, for each variant, we labeled each individual from the RYSA cohort with zero copies of the GWAS trait-increasing allele as “non-carriers” and those with one or two copies of the trait-increasing allele as “carriers,” and then labelled each individual cell in the cohort as coming from a “carrier” or a “non-carrier” individual. We subset the SAT snRNA-seq data to adipocytes and ran the FindMarkers function in Seurat²⁸ with logfc.threshold = 0 between cells labeled “carriers” and “non-carriers” while only testing the adipocyte MGSS gene set and adjusting the p-values for multiple testing using FDR<0.1. We repeated this for each of the seven SNPs. For each trans gene set, we built module scores using the AddModuleScore function in Seurat v4.3.0²⁸ in SAT adipocytes from the RYSA cohort to determine the average expression across all adipocytes between carriers and non-carriers.

Identification of adipogenesis open chromatin peaks by cell-type

To identify upregulated and shared peaks at the preadipocyte and adipocyte stages of human adipogenesis, we ran DESeq2⁵³ on the raw peak count data generated using ATAC-seq at the day 0 and 14 time points in the SAT adipogenesis experiment described above. Peaks with a log₂FC>0.1 and Bonferroni adjusted P (padj)<0.05 were classified as upregulated in preadipocytes, while the peaks with a log₂FC<-0.1 and padj<0.05 were classified as upregulated in adipocytes. The remaining peaks were classified as shared between the human preadipocytes and adipocytes.

Construction of abdominal obesity PRSs in UKB for the adipocyte marker genes affected by the seven regional SREBF1 variants in trans

To investigate the abdominal obesity risk of the regional variants in the trans gene sets, the adipocyte expression of which is up- or downregulated by the SREBF1 variant carrier status, we constructed partitioned PRSs using the cis regional SNPs of these trans gene sets. The PRSs were built separately for each up/downregulated trans gene set for the seven SNPs, resulting in 14 total unique PRS sets. As abdominal obesity risk differs by sex^42,43, we analyzed these partitioned abdominal obesity PRSs in all individuals, females, and males separately. Due to the small region sizes and to allow for fine-grained control of the SNPs comprising each PRS, we used a clumping and thresholding (C+T) approach with PLINK v1.9 for the PRS construction, with the same GWAS summary statistics for WHRadjBMI described above, as the base data. We partitioned the remaining 50%, not used for the GWAS, into two cohorts, a 30% test group (n=115,120) for performing the LD clumping and learning the optimal thresholding parameter, and a 20% validation group (n=76,758) to apply the learned PRS model and perform downstream analyses. We performed QC on the test and validation genotype data, similarly as described above.

To build the partitioned abdominal obesity PRSs, we first performed LD clumping on the full genotype data of the test groups with PLINK v1.9 using R² of 0.2 and window size of 250 kb, and built the PRSs from the independent SNPs passings p-value thresholds ranging from 5×10^-8 to 0.5 to empirically determine the p-value threshold that maximizes the incremental variance explained in WHRadjBMI in the test group. Next, for each trans gene set, we identified all C+T SNPs that landed in the adipocyte upregulated open chromatin peaks (described earlier), using the empirically determined p-value thresholds of p<0.05, p<0.005, and p<0.1 for all individuals, females, and males, respectively. To avoid direct SREBF1 cis effects on the abdominal obesity risk, we removed the seven SNPs in the SREBF1 region as well as any C+T SNPs in LD (R²>0.1) with any of these seven SNPs from these partitioned PRS analyses. We computed PRSs for the sets that contained >= 5 valid SNPs. As explained before, the incremental variance (R²) was calculated for each PRS and compared against 10,000 permuted PRSs built from the same numbers of C+T SNPs.

After computing individual level WHRadjBMI PRSs, we compared the polygenic risk, i.e., the PRS magnitudes, in the regional PRSs of the trans gene sets between the SREBF1 variant carriers and non-carriers. First, we classified individuals in the UKB PRS cohort into “carriers” or “non-carriers” in the same way as in the RYSA cohort. Next, we performed a two-sided Wilcoxon test between the WHRadjBMI PRSs of the carriers and non-carriers within each group while correcting for multiple testing using Bonferroni.

Study design

The summary of the study design is shown in Supplementary Figure 1. Briefly, by first identifying cell- and tissue-type expression profiles that drive functional differences between the two key human fat depots, subcutaneous adipose tissue (SAT) and visceral adipose tissue (VAT) and genetic enrichments of these SAT and VAT -specific genes for abdominal obesity, we discovered a transcription factor (TF) gene preferentially involved in numerous functional SAT pathways over VAT. We then further investigated the cell-type level regulatory mechanisms of this TF, and discovered that the regional abdominal obesity GWAS variants of this TF affect the adipocyte expression of tens of SAT adipocyte marker genes in trans in the independent RYSA cohort. Finally, we confirmed the identified trans effects of the TF variants by comparing the partitioned regional PRSs of the trans genes by their risk allele carrier status at the biobank level (Figure 1a).

• Download figure
• Open in new tab

Figure 1 Adipose single nucleus RNA-sequence data from the KOBS participants’ subcutaneous adipose tissue (SAT) and visceral adipose tissue (VAT) identify distinct sets of tissue-specific and shared cell-type marker genes.

(a) The study design visualized through cohort usage. The data from the individuals in the KOBS cohort were used to identify cell-type marker gene sets, while the data from the individuals in the RYSA cohort were used to identify trans gene sets by the carrier status of the SREBF1 risk alleles. Lastly, all polygenic risk scores were built using individuals from the UK Biobank. (b-c) The integrated SAT and VAT single nucleus RNA-seq data (n=6,732 nuclei) from 7 KOBS participants are visualized in a Uniform Manifold Approximation and Projection (UMAP) and colored by (b) fat depot of origin and (c) cluster cell-type identity, while mapping each cluster to one of the 16 identified cell-types. ASPC indicates adipose stem and progenitor cells; B, B cells; LEC, lymphatic endothelial cells; NK, natural killer cells; SMC, smooth muscle cells; and T, T cells. (d) A volcano plot depicts differential expression (DE) of the unique adipocyte marker genes by tissue. Genes are grouped into marker genes specific to SAT (MGSS), specific to VAT (MGSV), and shared between the SAT and VAT tissues (MGSBT) (see Methods), and then colored by these groups.

Marker gene analysis reveals distinct differences between SAT and VAT in cell-type level expression profiles

To elucidate cell-type level functional differences between SAT and VAT, the two main human fat depots¹, we first investigated SAT and VAT gene expression profiles at the cell-type resolution using SAT and VAT biopsies from the same individuals (Figure 1b-c). These dual tissue fat depot data from the same individuals, which reduce confounding by phenotypic differences between the donors of the SAT and VAT biopsies, enabled us to identify unique SAT and VAT cell-type marker genes, or genes that are differentially upregulated in only one cell-type, and also unique to each adipose tissue depot in their marker gene status (see Methods). Figure 1d shows the adipocyte cell-type markers distributed by tissue specificity. Supplementary Table 1 shows the identified unique cell-type marker genes for all SAT and VAT cell-types after integrating the single nucleus level data of the two tissues. As small differences often drive depot-specific gene and functional patterns⁴², we sought to study not only differences in expression between the depots but also shared transcriptional profiles, as this may elucidate the connection and shared characteristics between the tissues. Thus, from each cell-type marker set, we identified sets of unique marker genes specific to each tissue and shared across the two, resulting in 3 sets: marker genes specific to SAT (MGSS), marker genes specific to VAT (MGSV), and marker genes shared between tissues (MGSBT) (Supplementary Table 2). We focused our analyses on the three major shared cell-types between SAT and VAT, i.e., adipose stem and progenitor cells (ASPCs), adipocytes, and macrophages due to their functional pathway enrichments (see below) and highest number of nuclei available across the two tissues. The MGSS sets contained 244, 208, and 250 unique cell-type marker genes for ASPCs, adipocytes, and macrophages; the MGSV sets 131, 225, and 224; and the MGSBT sets 112, 155, and 9, respectively (Supplementary Table 2).

Shared and tissue-specific marker genes are enriched for genetic signals of obesity and related CMDs

As both adipose tissue depots are known to be impacted by obesity and related disorders^2,4,5,54, we first examined these marker gene sets for GWAS variant enrichments of obesity and related CMDs. Accordingly, we tested the variants in the cis regions (±500kb) of the three sets per cell-type in three cell-types, MGSS, MGSV, and MGSBT (nine total sets), for enrichment of BMI, abdominal obesity (using WHRadjBMI as a proxy), and type 2 diabetes (T2D) GWAS signals using GWAS summary statistics and the MAGENTA tool (see Methods). We observed significant enrichments (FDR<0.05) in the ASPC MGSBT for BMI, and in the adipocyte MGSS for T2D and WHRadjBMI (Table 1), while no enrichments were detected for any cell-types with the MGSV sets.

Next, we assessed each marker gene set for enrichment of polygenic risk for these obesity and related traits. We evaluated the predictive power of annotated polygenic risk scores (PRSs) for WHRadjBMI, BMI, and T2D from the cis regional variants of each set in the large UK Biobank (UKB) cohort. We first built PRSs for all individuals for each cell-type set and found that the regional PRS constructed for the adipocyte MGSS gene set was a significant predictor (p_R²<2.23×10^-308) of WHRadjBMI. Next, we followed up on these WHRadjBMI PRS findings using permutations. Of our 10,000 permutation scores, each similarly built from SNPs residing in the cis regions of the same number of randomly selected adipocyte expressed genes from either tissue, we observed only 0.82% to have an incremental R² greater than or equal to that of the adipocyte MGSS WHRadjBMI PRS (R²=1.20%; p_perm10,000=0.0082) (Figure 2; Supplementary Table 3). Due to this observed enrichment in variance explained of WHRadjBMI, which is known to be a sex-dimorphic trait^43,55, we also constructed these regional WHRadjBMI PRSs for females and males separately. We observed the enrichment in all individuals to be driven by females as the WHRadjBMI PRS constructed for the cis regional variants of the adipocyte MGSS genes was significantly enriched (R²=2.12%, p_perm10,000= 0.0101) in females, but not in males (Figure 2; Supplementary Table 3). In line with previous connections with visceral adipose tissue^56,57, we found that the variance explained for T2D was significantly enriched for the adipocyte MGSV gene set (Delta AUC=1.16%, p_perm1,000= 0.036). No other significant enrichments in variance explained were detected for BMI or T2D in our PRS analyses, and no other significant enrichments were detected for other cell-types for WHRadjBMI.

• Download figure
• Open in new tab

Figure 2 Regional polygenic risk scores (PRSs) constructed from the cis regional variants of the adipocyte MGSS gene sets show enrichment in variance explained for abdominal obesity in all individuals and females of the UK Biobank.

Lollipop plots visualize the predictive performance of each regional PRS relative to the 10,000 PRSs built using the same number of randomly selected adipocyte-expressed genes. Each horizontal lollipop represents a separate regional PRS, broken into three categories on the left (all individuals, females, and males). Dot size corresponds to the variance explained, R², of each PRS. The vertical dotted line represents the cutoff for significance (p_perm10,000<0.05).

After identifying the adipocyte MGSS WHRadjBMI PRS results as significantly enriched for variance explained for all individuals and females, we performed a partitioned heritability analysis to further corroborate these enrichments. We observed that the heritability (h²) of WHRadjBMI from the SNPs in the cis regions of the adipocyte MGSS set was significantly enriched in all individuals (h² enrichment=2.21, p=0.00150) and in females (h² enrichment=2.47, p=1.54× 10^-4), in line with the regional PRS analysis.

Taken together, the GWAS and PRS enrichments and significant heritability enrichment results highlight the adipocyte MGSS gene set as particularly important for the polygenic risk of abdominal obesity.

Identification of SREBF1 as the key SAT-enriched adipocyte marker gene

After identifying WHRadjBMI PRS enrichments in the adipocyte MGSS set, we explored the functional differences of the SAT and VAT gene sets. For each major cell-type, we tested for overrepresentation of genes in Gene Ontology (GO) categories using all three sets, MGSS, MGSV, and MGSBT (Supplementary Tables 4-10). There were no significant functional pathways for the macrophage MGSS and MGSBT sets. Among the adipocyte MGSV and MGSBT gene sets, we detected 170 and 9 significant pathways, respectively (FDR<0.05). Among the pathways of the adipocyte MGSV gene set, 103 (61%) and 98 (58%) pathways included the well-known insulin metabolism genes, INSR and IRS2, respectively. However, the strongest adipose tissue function-related pathway enrichments were observed in the adipocyte MGSS set (FDR<0.05), including regulation of fatty acid metabolic process, positive regulation of triglyceride metabolic process, and neutral lipid metabolic process amongst other breakdown-of-molecule centered processes (Supplementary Table 4).

In more detail, we found 23 significant functional pathways (FDR<0.05) for the adipocyte MGSS set (Figure 3a; Supplementary Table 4). Notably, we observed a strong presence of the master transcription factor (TF) of adipogenesis and fatty acid biosynthesis, SREBF1¹⁴, in these MGSS pathways. Of the 23 enriched functional pathways for adipocyte MGSSs (FDR< 0.05), 20 (87%) included SREBF1, making it the most prevalent pathway gene, and thus supporting its importance in the SAT adipocyte function. Consistent with the known role of SREBF1 as an adipose master TF^14,58,59, we also observed that 8 of the most prevalent pathway genes (genes present in 8 or more pathways) in the adipocyte MGSS gene set are predicted transcription targets of SREBF1 that have also previously been linked to obesity, including LEP⁶⁰ and FASN⁶¹ (Figure 3b).

• Download figure
• Open in new tab

Figure 3 Functional pathway enrichment analyses of cell- and tissue-type level expression profiles reveal distinctive differences in SAT and VAT pathways and highlight the adipose master transcription factor SREBF1 as the central pathway gene present in 87% of the adipocyte MGSS pathways.

(a) The top 10 most significant (FDR<0.05) functional pathways for adipocyte MGSS and MGSV are shown side by side, colored by adipocyte tissue-specific set. P-values for each pathway are shown on the x-axis, and each pathway is subsequently ranked by decreasing p-value. The point size of each dot corresponds to the enrichment ratio of the pathway. (b) Circos plots visualize the pathway memberships of the adipocyte MGSS genes that appear in >=8 adipocyte MGSS pathways and connect the genes to their respective pathways. Genes denoted with an asterisk are predicted target genes of SREBF1. Each gene is colored by the number of pathways it belongs to. We highlight SREBF1, an adipose tissue master TF and a key pathway gene that is present in 20 out of 23 pathways enriched in adipocyte MGSSs. (c) The cell-type level expression of SREBF1 across both fat depots is shown on the UMAP space, separately for SAT and VAT, and colored by normalized gene expression counts. d) The scaled expression counts over a 14-day adipocyte differentiation experiment of the longitudinally differentially expressed (DE) adipocyte MGSS genes with memberships in multiple enriched pathways are plotted and grouped by their longitudinal expression trajectory clusters.

Figure 3c shows how this key MGSS TF, SREBF1 exhibits a substantially higher adipocyte-enriched expression level in SAT than VAT. We confirmed the adipocyte unique marker gene status of SREBF1 in SAT snRNA-seq in the independent RYSA cohort (n=68) (average log₂ fold change=0.473, p_adj<2.23×10^-308), thus replicating its role as a significant SAT specific marker. Additionally, we replicated the lack of preferential expression of SREBF1 in VAT adipocytes, in an independent VAT snRNA-seq data from another cohort³⁴ comprising obese individuals (BMI>30) of European descent (n=5), which, similarly to the KOBS VAT samples, showed low, diffuse expression of SREBF1 in all VAT cell-types (Supplementary Figure 2).

Identification of temporally co-expressed adipocyte MGSS pathway genes during human adipogenesis

To link these functionally enriched MGSS genes to adipose tissue development and differentiation, we examined the temporal expression of the key pathway genes within the adipocyte MGSS set during human primary SAT preadipocyte differentiation (i.e., adipogenesis). We first identified 43 MGSS genes that contributed to more than one significantly (FDR<0.05) overrepresented pathway (Supplementary Table 4) and then found that the 42 of the 43 genes, expressed during SAT adipogenesis, were all longitudinally differentially expressed (DE) (FDR<0.05) during differentiation (Figure 3d). To link together genes with similar temporal transcription patterns, we clustered these longitudinally expressed genes using DPGP⁵¹ into 8 groups of temporally co-expressed genes (Figure 3d, Supplementary Table 11). Of note, we observed SREBF1 to cluster with two previously described obesity genes, LEP⁶⁰ and FASN⁶¹ (Figure 3d), further supporting the coordinated role of SREBF1 as a master TF with other key adipocyte marker genes.

Longitudinally differentially accessible peaks, co-accessible in the SREBF1 region during adipogenesis, harbor abdominal obesity GWAS variants

To elucidate the genetic regulatory role of SREBF1 in abdominal obesity, we first searched for open chromatin peaks with similar temporal accessibility patterns in the cis region of SREBF1 during human SAT adipogenesis (Figure 4a). We detected 129 peaks in the cis region (± 500 kb) of SREBF1 and found 120 (93%) to be longitudinally DA over the 14-day preadipocyte differentiation data (FDR<0.05), which we then clustered into 10 groups using DPGP⁵¹ (Supplementary Table 12). The clusters with peaks remaining after thresholding at a cluster assignment probability>0.9 are shown in Supplementary Figure 3. We focused on cluster 4, which contains 5 peaks with similar trajectories over the 14-day differentiation (Figure 4b).

• Download figure
• Open in new tab

Figure 4 Longitudinal trajectory analysis of open chromatin reveals co-accessible peaks in the cis region of SREBF1 that harbor WHRadjBMI variants.

(a) A regional schematic overview shows the proximity to SREBF1 of 5 co-accessible peaks that exhibit very similar open chromatin trajectory patterns. The peaks collectively harbor WHRadjGWAS variants, viewable through zoomed in loci of the chromosome. (b) The temporal accessibility patterns of the 5 peaks in the SREBF1 are plotted across the 14-day, 6 time-point SAT preadipocyte differentiation experiment. (c) The pairwise linkage disequilibrium (LD) (R²) in the UK Biobank (n=9,981) between the 7 WHRadjBMI GWAS SNPs residing within the co-accessible adipocyte peaks in the SREBF1 cis region is shown as a HaploView plot, colored by R².

To determine whether these five co-accessible peaks from cluster 4 in the SREBF1 region (Figure 4b) are involved in genetic risk of abdominal obesity, we next investigated them for overlaps with WHRadjBMI GWAS variants. Using the large, previously published GWAS summary statistics for WHRadjBMI from GIANT-UKB⁵², we identified 7 GWAS SNPs within these 5 peaks (Table 2). Of the 7 SNPs, rs4924755, rs7224725, and rs9944423 were significant (p<5×10^-8) in all individuals, driven by females (also significant in only females); rs35104205 was significant in all, females, and males; and rs12948060, rs4646347, and rs4646346 were significant only in males, reflecting the well-known sex effects of the WHRadjBMI trait (Table 2). Furthermore, rs4924755, rs7224725, rs9944423, and rs35104205 are also significant GWAS SNPs (p<5×10^-8) for the serum total triglycerides in all individuals from the European ancestry GLGC lipid GWAS⁶².

We observed tight LD patterns (R²>99%) among three of the 7 SNPs, rs12948060, rs4646347, and rs4646346, and separately among two SNPs rs7224725 and rs9944423, indicating 4 separate regional WHRadjBMI signals (Figure 4c). These results show that the 5 co-accessible peaks in the SREBF1 region contribute to the genetic risk of abdominal obesity.

Risk allele status for the abdominal obesity GWAS SNPs in the SREBF1 cis region affects SAT adipocyte expression of SREBF1 and more than hundred SAT adipocyte marker genes

We hypothesized that these seven abdominal obesity GWAS variants, which reside in the cis regulatory region of SREBF1, may influence the downstream expression of important adipocyte genes in trans. To address this hypothesis, we investigated whether the risk allele status of the seven WHRadjBMI GWAS SNPs in the SREBF1 cis region impacts SAT adipocyte expression of SREBF1 and downstream, the adipocyte expression of other adipocyte MGSSs, using independent SAT snRNA-seq data from another European bariatric surgery cohort, RYSA (n=68) (see Methods). We identified SREBF1 and 146 other adipocyte MGSS genes, for which the adipocyte expression is significantly affected by the risk allele carrier status of one or more of these seven WHRadjBMI GWAS variants in trans (see Methods) (Supplementary Table 13). Due to the tight LD (R²>99%) among rs12948060, rs4646347, and rs4646346, as well as among SNPs rs7224725 and rs9944423 (Figure 4c), their results are identical.

In more detail, we observed that the SAT adipocyte expression of SREBF1 was significantly impacted by each of the seven WHRadjBMI GWAS variants, being higher expressed in the risk allele carriers of rs12948060, rs35104205, rs4646347, and rs4646346, and lower expressed in the risk allele carriers of rs4924755, rs7224725, and rs9944423. The three SNPs in tight LD (rs12948060, rs4646347, and rs4646346) had the largest number of genes (n=115 genes) with significantly higher expression in the risk allele carriers (Supplementary Table 13), while the SNPs rs7224725 and rs9944423 in LD (R²>99) had the fewest (n=23 genes) (Supplementary Table 13). We found that the variant rs4924755 impacted the largest number of genes with lower expression in the risk-allele carriers (n=48 genes), including multiple functionally important adipose genes, e.g., FASN, AGPAT2, DGAT1, in line with the downregulation of SREBF1 in these carriers (Supplementary Table 13). The risk allele carrier status of this variant rs4924755 also associated with upregulation of numerous key adipose genes, such as PLIN1, PLIN4, and LIPE (Supplementary Table 13). We show the functional pathway enrichments of the gene sets affected by the risk allele carrier status in Supplementary Tables 14-16.

We also built module scores in the independent RYSA cohort using the SAT adipocyte expression of the up/down regulated trans genes by the risk allele status of the of the SREBF1 abdominal obesity GWAS SNPs (Figure 5a-b, Supplementary Figure 4). These average expression results of the trans genes further demonstrate the significant, wide-spread trans effects of the SREBF1 abdominal obesity GWAS SNPs on SAT adipocyte expression. Taken together, our results suggest that the WHRadjBMI GWAS SNPs in the cis region of SREBF1 affect SAT adipocyte expression of the major adipose tissue TF, SREBF1, and downstream of that, the adipocyte expression of 146 genes, comprising 70% of all adipocyte MGSS genes and including multiple key adipocyte MGSS genes.

• Download figure
• Open in new tab

Figure 5 Trans effects of SREBF1 abdominal obesity GWAS SNPs established using single cell level expression data and biobank investigation.

(a) Trans effects of abdominal obesity GWAS SNPs in the SREBF1 region on SAT adipocyte expression of the adipocyte MGSS genes in an allele-specific way is shown using module scores, constructed in the RYSA cohort (n=68) using the average SAT adipocyte expression of the upregulated trans genes by the risk allele status of the SREBF1 GWAS SNPs. For the module scores of the downregulated trans genes, see Supplementary Figure 4. The plots represent the module scores calculated for the upregulated trans genes by each SNP (or by a SNP representing a tight LD block, see Figure 4c). Thus, the module scores for the LD block of rs12948060, rs4646347, and rs4646346 and the LD block of rs7224725 and rs9944423 are shown using the data for rs12948060 and rs7224725. Each point represents a cell and its average gene expression of the trans genes of interest. (b) The differences between the carriers and non-carriers are also visualized through a two-sided Wilcoxon test and boxplot. (c) The trans gene sets are enriched for variance explained in abdominal obesity. Lollipop plots show the regional PRS results for each trans gene set (constructed in all individuals and females). Entries missing a lollipop did not either have enough SNPs to accurately build a PRS (<5 SNPs) or significantly predict the variance in WHRadjBMI (p_R²>0.05) after the regional PRS construction. (c) Differences in the magnitude of the regional PRSs between the carriers and non-carriers of the abdominal obesity risk allele for rs4924755 visualized through boxplots. Significance is calculated using a two-sided Wilcoxon test. Plots for the regional WHRadjBMI PRSs of the trans genes calculated in females and all individuals in UKB are displayed.

Adipocyte open chromatin variants in cis regions of the trans gene sets are enriched for variance explained in WHRadjBMI and their PRSs differ by the risk allele carrier status of the SREBF1 abdominal obesity GWAS SNPs

We hypothesized that the trans gene sets we identified had downstream effects on the polygenic risk of abdominal obesity by the risk allele carrier status of the SREBF1 abdominal obesity GWAS SNPs. Accordingly, we first examined the genetic risk contributions from the trans genes sets by constructing their regional PRSs for WHRadjBMI in UKB from the SNPs that fell in adipocyte open chromatin peaks in the cis regions of each trans gene set (see Methods). We found, through 10,000 permutations, that 12 out of the 21 PRSs which were built from variants from upregulated trans gene sets were significantly enriched (p_perm10,000<0.05) for incremental variance explained in abdominal obesity (Figure 5c, Supplementary Table 17), as were 6 out of the 10 PRSs built from variants from downregulated trans gene sets (Figure 5c, Supplementary Table 18), indicating that the regulatory regions of these trans regulated genes are important for abdominal obesity.

Next, we classified the UKB individuals as the carriers and non-carriers of the seven SREBF1 abdominal obesity GWAS risk alleles in the same way as in the RYSA adipocyte expression analysis and then compared the WHRadjBMI PRSs of the trans gene sets between the carriers and non-carriers in UKB to further assess the trans regulatory nature of the carrier status. We found that the SREBF1 cis regional variants exhibit significant differences in the magnitude of the regional WHRadjBMI PRSs of 10 out of the 21 upregulated trans gene sets between the carriers and non-carriers of the corresponding SREBF1 abdominal obesity GWAS SNP risk allele (Wilcoxon padj<0.05) (Supplementary Table 17), thus independently confirming the observed trans effects in RYSA adipocytes at the population level in UKB. Similarly, we found that 4 out of 10 downregulated trans gene sets also show significant PRS differences (Wilcoxon padj< <0.05) between the carriers and non-carriers (Supplementary Table 18). Figure 5d shows these significant differences between the regional WHRadjBMI PRSs of the trans genes by the risk allele carrier status of one of the 7 SREBF1 WHRadjBMI GWAS SNPs, rs4924755, and Supplementary Figure 5 illustrates similar significant WHRadjBMI PRS differences with several of the other SNPs. We also observed sex-specific differences in these PRS results (Figure 5c-d, Supplementary Tables 17-18), in line with the well-known sex differences in abdominal obesity. Overall, our results discover SREBF1 for SAT adipocyte function and genetic risk of abdominal obesity via variant-specific trans effects on numerous adipocyte trans genes.