Cell-Type Composition Affects Adipose Gene Expression Associations With Cardiometabolic Traits

Adipose tissue, including in the subcutaneous depot, affects development of cardiovascular and metabolic diseases, including obesity, type 2 diabetes (T2D), and dyslipidemia (1–3). Adipose gene expression associations with cardiometabolic traits can help interpret gene function. To understand potential mechanisms of genes on disease progression and development in adipose tissue, several studies have investigated the association between selected cardiometabolic traits and genes in 16–196 individuals (Supplementary Table 1). Larger sample sizes would increase power to detect associations and improve accuracy of effect size estimates.

A review of 25 adipose histology studies estimated composition as ∼60% adipocytes, 0.5–13% macrophages, and 0.5–2% endothelial cells (4). Cell types may have different gene expression profiles, and the variable proportions of cell types can influence energy homeostasis (3,5). This heterogeneity hinders determining whether gene expression differences between individuals are due to variation in gene activity within cell types, differing proportions of cells that express a gene, or both (3). Therefore, evaluating the transcriptomes of distinct cell populations can help elucidate how adipose heterogeneity influences metabolic conditions (3).

Although single-cell/single-nuclear RNA sequencing (scRNA-seq/snRNA-seq) technology is becoming more readily available, it is still costly (6), and the initial adipose snRNA-seq studies included relatively few samples (7–14). scRNA-seq cannot easily measure adipocytes due to their high lipid content; therefore, snRNA-seq is often used (6). snRNA-seq data can be used as a reference for deconvolution methods that estimate the cell-type composition of other tissue samples (7,13,15–18). Adjustment for cell-type composition can mitigate spurious associations between traits and gene expression (19). The cell-type proportion estimates can also help interpret trait-gene expression associations and potential gene mechanisms.

Associations between adipose gene expression and/or cell-type proportion with BMI have been described (7,13,20,21). BMI is associated with several cardiometabolic traits (22) and with gene expression levels (20,23), although the directions of causality are uncertain. An association between macrophage proportion and BMI has been reported, but the direction of effect between adipocyte proportion and BMI was inconsistent in three studies (7,13,20). Further, little is known about how cell-type composition and BMI affect associations between gene expression and cardiometabolic traits.

Here, we used snRNA-seq data from 16 adipose samples as a reference for cell-type deconvolution of 434 previously described and 425 new, a total of 859, subcutaneous adipose samples with bulk RNA-seq from the METabolic Syndrome in Men (METSIM) study. We tested the association between cell-type proportion and 22 cardiometabolic traits, including BMI, and determined that 17 traits were associated with a proportion of at least one cell type. We then evaluated inclusion of BMI and/or cell-type proportion as covariates in the best-fit models for a majority of trait-gene expression associations genome-wide. Using optimal linear models, we identified 2,664 significant associations across 1,252 genes and 14 traits.

The METSIM study is a population-based cohort of 10,197 men from Kuopio, Finland (24). All participants gave written informed consent. The study was approved by the Ethics Committee of the University of Eastern Finland and Kuopio University Hospital in Kuopio, Finland, and conducted in accordance with the Declaration of Helsinki. We studied 434 participants with subcutaneous adipose tissue collected near the umbilicus by needle biopsy, previously described (25), and 425 additional participants for which a physician obtained a surgical biopsy of ∼0.5–1 g of subcutaneous adipose tissue near the umbilicus. These samples were collected for research purposes and not secondary to other surgery. We evaluated 22 cardiometabolic traits measured closest to the time of the biopsy (Supplementary Table 2).

RNA-seq of needle biopsy samples to an average depth of 45 million paired-end 50-base pair (bp) reads has been described; the percentage of uniquely mapped reads averaged 82.2% (25).

