Detailed Functional Characterization of a Waist-Hip Ratio Locus in 7p15.2 Defines an Enhancer Controlling Adipocyte Differentiation.

We combined CAGE sequencing in human adipocytes during differentiation with data from genome-wide association studies to identify an enhancer in the SNX10 locus on chromosome 7, presumably involved in body fat distribution. Using reporter assays and CRISPR-Cas9 gene editing in human cell lines, we characterized the role of the enhancer in adipogenesis. The enhancer was active during adipogenesis and responded strongly to insulin and isoprenaline. The allele associated with increased waist-hip ratio in human genetic studies was associated with higher enhancer activity. Mutations of the enhancer resulted in less adipocyte differentiation. RNA sequencing of cells with disrupted enhancer showed reduced expression of established adipocyte markers, such as ADIPOQ and LPL, and identified CHI3L1 on chromosome 1 as a potential gene involved in adipocyte differentiation. In conclusion, we identified and characterized an enhancer in the SNX10 locus and outlined its plausible mechanisms of action and downstream targets.

Introduction

According to estimations by the World Health Organization, more than 1.9 billion adults were overweight and ∼13% of the world's adult population was obese in 2016. In addition, 41 million children younger than 5 years were overweight or obese. In parallel with this, dramatic increases of downstream consequences of obesity, such as insulin resistance (IR), type 2 diabetes (T2D), and cardiovascular disease are also expected over the next decade, which will have a devastating impact on global health. Once considered as a problem limited to high-income countries, obesity, IR, and associated conditions are now dramatically increasing also in low- and middle-income countries, especially in urban settings. Heart disease and stroke are the global leading causes of morbidity and mortality, being the underlying cause of death for ∼13 million people in 2010, or one in four deaths worldwide (Lozano et al., 2012). Hence, identifying molecular mechanisms underlying obesity and risk of cardiovascular disease and T2D are important public health priorities.

Genome-wide association studies (GWAS) have been tremendously successful in identifying loci associated with complex traits (Visscher et al., 2017). As a result of large meta-analyses of GWAS led by the GIANT consortium, more than 500 loci have been associated with body mass index (BMI), waist-hip ratio (WHR), and other obesity traits (Locke et al., 2015, Lu et al., 2016, Shungin et al., 2015, Turcot et al., 2018, Winkler et al., 2016, Pulit et al., 2019). However, the causal variants, genes, or downstream mechanisms have been established for very few of these loci; hence, the transformative potential of human genomics to unravel the molecular mechanisms underlying obesity, IR, and adipocyte biology remains to the realized.

In the present study, we characterized an active enhancer located in one of the first chromosomal regions to be associated with WHR (Heid et al., 2010), most strongly in women (Shungin et al., 2015, Winkler et al., 2015), and also associated with BMI (Monda et al., 2013), endometriosis (Nyholt et al., 2012, Sapkota et al., 2017), and kidney function (Pattaro et al., 2016) in several GWAS (Figure 1A). The GWAS associations have mostly been annotated as SNX10 based on expression quantitative trait loci (eQTL) analyses (Graff et al., 2017, Justice et al., 2017, Shungin et al., 2015), and also as NFE2L3, MIR148A, or CBX3 based on the gene nearest to the lead GWAS variant from the study in question. The enhancer is localized in an independent haplotype block tagged by the lead GWAS variant rs3902751 that spans 48 kb from rs10245353 to rs1451385 (Figure 1A).

Figure 1

Genetic Features of the 7p15.2 Region Tagged by GWAS Potentially Involved in Adipogenesis Centered around rs1451385

(A) Physical map of the region showing the position of annotated genes, MIR148A, and genetic polymorphisms (SNPs, annotated with their rs identifiers) associated with obesity and other related phenotypes. The region is delimited by the distal gene NPVF and the proximal gene SNX10. Genome-wide significant SNPs include lead variants associated with WHR (rs10245353, rs3902751, rs1534696, Shungin et al., 2015; rs1055144, rs7798002, rs1451385, Heid et al., 2010) and BMI (rs10261878; Monda et al., 2013). rs1451385 is located in a haplotype block of 48 kb (LD r2 > 0.8) that is delimited by the rs10245353 and rs1451385. rs10261878 (associated with BMI) is not included in this haplotype block; instead it is linked to a haplotype block that includes MIR148A. The SNP rs1534696 at SNX10 is located in an independent haplotype. SNPs indicated by # are located at transcriptional peaks as indicated by CAGE. The dotted line represents 17 SNPs with a weak association with BMI (-log p values >3.5; Speliotes et al., 2010).