We performed RNA extraction of frozen surgical biopsy specimens at The University of North Carolina Biospecimen Processing Core. After homogenizing in TRIzol and including glycogen, we quantified RNA concentration by Qubit and quality by Agilent TapeStation. For 504 samples with an RNA integrity number (RIN) ≥6, we performed RNA-seq at The University of North Carolina High Throughput Sequencing Core. We performed polyA+ selection and library preparation (KAPA Stranded mRNA), prepared dual-barcoded individual libraries, and sequenced to an average depth of 82.6 million paired-end 150-bp reads (Illumina NovaSeq 6000 S4Xp). We used the Illumina bcd2fastq pipeline (v.2.20.0) and dropped four samples due to poor sequencing quality. We assessed RNA-seq quality by applying FastQC (v.0.11.8) and MultiQC (v.1.7) (26,27). We filtered adaptor sequences and sequencing reads using Cutadapt (v.1.18), requiring a Phred quality score >20 and a minimum length of 50 bp (28). We used STAR (v.2.7.3a) two-pass alignment to align the reads to the hg19/GRCh37 reference, using the GENCODE v.19 annotation (29,30) and ENCODE alignment parameters (30). We retained duplicate reads and removed unpaired alignments. The uniquely mapped reads across samples averaged 74.3%. We used QTLtools mbv (v.1.1) to compare RNA-seq and genotype data concordance and retained 492 samples (31,32).

For subsequent analyses, we excluded surgical biopsy data for individuals who also contributed needle biopsy specimens or contributed two surgical samples. For the resulting needle (n = 434) and surgical biopsy specimens (n = 425), we quantified reads using QTLtools quan (v.1.1), excluding reads with >5 mismatches and requiring mapping quality >150 (29,31), and retained genes with ≥5 counts in ≥25% of samples. We normalized for library size by adjusting counts for trimmed mean of M (TMM) (33) using edgeR (34).

We used snRNA-seq raw counts and cluster information of 16 previously described subcutaneous adipose tissue samples from the Finnish Twin and CRYO studies (7 men, 9 women; average BMI 31.2) (35); these data derive from the same population and adipose depot as the bulk RNA-seq samples. Briefly, nuclei were isolated from frozen biopsy specimens, processed by 10X Chromium, and sequenced to ∼50,000 reads per nucleus. Reads were aligned to GRCh36 using the GENCODE v.19 gene annotation (29). After clustering, droplets were removed for contamination, <200 or >20,000 unique molecular identifiers, <200 genes, or >20% reads mapping to the mitochondrial genome. The resulting counts were log-normalized, the first 30 principal components (PCs) were used for clustering, and cell types were assigned using reference data sets. The data set includes 16 cell-type clusters: adipocytes, B cells, lymphatic endothelial cells, endothelial cells, fibroblasts, macrophages (M1–M5), mast cells, neutrophils, perivascular cells, preadipocytes, and T cells (natural killer [NK] and naïve) (Supplementary Table 3) (35). We converted Ensembl IDs to gene symbols to match the bulk and snRNA-seq formats; among duplicate gene symbols, we kept the one with the highest average gene expression. We removed snRNA-seq cell types that contained <1% of nuclei, which excluded lymphatic endothelial cells, retained the 11,730 genes with average gene expression ≥0.01 counts, and combined the five macrophage subtypes. We converted raw counts to expression sets in R (v.4.1.0) using Biobase (v.2.50.0) (36).

We estimated the cell-type proportions of the needle and surgical biopsy data sets separately using Multi-subject Single-cell (MuSiC) deconvolution and SCDC (v.1.0) (16,18). We used the music_prop() and scdc_prop() functions that leverage cell-type gene expression from a reference to estimate cell-type proportions. Additionally, we estimated cell-type proportions using the SCDC tree-based clustering method, for which we performed hierarchical clustering using Euclidean distance on the snRNA-seq data set. The resulting dendrogram (Supplementary Fig. 3) showed five clusters: adipocytes; fibroblasts, preadipocytes, endothelial cells, and mast cells; T-naïve cells, T-NK cells, and macrophages; neutrophils and B cells; and perivascular cells. We used these clusters and the SCDC_prop_subcl_marker() function to select markers and estimate cell-type proportions. For the remaining analyses, we used the SCDC tree-based clustering. We compared the proportions between the needle and surgical biopsy specimens using the prop.test() function in R.

We determined the first 10 gene expression PCs using QTLtools pca (v.1.1.0) on the gene expression counts of the needle and surgical biopsy specimens separately. We calculated Pearson correlations (r) between PCs, known technical covariates, and cell-type proportions using the cor() function in R. Technical covariates for both data sets included read deletion size, mean read insertion size, age, RIN, and sequencing batch. We inverse normal transformed 22 cardiometabolic traits and converted the cell-type proportions to z-scores using the scaled() function in R. We used the lm() function to test association between each cell-type proportion and trait. We confirmed associations using Dirichlet regression (37). We summed cell types with perfectly correlated estimated proportions (Pearson r = 1) to make five groups corresponding to the dendrogram (Supplementary Fig. 3), used the DirichReg() function, and modeled composition as the dependent variable. We plotted results using ggplot2 (38).

We inverse normal transformed gene expression and trait measurements and tested for association between 25,609 genes (needle) or 29,753 genes (surgical) with 21 traits (excluding BMI) using linear regression from the lm() function of R. We adjusted for the covariates listed above, except for age when testing for association with age. We performed the linear regression with four models adjusting for different covariates:

• 1. Trait ∼ gene expression + known covariates

• 2. Trait ∼ gene expression + known covariates + BMI

• 3. Trait ∼ gene expression + known covariates + Macro-phage +Endothelial + Perivascular + B cells

• 4. Trait ∼ gene expression + known covariates + BMI + Macrophage + Endothelial + Perivascular + B cells

To determine the best fit for each trait-gene expression association among these models, we used the AIC() function from the R stats package (39).

We performed a random-effects meta-analysis on each trait-gene association using the rma() function in the metaphor package (v.2.4.0) in R. We derived a Bonferroni significance threshold (P ≤ 2e−6) based on genes tested in both data sets (0.05/25,162).

We used our prior adipose analyses of gene expression quantitative trait loci (eQTL) and splice junction usage (sQTL) that colocalized with genome-wide association study (GWAS) signals (25,40). We extracted the e/sGene-GWAS pairs for traits measured in METSIM participants and used the Matsuda index and HOMA of β-cell function (HOMA-B) to correspond with the T2D GWAS because they were the most similar measures available. We extracted trait-gene expression associations (meta-analysis P ≤ 1e−3) with a corresponding e/sGene-GWAS pair. We created LocusZoom plots (41).

The full trait-gene expression associations are available at https://mohlke.web.unc.edu/data/. METSIM data are available in dbGaP phs000743.v3. The snRNA-seq data counts are available on Gene Expression Omnibus GSE236708.

We collected subcutaneous adipose tissue biopsy specimens and performed RNA-seq on 859 individuals from the METSIM study who have extensive cardiometabolic trait measurements (21,24,25). We evaluated specimens from 434 needle biopsies and 425 surgical biopsies (Supplementary Tables 2 and 4 and Supplementary Fig. 1). Owing to differences in biopsy type, trait values, and sequencing methods, we analyzed the biopsy groups separately and later meta-analyzed.

To determine the cell-type composition of adipose samples with bulk RNA-seq data, we performed cell type deconvolution using an independent snRNA-seq data reference obtained from the same adipose depot and geographical population (7,11). Of 12 cell types in the snRNA-seq data set, the most abundant were macrophages (25.7%), adipocytes (24.6%), and fibroblasts (13.7%) (Supplementary Tables 5 and 6 and Supplementary Fig. 2); the low percentage of adipocytes may reflect limitations of snRNA-seq (6) and is within the range of 16–59% adipocytes reported in other snRNA-seq studies (11,13). We compared three deconvolution strategies, MuSiC, SCDC, and SCDC tree-based. The SCDC tree-based method, which included clustering snRNA-seq cell types (Supplementary Fig. 3), estimated the highest median proportions of adipocytes, 0.79 and 0.71 for needle and surgical biopsy specimens, respectively (Fig. 1, Table 1, and Supplementary Figs. 4 and 5). The default SCDC and MuSiC methods yielded similar results to each other, and both unexpectedly estimated endothelial cells to be most abundant (Fig. 1). Based on similarity to cell-type proportions in a literature review (60% adipocytes, 0.5–13% macrophages, and 0.5–2% endothelial cells) (4) and other adipose deconvolution estimates ranging from 50 to 85% adipocytes (4,7,20), we used SCDC tree-based deconvolution for further analyses.