(B) Transcriptional signals during adipogenesis of mesenchymal stem cells visualized using the ZENBU browser (http://fantom.gsc.riken.jp/zenbu/). The variant rs1451385 is located in a transcriptionally active region of 400 bp delimited by CAGE25894547 and CAGE25894945. The directions of the transcription are indicated with green and purple arrows for forward and reverse directions, respectively. The intensities of the green and purple peaks correspond to the amount of CAGE-RNAs captured during the entire induction of adipogenesis.

(C) CAGE transcriptional signals colored by intensity on the black background and ordered from top to bottom from 1 h after induction of adipogenesis to day 14. The signals are aligned with the corresponding peaks in (B). A transcriptional promoter-like and an enhancer RNA (eRNA) region are predicted by FANTOM5 CAGE data and indicated with white filled and dotted arrows, respectively. The transcriptional like-promoter sequence is annotated as a POL2-binding site by ENCODE chromatin immunoprecipitation sequencing data, and the enhancer region is annotated as an enhancer according to the Regulatory Elements from ORegAnno at the UCSC Genome Browser in various cell lines.

(D) Numerical quantification of the CAGE signals as in (C) showing the three biological replicates for each time point. Maximal transcription is scored in both directions at 12 h post-induction.

(E) Chromosomal position at the short arm of chromosome 7 and conservation of the region flanking rs1451385 (marked by a red line) depicted by the vertebrate Multiz Alignment at the UCSC Genome Browser. Based on this phylogenetic conservation in mammals, an arbitrary sequence of 266 bp was chosen for further molecular characterization. The CAGE transcriptional signals are displayed as in (B).

In this study, we identified and characterized this active enhancer. Disrupting it causes the loss of differentiation capacity from precursor cells to mature adipocytes. Furthermore, combining CRISPR-Cas9 gene editing with global RNA sequencing (RNA-seq), we identified and characterized potential downstream targets of the enhancer. We aimed to disentangle the molecular mechanisms behind one of the first and strongest GWAS signals associated with obesity-related traits.

Results

In Silico Analyses Highlight a Regulatory Region Tagged by GWAS Potentially Involved in Adipocyte Biology

Using the GWAS catalog (www.ebi.ac.uk/gwas/), we identified all lead variants associated with obesity-related traits (WHR, BMI, lipid levels, glucose-related traits) at the time of study initiation (October 2015). Next, we selected all common single nucleotide polymorphisms (SNPs) in high linkage disequilibrium (LD; r2≥0.8) with the lead variants and searched for overlaps with annotated enhancers using FANTOM5 CAGE data from adipogenesis (Arner et al., 2015, Ehrlund et al., 2017) favoring (1) enhancers with high gene expression during induction of adipogenesis from mesenchymal stem cells and (2) enhancers showing higher expression levels in adipocytes and preadipocytes than in other cell types. This approach identified 14 putative adipocyte-specific enhancers with supporting evidence from GWAS of obesity-related traits. Of these, nine had weak human genetics support (non-genome-wide significant association and/or reported in a small GWAS without further replication), two were associated with proinsulin and low-density lipoprotein cholesterol, respectively (making them less attractive for studies in adipocytes), and two were upstream of NEGR1 (which is an already well-studied gene). In contrast, the enhancer in 7p15.2 had strong evidence for involvement in adipogenesis from FANTOM5, represented a very strong GWAS signal associated with WHR adjusted for BMI, and had supporting evidence from transcription-binding motifs and chromatin immunoprecipitation sequencing from ENCODE, but little was known about its function.

In this GWAS locus, the FANTOM5 enhancer data highlighted a bidirectional transcription start site active during adipocyte differentiation flanking rs1451385, suggesting that this SNP is located in a functional element of relevance for adipogenesis, as well as fat distribution (Figures 1B–1D). Furthermore, data from ENCODE and the Roadmap Epigenomics Project show that rs1451385 colocalizes with a DNaseI hypersensitive site and with dynamic changes of chromatin acetylation of histone 3 lysine 27 (H3K27ac), suggesting an active regulatory element in this region. The sequence comprising rs1451385 co-immunoprecipitates with the following transcription factors: CTCF (CCCTC-binding factor), USF1 (upstream transcription factor 1), IRF1 (interferon regulatory factor 1), and POLR2A (RNA polymerase II subunit A). USF1 is a transcription factor controlling expression of several genes involved in lipid and glucose homeostasis (Putt et al., 2004) that has been linked to familial combined hyperlipidemia (Pajukanta et al., 2004). USF1 binds to E-box motifs (5′-CACGTG-3′), and it is positioned about 23 nucleotides from a pyrimidine-rich region (Massari and Murre, 2000). The rs1451385 polymorphism is located 2 bp upstream of the E-box motif consensus sequence (5′-ACACGTGA-3′), which is located 23 nucleotides upstream of a 31-nucleotide-long pyrimidine rich-region (23 pyrimidines of a total of 31 nucleotides). Furthermore, analyses of nucleotide-binding sequences using PROMO (Messeguer et al., 2002) show differences in putative binding sites for several transcription factors. Specifically, the T allele of rs1451385, which is in perfect LD with the risk allele (A) at rs3902751 associated with increased WHR (Shungin et al., 2015), confers an extra VDR (vitamin D receptor)-binding site, a FOXP3 site, and a PXR-1:RXR-alpha (a nuclear receptor involved in metabolism sensing)-binding site. Finally, rs1451385 is a borderline significant eQTL (p = 6.7 × 10−6) in subcutaneous adipose tissue in GTEx (www.gtexportal.org/home) for AC003090.1 (ENSG00000223561.2), a long intergenic noncoding RNA (lincRNA) with unknown function, located 100 kb from rs1451385. However, due to its low expression in adipose tissue, it is unlikely that this lincRNA plays an important role in adipogenesis. There were no other eQTLs for subcutaneous or visceral fat in this locus.

The Sequence Comprising the GWAS Signal Shows High Enhancer Activity

We used a luciferase reporter assay to test the regulatory function of the region surrounding rs1451385. Based on phylogenetic conservation (Figure 1E), we cloned a 266-bp fragment (from now on called Enh#385) into the reporter vector pGL4.10 (Promega) with a minimal promoter upstream of the luciferase gene. The two variants (C/T) of rs1451385 were cloned in both directions in relation to the luciferase gene. As a positive control, we chose a human 411-bp enhancer from chromosome 12 that consistently increases luciferase expression 3-fold in a variety of cells lines (Cavalli et al., 2016).

We tested the four variants of Enh#385, the positive and empty controls in several cell lines. In all tested cell lines, we found a strong luciferase activity indicating a functional DNA fragment (Figures 2A–2D). The activity varied depending on the cell line and the Enh#385 variant used. The positive control enhancer increased luciferase activity to about 4-fold on average in all the cell lines tested, whereas Enh#385 increased luciferase activity up to 300-fold when compared with the empty vector. We observed an overall higher activity of the T allele, which is on the same haplotype as the WHR-increasing variant (Shungin et al., 2015).

Figure 2

Enhancer Tagged by GWAS Signal Shows High Activity and Response to Insulin

(A–D) rs1451385 shows allelic differences and is located in a strong enhancer. DNA fragments of 266 bp (Enh#385) containing the two allelic variants of rs1451385 (C or T) were tested for enhancer activity in (A) HepG2, (B) HEK293, (C) HeLa, and (D) SGBS human cell lines in transient reporter assays using 12–14 biological replicates. Both forward (f) and reverse (r) directions were tested along with a positive control (Cavalli et al., 2016) and empty vector as negative control.

(E and F) The enhancer is inducible by insulin and isoprenaline, compounds that have profound effects on adipocytes and are involved in insulin resistance. Reporter assays comparing luciferase expression after treatment of the indicated cell lines with insulin (100nM) and/or isoprenaline (3 μM) for 24 h. Data are presented as means ± SD. **p < 0.01, ***p < 0.001 comparing the C allele to T allele cloned in the same direction using a two-tailed unequal variance Student's t test. All the Enh#385 construct variants increased luciferase expression when compared with the empty vector (p < 0.01), and the positive control consistently showed a significant increase in expression when compared with the empty vector (p < 0.05). (E) Significant induction with insulin and isoprenaline was observed in HepG2 cells for all constructs when compared with non-treated cells (p < 0.001). (F) In SGBS cells, significant induction was only achieved with the constructs 385-C-f and 385-T-r (p < 0.05; Figures 2E and 2F).

The Enhancer Enh#358 Is Potentiated by Insulin and Isoprenaline

Next, we tested the transcriptional response of Enh#385 to insulin and isoprenaline in the reporter assay. Such a response would indicate a role of the enhancer in adipocyte biology and/or IR given that both compounds have profound effects on adipocytes and are involved in IR. Isoprenaline, also called isoproterenol, is a non-selective β-adrenoreceptor agonist that increases intracellular cAMP activity, stimulates lipolysis, and inhibits insulin-stimulated glucose transport, whereas insulin stimulates glucose uptake via Glut4 translocation and activation.

In both human hepatocytes (HepG2) and preadipocytes (SGBS), we observed a strong induction of luciferase activity with both substances (Figures 2E and 2F). The induction was suppressed with simultaneous treatment in HepG2 (Figure 2E). This suppression might reflect the opposite physiological effect of the lipogenic insulin and the lipolytic isoprenaline. At a molecular level, this type of suppression could be explained by competitive binding mechanisms, squelching, or transcriptional interference (Kamei et al., 1996, Levine and Manley, 1989, Manna and Stocco, 2007, Step et al., 2014, Yang-Yen et al., 1990, Zhang and Teng, 2001), resulting in insulin-induced factors competing with isoprenaline-induced factors for binding to Enh#385 (Figure S1). The formation of heterogeneous complexes could block the enhancer activity; a similar model has been proposed for scaffold protein complexes (Ferrell, 2000).

Molecular Dissection of the Enhancer

To define the boundaries and functional components of the enhancer, we used a luciferase assay to assess 20 independent constructs corresponding to different parts of Enh#385 with the T allele in forward direction (Figure 3A). Owing to poor transfection efficiency of SGBS cells, we restricted these experiments to HepG2 cells. To assess potential biotechnological applications of Enh#385, we included the enhancer of human CMV (cytomegalovirus), which is the most powerful enhancer currently used in commercial applications as a comparison. Based on phylogenetic conservation, we defined six sequential blocks (Figure 3A). Luciferase expression showed construct Enh#385-17 as having the highest enhancer activity (Figure 3B). The intensity of this construct was twice that of the original Enh#385 (Enh#385-11) but only 36% of the CMV enhancer. Furthermore, we tested the inducible properties of the most potent constructs (Figure 3C). Enh#385-11 showed the highest response to insulin and isoprenaline, whereas the CMV enhancer was not affected by either substance. In summary, our results indicate the presence of a core enhancer region centered around rs1451385.

Figure 3

Molecular Dissection of the Enhancer

(A) Restriction map of the constructs used for the luciferase assay. The original Enh#385 sequence and the phylogenetic blocks are marked at the bottom. The SNP rs1451385 is marked as a white square centered in Block5. CAGE transcription signals are indicated with green and purple arrows as in Figures 1B and 1C.

(B) Comparison of luciferase activity between constructs containing variations of Enh#385 and the CMV enhancer (E-CMV). The luciferase activity is displayed as percentage activity of E-CMV. The length of the fragments varies between the 62 bp of the construct Enh#385-20 to the 404 bp of the construct Enh#385-16; the E-CMV is 382 bp in length.

(C) Comparison of the inducible properties of the most potent enhancers from (B) and normalized to the luciferase activity of untreated parallel cultures.

CRISPR-Cas9 Mutation of the Enhancer Impairs Lipid Accumulation in Differentiating SGBS Cells

To study the function of the enhancer, we generated two independent mutations in undifferentiated preadipocytes using CRISPR-Cas9 (Figure S2).

After supplementation with adipogenic agents, we compared adipocyte differentiation of wild-type (WT) and mutated cultures after knockout (KO) of Enh#385 (KO-385-V11 and KO-385-V3). The WT cells differentiated to mature adipocytes in a normal fashion as evidenced by accumulation of lipid droplets, whereas the mutated cells showed highly reduced differentiation capacity (Figure 4A). Quantification of lipid droplets showed a 5-fold reduction of droplets in Enh#385 KO cells after 19 days of differentiation when compared with cells edited with a guide targeting a random intergenic region (Ctrl-24) or WT cells (Figure 4B). Additional experiments with independent CRISPR-Cas9 transductions provided similar results when measuring lipid content with oil red O staining (Ramirez-Zacarias et al., 1992) (Figure S3) or staining with the fluorophore BODIPY 493/503 (Majithia et al., 2014, Warnke et al., 2011) (data not shown).

Figure 4

Mutations of the Enhancer Impair Lipid Accumulation in Differentiating SGBS Cells

(A) Representative photographs of wild-type (WT) and mutated (KO-385-V3) differentiated adipocytes cultures. The faded areas in the pictures on the left are enlarged on the right, where the oil droplets are visible in the differentiated adipocytes. Scale bars, 200 μm.

(B) Quantification of adipogenesis by measuring oil droplets in the cell culture normalized by the total cell number after 19 days of differentiation. The KO-385-V11 and KO-385-V3 are cultures targeted with CRISPR-Cas9 upstream and downstream of the rs1451385, respectively, whereas Ctrl-24 targeted a random intergenic region. Data are presented as means ± SD.

(C) Quantitative PCR of WT and enhancer mutants for insulin receptor (INSR), phosphoenolpyruvate carboxykinase 1 (PCK1), lipoprotein lipase (LPL), and adiponectin (ADIPOQ). Data are presented as means ± SD.

Markers of Adipocyte Differentiation Are Reduced in Cell Cultures with Disrupted Enhancer

The impairment in differentiation was confirmed at the gene expression level by RT-qPCR of three adipocyte markers with divergent functions: the phosphoenolpyruvate carboxykinase 1 kinase (PCK1), involved in glucose metabolism; the lipoprotein lipase (LPL), involved in metabolism of fat; and the secreted adipokine adiponectin (ADIPOQ), involved in the control of fat metabolism and insulin sensitivity. In addition, insulin receptor (INSR) was included in this hypothesis-driven experiment owing to the potent induction by insulin of Enh#385. After 8 days of differentiation measured by RT-qPCR, all three adipocyte markers were strongly reduced and INSR was slightly reduced in the KO-385-V3 culture (Figure 4C; Figure S4); these results were subsequently confirmed by RNA-seq (next section).

RNA Sequencing Highlights Several Genes Showing Differential Expression after Mutation of the Enhancer

To understand the downstream consequences of the editing of Enh#385, we compared gene expression between WT and Enh#385-mutated cells after 8 days of induced differentiation by global RNA-seq. We included a culture transduced with a single-guide RNA targeting an unrelated intergenic region as a comparison to cope with the effect of transduction and the ectopic expression of the Cas9 nuclease. We focused on the most differentially expressed genes across the whole transcriptome to address downstream effects of Enh#385 and to help disentangle the role of this element on adipocyte differentiation (Figure 5A; Table S1). Among the top differentially expressed genes, we selected 20 genes for technical and biological validations (Table 1; Figure 5A). Technical validations were done by RT-qPCR (i.e., using a different technique than RNA-seq). The technical replication was high: differential gene expression assessed by RT-qPCR was confirmed for 19 of 20 genes when using cells mutated with the same CRISPR-Cas 9 vector that was used in cell cultures undergoing RNA-seq (downstream mutant [M(ds)] KO-385-V3). One gene, HYDIN, was impossible to consistently amplify by RT-qPCR, probably owing to its low expression in SGBS cells. Although the two methods correlated very well, we noted that especially for genes with low expression, the RNA-seq technology was more accurate than RT-qPCR. Two types of biological validations were done: first, using an independent mutation generated with a different CRISPR-Cas9 vector (KO-385-V11) targeting the sequence just upstream of rs1451385 (hence, annotated M(us) in Table 1) and second, using an innocuous lentivirus containing only GFP to assess a potential effect of the transduction protocol on gene expression (annotated as GFP in Table 1). Biological validation was confirmed for 11 of the remaining 19 genes when using the independent CRISPR-Cas9 upstream mutation, but when we evaluated the potential effects of the lentiviral transduction protocol per se on gene expression, we found that HEPH, NDN, and RCAN2 were downregulated in transduced cells (without the Enh#385 mutation), suggesting that the differential expression of these genes was, at least in part, due to the transduction protocol, rather than due to the mutation in Enh#385 (Table 1 and Figure 5B).

Figure 5

Differential Gene Expression in Differentiated Adipocytes

(A) Volcano plot comparing gene expression in differentiated adipocytes between wild-type and Enh#385 mutated cell cultures. Genes with significant differential expression (DE; padj < 0.05 and abs(log2 fold change) > 1) are in green, whereas those with abs(log2 fold change) > 1 are in orange. The 20 genes selected for subsequent validation by RT-qPCR are in blue, and the three genes subject to further investigation (CHI3L1, RELN, and NDN) are in red. Positive values mean that those genes are upregulated in mutants (RELN and PPAPDC1A), whereas negative values mean that they are downregulated. These results correspond to those shown in Table 1 (column 3).

(B) Comparison of relative gene expression from qPCR between differentiated SGBS adipocytes and undifferentiated SGBS preadipocytes for the 20 differentially expressed genes selected for validation with RT-qPCR using three biological replicates (Table 1). Genes wherein the expression was significantly changed more than 2-fold by the transduction protocol itself are marked with #. Adipocyte-specific genes, wherein the expression changed more than 2-fold between undifferentiated and differentiated cells, are marked with §. RCAN2 was both sensitive to transduction and adipocyte-specific and hence marked with both § and #. Data are presented as means ± SD.

Table 1

Validation and Filtering of Candidate Genes for Further Studies among 20 Top-Ranked Differentially Expressed Transcripts

		Differentiated SGBS Cells									Undifferentiated Cells						Comment
		RNA-Seq			RT-qPCR Validation						RT-qPCR
		M(ds)/WT			M(ds)/WT		M(us)/WT		GFP/WT		M(ds)/WT		M(us)/WT		GFP/WT
Gene Symbol	Full Gene Name	Log2FC	p Value*	Adjusted p	Log2FC	p Value*	Log2FC	p Value*	Log2FC	p Value*	Log2FC	p Value*	Log2FC	p Value*	Log2FC	p Value*
GAS7	Growth arrest-specific 7	−3.45	4.2 × 10−9	3.8 × 10−13	−3.19	1.2 × 10−2	0.67	2.1 × 10−1									Not validated M(us)
CHI3L1	Chitinase 3-like 1 cartilage glycoprotein-39	−2.03	1.9 × 10−8	5.7 × 10−6	−3.46	4.0 × 10−3	−1.57	8.0 × 10−3	0.88	3.9 × 10−1	−2.12	5.4 × 10−3	−1.05	2.4 × 10−2	−0.12	7.4 × 10−1	Putative causal gene
HEPH	Hephaestin	−2.87	2.3 × 10−7	1.8 × 10−5	−5.47	4.9 × 10−3	−2.24	1.4 × 10−2	−2.82	1.6 × 10−2							Downregulated by transduction
NDN	Necdin melanoma antigen MAGE family member	−2.15	1.9 × 10−6	1.4 × 10−4	−4.71	3.8 × 10−4	−1.88	4.2 × 10−3	−2.70	9.8 × 10−4							Downregulated by transduction
H2AFY2	H2A histone family member Y2	−2.62	2.0 × 10−6	3.4 × 10−2	−2.70	4.5 × 10−2	0.92	3.6 × 10−2									Not validated M(us)
MID2	Midline 2	−2.67	1.1 × 10−5	1.4 × 10−5	−4.35	3.4 × 10−4	0.60	6.8 × 10−2									Not validated M(us)
RCAN2	Regulator of calcineurin 2	−2.63	1.5 × 10−5	1.1 × 10−7	−4.19	2.5 × 10−2	−3.20	3.9 × 10−2	−3.35	3.1 × 10−2							Downregulated by transduction
RELN	Reelin	1.82	1.7 × 10−5	3.9 × 10−5	1.41	4.2 × 10−2	1.92	1.5 × 10−3	0.44	4.2 × 10−1	1.48	2.4 × 10−2	1.06	3.4 × 10−2	0.83	7.8 × 10−2	Putative causal gene
SGCD	Sarcoglycan delta 35 kDa dystrophin-associated glycoprotein	−3.43	2.4 × 10−5	5.5 × 10−5	−2.03	4.4 × 10−2	0.10	8.9 × 10−2									Not validated M(us)
ROR2	Receptor tyrosine kinase-like orphan receptor 2	−2.91	2.8 × 10−5	3.9 × 10−5	−1.39	5.2 × 10−2	0.63	5.0 × 10−2									Not validated M(us)
MYH2	Myosin heavy chain 2 skeletal muscle adult	−2.77	2.8 × 10−5	4.1 × 10−5	−1.65	3.9 × 10−2	0.47	3.7 × 10−1									Not validated M(us)
ITIH5	Inter-alpha-trypsin inhibitor heavy chain family member 5	−1.32	2.9 × 10−5	1.9 × 10−1	−3.17	4.6 × 10−5	−5.01	1.0 × 10−5	−1.56	8.1 × 10−3	Not expressed						Adipocyte specific
HYDIN	HYDIN axonemal central pair apparatus protein	−3.22	3.2 × 10−5	6.3 × 10−2	ND
PCK1	Phosphoenolpyruvate carboxykinase 1 soluble	−1.90	4.1 × 10−5	3.7 × 10−4	−2.03	7.7 × 10−3	−2.94	2.9 × 10−3	0.70	4.5 × 10−1	Not expressed						Adipocyte specific
MLXIPL	MLX-interacting protein-like	−1.67	5.3 × 10−5	8.6 × 10−1	−1.51	2.7 × 10−2	−0.73	5.3 × 10−2	−0.62	2.8 × 10−1	Not expressed						Adipocyte-specific
PLA2G2A	Phospholipase A2 group IIA platelets synovial fluid	−1.51	6.1 × 10−5	2.9 × 10−3	−1.32	2.2 × 10−3	−2.81	3.1 × 10−4	1.25	8.4 × 10−3	−0.17	8.0 × 10−1	−0.42	4.1 × 10−1	0.34	4.2 × 10−1	Adipocyte specific
PPAPDC1A	Phosphatidic acid phosphatase type 2 domain-containing 1A	1.31	7.9 × 10−5	2.7 × 10−1	2.18	4.7 × 10−4	1.45	3.8 × 10−2	0.55	3.0 × 10−1	−0.05	9.2 × 10−1	−0.75	5.3 × 10−2	−0.45	1.8 × 10−1	Preadipocyte specific
LPL	Lipoprotein lipase	−1.40	1.3 × 10−4	4.1 × 10−2	−2.04	2.5 × 10−9	−2.90	2.1 × 10−11	−0.07	4.7 × 10−1	Not expressed						Adipocyte specific
DYRK1B	Dual-specificity tyrosine-Y-phosphorylation regulated kinase 1B	−1.61	1.8 × 10−4	2.7 × 10−1	−2.43	8.9 × 10−6	−0.15	5.7 × 10−1									Not validated M(us)
INSR	Insulin receptor	−0.83	8.1 × 10−3	2.7 × 10−1	−0.89	2.6 × 10−2	0.56	8.4 × 10−2									Not validated M(us)

Genes were selected for confirmatory RT-qPCR analyses based on differential expression (DE) in RNA-seq analyses, excluding genes wherein DE was caused by the transduction procedure and genes with very low gene expression (see Transparent Methods supplemental file for details). The criteria for filtering out false-positives and selecting candidate genes for downstream experiments based on the RT-qPCR validation were three-fold: (1) A technical validation using RNA from the same cells as those used for the global transcriptome analysis requiring significant differentially expressed in the same direction as in RNA-seq for successful validation (columns 6 and 7 [M(ds)/WT]). (2) A biological validation using RNA from an independent mutation (columns 8 and 9 [M(up)/WT]). (3) A control of potential transduction effects in which the candidate genes should not show differential expression comparing non-transduced WT cells and cells transduced with a lentivirus expressing only GFP (columns 10 and 11 [GFP/WT]). Validated genes were further interrogated for possible involvement in adipogenesis by studying their expression in undifferentiated cells. Genes not expressed in preadipocytes or with a >2-fold difference in expression between the differentiated and undifferentiated state were considered adipocyte specific (Figure 5B), and likely a consequence of the phenotype and not causative genes. RT-qPCR data of three independent experiments; p values comparing expression in mutated cells with WT cells (as log2 fold change [FC]) were calculated by Student's t test with Benjamini-Hochberg correction (columns 7, 9, 11, 13, 15, and 17). * P-values are from two-sample Welch t-tests, unadjusted for multiple testing.

Next, to investigate a potential role in adipogenesis of the eight differentially expressed genes passing both technical and biological validation, we examined their expression in undifferentiated preadipocytes (Table 1). We assumed that only genes that were differentially expressed already in preadipocytes could be responsible for the impairment in adipogenesis, whereas a large number of differentially expressed genes represent the consequence of altered adipogenesis of the mutated cultures rather than causal genes. After comparisons of gene expression by RT-qPCR between undifferentiated downstream mutants, upstream mutants, GFP-transduced cells, and WT cells, we found that only chitinase-3-like protein 1 (CHI3L1) on chromosome 1 and reelin (RELN) on chromosome 7 were differentially expressed in undifferentiated preadipocytes (at least 2-fold difference; Table 1). Several genes with strong prior evidence for involvement in metabolic disease, such as ITIH5, PCK1, MLXIPL, and LPL, were only expressed in mature adipocytes, and hence unlikely to be causal in driving the decrease in adipogenesis, but rather be a result of the difference in cell fate between the Enh#385 KO and WT cells (Table 1; Figure 5B). In subsequent experiments using global gene expression analyses, we were able to extend our search of causal genes and identified altered pathways important for adipogenesis (see below, Table S4 and Figure 7).

Figure 7

A Schematic Illustration of Key Adipogenic Factors that We Propose as Being Regulated by the Enhancer (Enh#385)

Arrows with dotted lines annotate predictions based on the present work, whereas solid arrows indicate associations experimentally demonstrated and reported in prior literature. Arrowhead lines indicate transcriptional activation; blunt-head lines indicate inhibition. Gradient arrowheads after the name of the gene indicate the direction of the expression measured in our mutated cultures when compared with wild-type cell cultures. The upstream regulation of CHI3L1 is uncertain and indicated with a question mark.

Ablation of CHI3L1 Decreases Adipogenesis in Differentiating SGBS Cells

To evaluate the direct role of the two consistently differentially expressed genes on adipogenesis, we used CRISPR-Cas9 to disrupt CHI3L1 and RELN in SGBS cells (Figure S9). After induction of adipogenesis, the CHI3L1-KO cells showed a significant reduction in lipid accumulation, whereas RELN-KO cells did not differ when compared with controls with regard to differentiation capacity (Figure 6). The phenotype similarity obtained by direct mutation of CHI3L1 and reduction of CHI3L1 via mutation of Enh#385 indicates that CHI3L1 and the enhancer are linked by functional interaction circuits. RELN was upregulated in Enh#385 mutant cells; therefore an increase of lipid accumulation in RELN-KO cells could be expected. We did not detect any differences in lipid accumulation between the WT and RELN-KO cells, but we cannot rule out that the overexpression of this gene could reduce the rate of adipocyte differentiation.

Figure 6

Adipocyte Differentiation Capacity of Wild-Type and Mutant SGBS Cells Estimated as Accumulation of Neutral Lipids

(A) Representative images of BODIPY and Hoechst staining, in green (lipids) and blue (nuclei), respectively, of differentiating SGBS cells. From left to right: wild-type (WT), CRISPR-edited CHI3L1 (CHI3L1-KO), and CRISPR-edited RELN (RELN-KO) cells.

(B) Quantification of neutral lipids after 8 days of differentiation.

(C) Quantification of neutral lipids after 12 days of differentiation. Both at 8 and 12 days after induction of differentiation, the lipid accumulation was significantly reduced in CHI3L1 mutants (p < 0.001), whereas RELN mutants were no different from WT. The values in (B and C) represent the mean ± SD of a minimum of three cell cultures.

(D) Relative CHI3L1 expression comparing wild-type cultures and CRISPR-edited SGBS cells (knocking down CHI3L1).

(E) Relative RELN expression comparing wild-type cultures and CRISPR-edited SGBS cells (knocking down RELN). The data in (D and E) are presented as mean ± SD.

CHI3L1, also known as YKL-40, is a secreted glycoprotein coupled with stress-induced cellular responses (Ling and Recklies, 2004). Protein levels of CHI3L1 have been associated with several pathogenic processes including schizophrenia, asthma, obesity, and cancer, but the biological function of YKL-40 in specific tissues is largely unknown (Kyrgios et al., 2012, Ober et al., 2008, Zhao et al., 2007). Experimental data have shown that mice deficient in Chi3l1 develop less visceral obesity and have smaller adipocytes; in contrast, its overexpression induced adiposity (Ahangari et al., 2015). The precise molecular mechanism of CHI3L1 on adipogenesis has not been determined, but a direct role on the extracellular matrix through type I collagen has been proposed (Iwata et al., 2009, Mariman and Wang, 2010). It is possible that CHI3L1 is required for the remodeling of the extracellular matrix that precedes the adipocyte differentiation process. In a mouse model of osteomyelitis, partial restoration of osteogenesis has been achieved by suppressing Chi3l1 (Chen et al., 2017). This suggests that CHI3L1 might play a key role in cell fate decision (Figure 7).

Enh#385 Regulates Genes Involved in Adipocyte Biology and Nearby Genes

We specifically examined the expression of genes linked to adipocyte biology in our RNA-seq data from differentiated cells. We confirmed the reduction in expression of the markers previously measured by RT-qPCR and observed reduction of several additional adipocyte markers. For example, LEP, PPARGC1A, PLIN1, CD36, and SLC2A4 (GLUT4) were reduced in mutant differentiated cultures. We noted that expression of some key proteins involved in fat synthesis (FASN), transport of lipoproteins (APOE), and the receptors for ADIPOQ (ADIPOR1 and ADIPOR2) were unaffected by the mutations of the enhancer (Table S2).

Next, we explored whether there were any putative cis-regulated genes among the differentially expressed genes. We scanned the chromosomal region surrounding rs1451385 (±5Mb) for differentially expressed genes. All genes annotated in the UCSC Browser (https://genome.ucsc.edu/), including protein-coding genes and non-coding RNA genes within this interval, were identified and expression levels were compared between mutant and WT cells. This region included a microRNA (miRNA) (MIR148A) and a few genes that are good candidates for involvement in metabolic disease based on prior literature, such as IGF2BP3 (insulin-like growth factor 2 mRNA-binding protein 3) at chr7:23,349,828 (7p15.3) and CREB5 (cAMP responsive element-binding protein 5) at chr7:28,865,511 (7p14.1). The only transcript with a significant differential expression in this large region was NFE2L3 (2.3-fold reduction, p = 0.009). However, as NFE2L3 is expressed at very low levels in SGBS cells, and because we did not find any differential expression in similar experiments using undifferentiated SGBS cells (see below), we did not consider this association further. These results suggest that the enhancer does not act in cis on protein-coding genes, but it should be noted that our global RNA sequence method did not allow quantification of short RNAs, such as miRNAs.

Causative Genes Regulated by Enh#385 Uncovered by Differential Gene Expression in Undifferentiated Cells

To extend the list of potential downstream causal genes regulated by Enh#385, we performed another transcriptomics analysis by RNA-seq, this time in undifferentiated SGBS preadipocytes. We compared the two independent mutants of Enh#385 (KO-385-V3 and KO-385-V11) with WT cells and a GFP transduction control (Figure S5; Table S3). In these analyses, we confirmed the downregulation of CHI3L1 (fold change, 0.49; padj, 0.0002) and the upregulation of RELN (fold change, 1.55; padj, 0.041) in mutated cells corroborating a potential causative role in driving the decreased adipogenesis noted after mutation of Enh#385. However, RELN was also upregulated in GFP-transduced cells, which may indicate that this differential expression was, at least partly, driven by the transduction procedure. In addition, NDN (Necdin; fold change, 0.42; padj, 0.018), associated with the Prader-Willi syndrome in humans, characterized by severe obesity, showed significant DE in undifferentiated preadipocytes, indicating a role as effector transcript downstream of Enh#385. This gene was discarded in the first qPCR validation analysis due to concerns regarding unspecific differential expression induced by the transduction protocol. Comparison of differentiated and undifferentiated WT cells revealed several adipocyte-specific genes (Figure 5B). These genes show higher (or even exclusive) expression in differentiated adipocytes; hence, it is likely that their differential expression in the Enh#385 mutants is a consequence of the lack of adipogenesis, rather than being involved in the causative process.

The most differentially expressed gene was H19, a maternally imprinted lincRNA that interacts with insulin-like growth factor 2 (IGF2), and has been suggested to be a tumor suppressor (Bartolomei et al., 1991, Bell and Felsenfeld, 2000). Defects in the H19/IGF2 imprinting have been associated with Beckwith-Wiedemann syndrome (DeBaun et al., 2002), Silver-Russell syndrome (Bartholdi et al., 2009, Bliek et al., 2006), and Wilms tumor 2 (Steenman et al., 1994). Other top differentially expressed genes included several plausible candidates for involvement in adiposity and/or IR, such as ADH1B, IL8, ACE, IL1B, and IL6 (Figure S5; Table S3).

Next, we focused on genes linked to adipocyte biology and adipogenesis based on prior literature (Table S4). We found that PPARG was reduced (0.7-fold) in mutant cells, whereas WNT10B, encoding a molecular switch that inhibits adipogenesis, was significantly increased (2.5-fold) (Christodoulides et al., 2006, Ross et al., 2000). It has been shown that WNT10B inhibits PPARG and promotes the expression of the osteoblastogenic transcription factor RUNX2 (Bennett et al., 2005). Consistent with a potential pivotal role of WNT10B in the observed phenotype, we found an increase (1.3-fold) of RUNX2 in the mutated cells. The anti-adipogenic transcription factor KLF2 has also been suggested as contributor in this pathway (Banerjee et al., 2003). Again, consistent with this hypothesis, KLF2 was significantly upregulated (1.4-fold) in our enhancer mutants, thus blocking adipogenesis. The increase of these two important anti-adipogenic regulators, WNT10B and KLF2, acting upstream of PPARG suggests that Enh#385 might function very early in the initiation of the adipogenesis program.

Discussion

Our study offers a detailed characterization of a GWAS locus (previously annotated as SNX10) associated with WHR and other obesity-related traits in a series of functional experiments. Our main conclusions are severalfold. First, we identified an enhancer active during adipogenesis that responds strongly to external metabolic stimuli, such as insulin and isoprenaline. The known WHR-raising allele (A at rs3902751) (Shungin et al., 2015) is on the same haplotype as the allele being associated with higher enhancer activity (T at rs1451385), suggesting that higher activity of this enhancer is associated with higher WHR. Second, mutation of the enhancer using CRISPR-Cas9 in preadipocytes dramatically impairs adipocyte differentiation. RNA-seq followed by replication using RT-qPCR in mature adipocytes demonstrated that several adipocyte genes were downregulated in cells with disrupted Enh#385 and highlighted CHI3L1 as a potential downstream target, also downregulated. Third, disruption of CHI3L1 caused a decrease in adipogenesis.

Independently of the allele tested, the enhancement of luciferase expression was higher with Enh#385 cloned in forward direction (Figure 2), suggesting that the cloned DNA fragments may have a significant promoter-like activity. Consistent with this, the signal from the CAGE-sequencing experiment was generally stronger in the forward direction throughout adipogenesis (Figure 1D). However, given the long physical distance to the nearest gene, the regulatory element is most likely to be an enhancer with promoter-like activity. This complex nature of the regulatory element has been observed in a previous study, wherein the sequence of Enh#385 was included in a DNA fragment (976 bp) functionally classified as weak anti-repressor element (Kwaks et al., 2003).

CHI3L1 could be a downstream effector of Enh#385 having a role in adipocyte differentiation, but the mechanistic circuits between the Enh#385 and CHI3L1 function are unclear, especially as they are not located on the same chromosome. We hypothesized that the mediator of the link from Enh#385 to CHI3L1 and other downstream genes involved in adipogenesis could be MIR148A, a non-protein-coding gene encoding an miRNA localized 95 kb from the enhancer (Figure 1). miRNAs are known to be master regulators of groups of proteins (Stefani and Slack, 2008), and indeed, several miRNAs have been reported to regulate adipogenesis and lipid metabolism (Ahn et al., 2013, Esau et al., 2004, Yang et al., 2011). Indeed, our analysis of enhancer-promoter interactions by HiChiP in a human coronary artery smooth muscle cell line shows connections between the Enh#385 and MIR148A (Figure S6).

Searching for MIR148A targets using TargetScan (http://www.targetscan.org), we found two conserved sequences at the 3′ UTR of WNT10B and two sequences of KLF2 that acts as binding targets for miR-148a (Figure S7). Consistent with this, we observed a significant increase of WNT10B in Enh#385 mutants (Table S4). There is also support in prior literature for involvement of MIR148A and WNT10B in adipogenesis. In mice, miR-148a is upregulated during adipogenesis and downregulated in mature adipocytes of obese animals (Xie et al., 2009). In a recent report using differentiating human adipose-derived mesenchymal stem cells, it was shown that MIR148A induces adipogenesis by suppressing its target gene, WNT10B, an endogenous inhibitor of adipogenesis (Shi et al., 2015). It is tempting to speculate that the increase of WNT10B and KLF2, two anti-adipogenic factors, in our mutants is the result of the reduction of MIR148A (Table S4). This scenario is consistent with the observed increase of RUNX2 and decrease of PPARG in our DE data and with the lack of adipogenesis in the Enh#385 mutant cells. A schematic illustration of this proposed model is shown in Figure 7.

Experimentally, we were able to measure an increase of expression of MIR148A between 0 and 8 days of differentiation, but we could not detect significant differences in expression of MIR148A between WT and Enh#385 KO cells before or 8 days after induction of differentiation. It is possible that a role of MIR148A in regulating adipogenesis is transient, and can only be recorded in a short time window during the differentiation process (and we analyzed differential expression at just two time points). Supporting this interpretation, CAGE-sequencing data (from the FANTOM consortium) measuring the activity of the MIR148A promoter during MSC adipocyte differentiation show only a consistent increase in activity at 4 and 14 days, but not at the rest of the 17 time points, indicating a fluctuating and transient regulation of MIR148A (Figure S8). The inconsistent results in detection of MIR148A in mutant cells could be also due to the negative regulatory feedback loop between MIR148A and one of its targets, the methylase DNMT1 (Figure 7). The methylase downregulates the expression of MIR148A by hypermethylation of its promoter (Hong et al., 2018, Long et al., 2014).

The consistent downregulation of the imprinted gene NDN (Necdin) leads us to hypothesize that imprinting by methylation of NDN could be mediated by the DNMT1 (DNA methyltransferase 1)-MIR148A circuit. NDN is a particularly interesting gene for adipocyte biology as it is one of the candidate genes in the chromosomal region 15q11-15q13, which results in Prader-Willi syndrome when the paternal copy is deleted. In addition to mild to moderate intellectual impairment and behavioral problems, obesity and T2D are the most common symptoms of Prader-Willi syndrome, owing to insatiable appetite and chronic overeating (MacDonald and Wevrick, 1997). DNMT1 is a predicted target of MIR148A, and its expression is elevated in adipocytes of obese individuals (Braconi et al., 2010). DNMT1 maintains methylation pattern during adipocyte differentiation, and its silencing accelerated adipogenesis (Londono Gentile et al., 2013). It has been proposed that DNMT1 and MIR148A are regulated by a negative feedback loop that could explain the discrepancy in the published reports (Xu et al., 2013). In fact, the promoter of MIR148A is in proximity of CpG islands, and it is silenced by hypermethylation (Hanoun et al., 2010). Our differential expression analysis shows that DNMT1 expression is only slightly increased in mutant cultures; however, it has been reported that the targeting of this gene by MIR148A occurs exclusively at the protein level (Pan et al., 2010), which makes the differences in mRNA a less ideal indicator of enzyme methylation activity. Imprinted genes (identified using the Genomic Imprinting website: http://www.geneimprint.com) were generally not affected by the mutations of Enh#385. However, in addition to NDN, the adipogenesis-linked and imprinted transcript H19 was the most differentially expressed gene in preadipocytes, showing a consistent downregulation in mutated cells (Han et al., 2018, Huang et al., 2016). It is plausible that only loci involved in adipogenesis are accessible to DNMT1, so it should not be expected that all imprinted genes and genes regulated by CpG methylation are downregulated in Enh#385 mutants. The increased activity of this enzyme may lead to the repression of active adipogenic loci controlled by methylation. Our expression analysis shows that the imprinted gene NDN and genes regulated by methylation such as DLK2 (Delta like non-canonical Notch ligand 2) and EBF2 (Early B-cell factor 2), both promoting adipogenesis (Jimenez et al., 2007, Nueda et al., 2007), are consistently downregulated in cells with disrupted Enh#385 (Tables S2 and S4). Thus, despite being limited to a few genes, this methylation link reinforces our model in which MIR148A is the most likely target of the enhancer Enh#385 (Figure 7).

GWAS have been remarkably successful in discovering loci associated with complex traits, but to unlock the transformative potential of these findings, there is a strong need for detailed studies of the molecular mechanisms underlying these signals. We have performed such a study that improves the understanding of adipocyte biology and fat distribution. Although the exact molecular events that link genetic variation in Enh#385 (previously annotated as the SNX10 locus) with the impairment of adipocyte differentiation still need to be characterized in detail, one potential explanation for link of the enhancer with adipogenesis could be cis-acting regulation through the nearby MIR148A. We hypothesize that allelic variants that reduce the activity of Enh#385 could lead to the downregulation of MIR148A. This reduction of MIR148A may produce an increase of the cell fate determinant WNT10B and the subsequent inhibition of adipogenesis. These initial and complex alterations of gene activity could lead to the reduction of CHI3L1 expression, which, as we demonstrate in this work, is essential for the adipocyte maturation. Further studies are needed to establish these potential mechanisms linking Enh#385 to adipocyte differentiation.

Limitations of the Study

The model we propose in which the enhancer Enh#385 acts on MIR148A is based on indirect observations and is not directly proved in this study. Hence, it should be viewed as hypothesis generating, and we do not discard the possibility that the enhancer also acts on other genes or genetic elements.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.