We compared the estimated cell-type proportions by biopsy type. Adipocyte proportions ranged from 0.36 to 1.00 for needle biopsy specimens and from 0.32 to 0.89 for surgical biopsy specimens (Table 1 and Fig. 1). Overall, the surgical biopsy specimens showed less variation across samples (Fig. 1 and Supplementary Figs. 5 and 6). The median cell-type proportions did not differ by biopsy type; however, the maximum proportions for adipocytes, neutrophils, and B cells were significantly (P < 0.05) higher in the needle biopsy specimens (Table 1). Among the needle biopsy specimens, 67 samples with low (<0.65) adipocyte proportions had higher neutrophil proportions (median 0.06 vs. median 0.0 in 367 remaining samples; P = 0.04), which could reflect higher blood content, as observed previously (25) (Supplementary Fig. 7). In contrast, surgical biopsy specimens had a median neutrophil proportion of 0. Thus, the adipose surgical biopsy specimens have more consistent cell-type composition than needle biopsy specimens, potentially due to mechanical stress during needle sampling, which can cause tissue breakdown and uneven cell composition.

To assess the roles of cell-type proportion and technical factors in gene expression variation across individuals, we measured the correlation between cell-type proportions, known technical covariates, and gene expression PCs. Cell-type proportion estimates were correlated with technical factors, sequencing batch, BMI, and gene expression PCs (Supplementary Fig. 8). These correlations suggest that cell-type proportions could at least partly explain some gene expression PCs, consistent with an effect of variation in cell-type proportion on gene expression.

We evaluated the effect of cell-type proportion estimates on cardiometabolic traits, including BMI, using linear regression. For both biopsy types, we observed the strongest effects for the Matsuda index, BMI, and waist circumference (WC) associated with adipocytes, macrophages, and T cells (Fig. 2, Supplementary Fig. 9, and Supplementary Table 7). HDL-cholesterol, fat-free mass percentage, adiponectin levels, and the Matsuda index were associated with higher adipocyte proportion, while the remaining traits with significant associations were associated with lower adipocyte proportion (Fig. 2, Supplementary Fig. 9, and Supplementary Table 7). While most associations and directions of effect between cell-type proportions and traits were shared between biopsy types, age, HbA_1c, total cholesterol, and LDL-cholesterol were only significantly associated with macrophages and T cells in needle biopsy specimens (Fig. 2, Supplementary Fig. 9, and Supplementary Table 7).

Given that BMI has shown association with expression levels of most genes (25), we evaluated the relationship between BMI and cell-type proportion. BMI was positively associated with macrophage proportion (needle r = 0.35, surgical r = 0.45), consistent with an increase in macrophage infiltration with higher BMI (Fig. 2, Supplementary Table 7, and Supplementary Figs. 8–10). In contrast, BMI was negatively associated with adipocyte proportion (needle r = −0.25, surgical r = −0.41) (Fig. 2, Supplementary Table 7, and Supplementary Figs. 8–10). Models that account for the compositional nature of the cell-type proportions showed largely similar results (Supplementary Fig. 11 and Supplementary Table 8). Overall, cell-type proportion was associated with 17 of 22 cardiometabolic traits, suggesting that specific gene expression associations with these traits may differ based on cell-type composition.

Initially, we measured the association between 21 cardiometabolic traits (excluding BMI) and adipose gene expression on a genome-wide scale without adjusting for BMI or cell-type proportion. We found 4,752 of 25,509 genes (18.6%) tested in the needle biopsy specimens had a significant association (P < 2e−6) with at least one trait, as did 6,558 of 29,753 genes (22.0%) tested in the surgical biopsy specimens (Supplementary Table 9). Overall, we observed 25,145 and 35,213 trait-gene expression associations in the needle and surgical biopsy specimens, respectively.

Next, we examined models incorporating BMI and adipose cell-type proportions as covariates. We identified trait-gene expression associations with 21 traits using four linear models (models 1–4, Research Design and Methods). For needle biopsy specimens, we observed 25,145 significant (P < 2e−6) associations with model 1 and 472 with model 4 (Supplementary Table 9); 25,942 trait-gene associations were significant in at least one model. For surgical biopsy specimens, we observed 35,213 associations with model 1 and 1,094 with model 4 (Supplementary Table 9); 36,343 trait-gene associations were significant in at least one model. For all associations significant in both models, the effect directions were consistent. We used Akaike Information Criteria (AIC) values to select the model with the best fit for each trait. The best-fitting model nearly always included BMI: 25,483 of 25,942 trait-gene associations (98.2%) in needle biopsy specimens and 36,006 of 36,343 trait-gene associations (99.1%) in surgical biopsy specimens (Fig. 3). Only for adiponectin and age did best-fitting models not include BMI for a majority of genes, and glomerular filtration rate (GFR) showed no significant associations in either data set (Supplementary Table 10). Furthermore, model 4, which adjusted for both BMI and cell-type proportions, was the best fitting model for 20,536 (79.2%) and 29,540 (81.3%) of all trait-gene associations for needle and surgical biopsy specimens, respectively (Fig. 3 and Supplementary Table 10). Although excluding BMI and cell-type proportions from models resulted in more ostensibly significant findings, AIC values indicated that including these covariates provided better models. Thus, both BMI and cell-type proportion affected associations between adipose gene expression levels and cardiometabolic traits.

We meta-analyzed the trait-gene expression associations from the needle and surgical biopsy data sets using the best fitting model adjusting for BMI and cell-type proportions for all traits, except age and adiponectin that we adjusted only for cell-type proportions (Supplementary Table 11). We identified 2,664 significant associations across 1,252 genes and 14 traits (Table 2 and Supplementary Table 11). Of the 25,162 genes tested, 1,252 have a significant association with at least one trait (Table 2). Overall, the meta-analyses showed at least 2.4 times more significant associations compared with either of the individual studies.

The traits with the most significant associations were WC, Matsuda index, adiponectin, fasting insulin, HOMA-B, and waist-to-hip ratio (WHR) (Table 2, Fig. 4, and Supplementary Table 11). For each of these traits, we ranked the genes by absolute effect sizes; of the 10 genes with the strongest effects, 15 were shared among two or more traits: ITIH5, SPX, ALPK3, AGPAT9, HADH, GPD1L, UCHL1, AZGP1, GPR146, NIPSNAP3B, NAALAD2, MSC, SLC24A, STBD1, and FAM47E (Fig. 4 and Supplementary Table 6). Higher expression of SPX was associated with a higher Matsuda index and lower fasting insulin and HOMA-B levels (Fig. 4 and Supplementary Table 6). SPX encodes a hormone involved in cardiovascular and renal function, and one study showed that SPX levels are lower in obese children and negatively correlated with insulin resistance (42,43). In addition, higher AGPAT9 (GPAT3) was strongly associated with a higher Matsuda index and lower WC and WHR (Fig. 4 and Supplementary Table 6). GPAT3-depleted mice showed less body weight gain and adiposity and increased energy expenditure on a high-fat diet (44). These results highlight genes for which differential adipose expression is associated with several cardiometabolic outcomes.

We compared trait-associated genes identified here with genes previously implicated by adipose eQTL or sQTL colocalized with cardiometabolic trait GWAS signals (25,40). Among 235 colocalizations with adipose eQTL (25), 25 (10.6%) also had a nominally significant trait-gene expression association, at an arbitrary P-value threshold of P ≤ 1e−3, despite a sample size of just 859 participants (Supplementary Table 12); 80% of these had consistent effect directions. The traits with the most trait-gene associations corresponding to colocalized GWAS signals include WHR and T2D. Of the 46 eQTL signals colocalized with a GWAS signal for WHR, 6 (13%) also showed association between gene expression level and WHR (Supplementary Table 12 and Supplementary Fig. 12). Of the 36 eQTL signals colocalized with a GWAS signal for T2D, 4 (11%) also showed association between gene expression level and the Matsuda index (Supplementary Table 12 and Supplementary Fig. 12). Among 110 genes implicated by sQTL-GWAS colocalization, 3 showed nominally significant trait-gene expression association (Supplementary Table 13). PABPC4 is an example of an eQTL gene colocalized with a T2D GWAS signal that also showed association between PABPC4 expression and the Matsuda index (Supplementary Tables 6 and 12 and Fig. 5). These genome-wide adipose gene expression level associations with cardiometabolic traits provide further support for the validity of GWAS-eQTL colocalization results.

We estimated the cell-type composition of 859 adipose tissue samples with bulk RNA-seq, of which 425 are first reported here. We observed considerable variation in cell-type proportions between individuals, including that needle biopsy specimens showed more variation in cell-type proportions than surgical biopsy specimens. Cell-type proportion estimates were correlated with gene expression, BMI, and 17 of 21 other cardiometabolic traits (Fig. 2). We identified optimized models of association between traits and gene expression and observed that even after adjusting for BMI and technical covariates, the optimal models also included cell-type proportions as covariates for 19 traits. Using models that accounted for cell-type proportions, we observed 2,664 significant associations between adipose expression level of 1,252 genes and 14 cardiometabolic traits in 859 individuals.

The directions of association for many pairs of traits and cell types reflected known mechanisms, although a wider range of trait measurements, especially age, could provide more accurate correlations. For example, we observed a positive association between WHR and adipocyte proportion, consistent with findings that a higher adipocyte proportion corresponds to more fat around the stomach than the hips (20). We also observed a positive association of BMI with macrophage proportion, consistent with more macrophage infiltration at higher BMI (45,46), and a negative association with adipocyte proportion, consistent with a larger volume of adipocytes in individuals with higher BMI (47). These results are consistent with a previous study of 100 lean and obese individuals with adipose tissue deconvoluted using Bisque and snRNA-seq of 6 samples (7) and a study of 331 METSIM needle biopsy specimens deconvoluted using MuSiC and snRNA-seq of 13 samples (13), but not consistent with a study of 652 women with adipose tissue deconvoluted using CIBERSORT and four independent cell type references (20). The inconsistency could be due to different deconvolution methods or the choice of references. Understanding the effects of cell-type proportions on cardiometabolic traits can reveal potential unknown biological mechanisms and elucidate cell-type contributions in disease.

We sought linear models that were broadly applicable to trait-gene expression associations with 21 traits genome-wide. Using AIC values to determine the optimal linear model for the largest number of genes and traits, we found that including BMI as a covariate improved model fit for >95% of genes in 18 of 21 traits. We considered stepwise linear regression with each cell-type proportion to determine the optimal model for every trait-gene pair tested; however, with cell-type estimates that are relative and, in some cases, perfectly correlated, interpreting which cell types contribute the most to particular traits would be challenging. Instead, identifying broadly optimal models for most trait-gene associations allowed us to mitigate the effects of variation in cell-type proportions across individuals. In general, our findings suggest that many trait-gene expression associations may be influenced by cell-type composition and should be interpreted carefully when adjustment for cell-type proportion is not available.

We used the data-rich METSIM study, which has many individuals with both gene expression and cardiometabolic traits measured. Trait-gene expression associations provided support for GWAS-QTL colocalized signal pairs and potential directions of effect. While mediation methods could provide evidence of causal effects of variants on genes to affect traits, the hundreds of samples here have low power compared with the thousands of samples used to discover the GWAS associations. This discrepancy between sample sizes could partially account for the limited overlap between GWAS-eQTL colocalizations and trait-gene associations. Nevertheless, 14 of the 25 trait-gene associations with a corresponding GWAS-eQTL colocalization also previously showed evidence of mediation (25).

Our strategy to include cell-type proportions when identifying trait-gene expression associations has limitations. Although snRNA-seq can assay adipocytes in adipose better than droplet-based methods of scRNA-seq (48), nuclear gene expression levels may differ from cytoplasm levels. The 24% adipocytes in the snRNA-seq data set may still be an underestimate for this reason, and the reference from only 16 individuals may not fully represent gene expression in the 859 individuals. Also, cell type estimates differed greatly by deconvolution method, although the SCDC tree-based had adipocyte proportions similar to previous estimates (4,7,20). Each deconvolution method produced relative, not absolute, cell-type proportions, and clustered cell types were perfectly correlated; using these proportions as covariates affected decisions to include a cell type in the best fitting model. Additionally, these cell-type proportions were estimated in men and may differ in women. These results highlight the need for further development of deconvolution methods or larger snRNA-seq studies.

In summary, we report associations between traits and adipose gene expression in 859 individuals adjusted for cell-type contributions. The emerging studies that have measured gene expression by cell type provide useful information about the proportions of each gene’s expression detected in each cell type. Ultimately, gene expression measured by cell type in larger studies of clinically characterized individuals is needed to interpret which genes act in specific cell types to influence cardiometabolic traits.

Acknowledgments. The authors thank the METSIM study investigators and participants for providing the subcutaneous adipose tissue samples and clinical trait measures that made this study possible. The authors thank Jiawen Chen (University of North Carolina Department of Biostatistics) and Yun Li (University of North Carolina Departments of Biostatistics and Genetics) for advice on deconvolution methods and models.

Funding. This study was supported by National Institutes of Health grants F31HL154730, T32GM007092 (to S.M.B.), R01DK079357, R01DK072193, UM1DK126185, R01DK132775 (to K.L.M.), R01DK062370 (to M.B.), T32HL129982 (to J.D.R.), and R01HG010505 (to P.P); Academy of Finland (272376, 266286, 314383, 335443), Finnish Medical Foundation, Finnish Diabetes Research Foundation, Novo Nordisk Foundation (NNF10OC1013354, NNF17OC0027232, NNF20OC0060547), Gyllenberg Foundation, Sigrid Juselius Foundation, Government Research Funds for Helsinki University Hospital and the University of Helsinki to K.H.P.; Academy of Finland (338417), Finnish Diabetes Research Foundation, Orion Foundation, Paavo Nurmi Foundation, Helsinki University Hospital Research Funds, Paulo Foundation, Finnish Medical Foundation, and Ida Montin Foundation to S.H.; and the Academy of Finland (321428), Sigrid Juselius Foundation, Finnish Foundation for Cardiovascular Research, and Centre of Excellence of Cardiovascular and Metabolic Diseases (0245896-3) supported by the Academy of Finland to M.L.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. S.M.B. drafted the manuscript. S.M.B., A.O., H.J.P., S.V., C.P., S.D., P.V.B., and H.M.S. generated data. S.M.B., J.D.R., M.A., S.H., J.M.V., and M.N.N. analyzed data. S.M.B., J.D.R., M.I.L., and K.L.M. interpreted the results. S.M.B. and K.L.M. designed the study. A.O., M.A., S.H., B.W.v.d.K., L.F.S., H.M.S., M.B., J.K., K.H.P., P.P., and M.L. provided resources. All authors reviewed and edited the manuscript. K.L.M. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.