Intergenic GWAS SNPs are key components of the spatial and regulatory network for human growth

Meta-analysis of genome-wide association studies has resulted in the identification of hundreds of genetic variants associated with growth and stature. Determining how these genetic variants influence growth is important, but most are non-coding, and there is little understanding of how these variants contribute to adult height. To determine the mechanisms by which human variation contributes to growth, we combined spatial genomic connectivity (high-throughput conformation capture) with functional (gene expression, expression Quantitative Trait Loci) data to determine how non-genic loci associated with infant length, pubertal and adult height and contribute to gene regulatory networks. This approach identified intergenic single-nucleotide polymorphisms (SNPs) ∼85 kb upstream of FBXW11 that spatially connect with distant loci. These regulatory connections are reinforced by evidence of SNP-enhancer effects and altered expression in genes influencing the action of human growth hormone. Functional assays provided evidence for enhancer activity of the intergenic region near FBXW11 that harbors SNP rs12153391, which is associated with an expression Quantitative Trait Loci. Our results suggest that variants in this locus have genome-wide effects as key modifiers of growth (both overgrowth and short stature) acting through a regulatory network. We believe that the genes and pathways connected with this regulatory network are potential targets that could be investigated for diagnostic, prenatal and carrier testing for growth disorders. Finally, the regulatory networks we generated illustrate the power of using existing datasets to interrogate the contribution of intergenic SNPs to common syndromes/diseases.

Introduction

The causes of clinically relevant deviations in growth and growth hormone (GH) levels remain largely unexplained but often are associated with known syndromes, either as a direct cause or as a symptom of a larger problem (1,2). Normal variation in human growth primarily results from changes in the production of, secretion of and sensitivity to GH. This variation is further compounded by fluctuations associated with the hormones released from the hypothalamus (e.g. GH-releasing hormone or somatostatin to regulate GH secretion), thyroid (e.g. thyroid hormone to alter GH response in the hypothalamus), liver (e.g. IGF-1 to inhibit GH secretion) and gastrointestinal tract (e.g. ghrelin to stimulate GH secretion) (1,3). There is a wide variability in growth response to GH treatment in short children (e.g. 3-M Syndrome (4)) and no gold standard to predict treatment responsiveness (5). Moreover, most patients do not achieve clinical treatment goals using standard treatments, irrespective of the treatment method that is chosen (1).

For individuals predisposed to clinically relevant deviations in height (short stature or overgrowth), lack of (or resistance to) treatment can be costly. These costs are both financial and associated with an increased lifetime risk of diseases such as cancer or diabetes (1,2). Moreover, unrestrained GH and IGF1 levels (i.e. acromegaly) are associated with a 3-fold higher mortality (1). Finally, many growth syndromes are associated with mental disabilities (6). Therefore, we need a more detailed understanding of the genetic regulation of stature to enable clearer diagnoses and better inform the management (e.g. predict GH treatment failure) of children with growth disorders (2).

Rapid increases in our understanding of the genetics of stature have been obtained using genome-wide association studies (GWAS), which test for statistical associations between genotype (i.e. common gene variants or single-nucleotide polymorphisms (SNPs)) and phenotype. Three recent GWAS meta-analyses have investigated the regulation of human growth at different biological time-points during infancy (7), puberty (8) and adulthood (9). Many of the SNPs that are associated with growth do not occur in exons or promoters, making it difficult to assign biological relevance to these intergenic SNPs (10,11).

Fine-mapping frequently validates intronic and intergenic SNPs as leading risk-associated variants (12), and it is therefore important to identify their function. Intergenic SNPs are typically assigned a theoretical pathogenicity based on their proximity to nearby genes (13). For example, in the analysis by Wood et al. (9), SNPs that were within a 1 Mb window were combined into a single locus that was reported as being associated with adult height. More recently, van der Valk et al. incorporated cis-expression Quantitative Trait Loci (eQTLs) data into their analysis to refine the SNP-gene attributions they made for three novel SNPs associated with pubertal height. However, in this same meta-analysis, van der Valk et al. (7) also classified other intergenic variants based on the functions of nearby genes in the absence of additional supportive (e.g. eQTL) evidence.

Without additional supportive evidence, proximity-based SNP-gene attributions are problematic. This is particularly true of SNPs associated with stature, which are significantly enriched for regulatory motifs (9). The problems arise because genetic variants associated with regulatory elements are often connected to the activity of spatially associated (physically coupled) distal genes rather than with nearby genes (14,15). For example, SNPs that are strongly associated with diabetes and obesity lie in 89 common intronic variants of the FTO gene (16). However, the lead SNP (rs9930506) has been shown to participate in spatial associations that affect regulation of the downstream IRX3 gene (Fig. 1a (17)). A similar situation has been observed for the rs1736020, rs1297265 and rs2823286 SNPs within an intergenic region of chromosome 21. In this instance, these SNPs interact with the NRIP1 gene, located 380 kb away and not the closer anti-inflammatory gene USP25 (Fig. 1b (18)). Collectively, these examples highlight the potential for drawing misleading conclusions from automatically associating the SNP with the closest gene.

Figure 1.

SNPs associated with (a) diabetes and (b) inflammatory bowel disease do not interact with the closest gene. (a) rs9930506 has been shown to spatially interact with and regulate the distal IRX3 gene (17). (b) Three intergenic SNPs located closer to the downstream USP25 and associated with inflammatory bowel disease are instead shown to interact and regulate the distal upstream NRIP1 gene (18).

There is no a priori reason why spatial regulation requires SNPs and genes to be adjacent (in cis-proximity, <1 Mb in distance) within the chromosome sequence. In fact, spatial associations between regulatory elements and promoters on different chromosomes contribute to gene regulation (18), as captured by proximity ligation (Chromosome Conformation Capture, Hi-C, and other related methodologies; e.g. (18–20)). Therefore, since the majority of SNPs (i.e. those that are intergenic) are, by default, regulatory in nature (i.e. not causing direct alterations in DNA coding potential), it is probable that the effects of SNPs are not simply directed to genes that are close in the linear sequence (14,21). Rather, long-range (trans, >1 Mb in distance or interchromosomal) promoter contacts are often associated with the tag SNPs, indicating these interactions as a more likely mechanism to explain the contribution that the genetic variants make to the phenotype (18).

Growth and development are particularly important processes throughout which long-range spatial regulatory interactions between promoters and enhancers are required to obtain precisely regulated patterns of gene expression (15,22,23). Modulation of this activity by genetic polymorphisms within the enhancer regions (i.e. short genomic regions that are bound by transcription factors) can be associated with expression changes in the target gene(s) (i.e. an expression quantitative trait locus [eQTL]). These trans eQTLs can arise by two main modes of action: (1) alterations in enhancer activity directly changing the regulation of physically associated gene(s) in trans and (2) alterations to transcription factor binding causes changes to the regulation of gene(s) in cis that are epistatic to the genes in trans. Mechanism (2) includes direct or indirect changes to the expression of genes caused by enzymatic feedback/forward cycles. It is likely that these enhancer-promoter interactions are highly tissue specific, making them capable of contributing to both spatial and temporal control of gene expression (24). Thus, SNPs that alter enhancers are likely to play significant roles in diverse processes and tissues across time.

We previously incorporated spatial relationships into an analysis linking SNPs associated with diabetes to potential functional outcomes (21). Here, we hypothesized that the inclusion of evidence of trans-interactions into an analysis of growth-associated SNPs (infant length (7), pubertal (8) and adult height (9)) would provide an improved explanation for genetic regulation of stature.

Results

We screened 722 significant GWAS meta-analysis SNPs from three meta-analyses (86 underlying GWAS studies; Supplementary Material, Table S1, Fig. 2) for their spatial connections using existing proximity ligation data from 16 different cell types and 73 regulatory factors motifs from ENCODE (7–9). Six-hundred seventy of the 722 SNPs are intronic (360) or intergenic (310), meaning they do not directly affect coding regions. Additionally, Haploreg identified 490 of the 722 SNPs (Supplementary Material, Table S2) as falling in regions that were significantly enriched for enhancer activity in a broad list of ENCODE tissue types (Table 1). Therefore, it is almost certain that a majority of the SNPs that are associated with growth are regulatory in nature, prompting us to study these SNPs further.

Figure 2.

Genome3D Pipeline for GWAS SNP spatial analysis. GWAS SNPs are filtered through spatial and eQTL data to determine SNPs with high likelihood of altering long-distance enhancers. SNPs are then validated as enhancers by luciferase activity testing in in vitro human cell cultures.

Table 1.

Enrichment for histone modifications associated with enhancer activity surrounding the 722 growth associated SNPs

ENCODE tissue type	Enhancer activity within set of 722 SNPs		Binomial P
	Observed	Expected
E017 LNG.IMR90 (IMR90 fetal lung fibroblasts cell line)	92	40.1	<1.00E-06
E006 ESDR.H1.MSC (H1 derived mesenchymal stem cells)	75	34.8	<1.00E-06
E026 STRM.MRW.MSC (bone marrow derived cultured mesenchymal stem cells)	75	30.7	<1.00E-06
E049 STRM.CHON.MRW.DR.MSC (mesenchymal stem cell derived chondrocyte cultured cells)	84	38.3	<1.00E-06
E025 FAT.ADIP.DR.MSC (adipose derived mesenchymal stem cell cultured cells)	107	46.6	<1.00E-06
E023 FAT.MSC.DR.ADIP (mesenchymal stem cell derived adipocyte cultured cells)	82	35.2	<1.00E-06
E052 MUS.SAT (muscle satellite cultured cells)	72	34	<1.00E-06
E055 SKIN.PEN.FRSK.FIB.01 (foreskin fibroblast primary cells skin01)	90	38.2	<1.00E-06
E056 SKIN.PEN.FRSK.FIB.02 (foreskin fibroblast primary cells skin02)	65	28	<1.00E-06
E061 SKIN.PEN.FRSK.MEL.03 (foreskin melanocyte primary cells skin03)	69	33.1	<1.00E-06
E028 BRST.HMEC.35 (breast variant human mammary epithelial cells (vHMEC))	70	37.4	<1.00E-06
E027 BRST.MYO (breast myoepithelial primary cells)	94	47.1	<1.00E-06
E054 BRN.GANGEM.DR.NRSPHR (ganglion eminence derived primary cultured neurospheres)	47	20.8	<1.00E-06
E053 BRN.CRTX.DR.NRSPHR (cortex derived primary cultured neurospheres)	61	29	<1.00E-06
E063 FAT.ADIP.NUC (adipose nuclei)	74	38.1	<1.00E-06
E108 MUS.SKLT.F (skeletal muscle female)	87	35.5	<1.00E-06
E107 MUS.SKLT.M (skeletal muscle male)	92	35.3	<1.00E-06
E089 MUS.TRNK.FET (fetal muscle trunk)	76	38	<1.00E-06
E090 MUS.LEG.FET (fetal muscle leg)	100	51.7	<1.00E-06
E083 HRT.FET (fetal heart)	86	42.2	<1.00E-06
E104 HRT.ATR.R (right atrium)	54	23.7	<1.00E-06
E111 GI.STMC.MUS (stomach smooth muscle)	45	19.7	<1.00E-06
E092 GI.STMC.FET (fetal stomach)	69	34	<1.00E-06
E088 LNG.FET (fetal lung)	80	36.1	<1.00E-06
E080 ADRL.GLND.FET (fetal adrenal gland)	91	44.3	<1.00E-06
E066 LIV.ADLT (liver)	66	33.8	<1.00E-06
E098 PANC (pancreas)	58	25.5	<1.00E-06
E114 LNG.A549.ETOH002.CNCR (A549 EtOH 0.02pct lung carcinoma cell line)	72	27.3	<1.00E-06
E118 LIV.HEPG2.CNCR (HepG2 hepatocellular carcinoma cell line)	86	40.3	<1.00E-06
E119 BRST.HMEC (HMEC mammary epithelial primary cells)	72	35.5	<1.00E-06
E120 MUS.HSMM (HSMM skeletal muscle myoblasts cells)	65	31.4	<1.00E-06
E121 MUS.HSMMT (HSMM cell derived skeletal muscle myotubes cells)	63	28.3	<1.00E-06
E122 VAS.HUVEC (HUVEC umbilical vein endothelial primary cells)	58	26.7	<1.00E-06
E126 SKIN.NHDFAD (NHDF-Ad adult dermal fibroblast primary cells)	83	35.6	<1.00E-06
E128 LNG.NHLF (NHLF lung fibroblast primary cells)	68	26.4	<1.00E-06
E129 BONE.OSTEO (osteoblast primary cells)	86	31.8	<1.00E-06
E076 GI.CLN.SM.MUS (colon smooth muscle)	55	27.6	1.00E-06
E097 OVRY (ovary)	49	23.3	1.00E-06
E040 BLD.CD4.CD25M.CD45RO.MPC (primary T helper memory cells from peripheral blood 1)	49	24	2.00E-06
E078 GI.DUO.SM.MUS (duodenum smooth muscle)	40	17.8	2.00E-06
E091 PLCNT.FET (placenta)	76	44.2	3.00E-06
E057 SKIN.PEN.FRSK.KER.02 (foreskin keratinocyte primary cells skin02)	62	33.8	4.00E-06
E068 BRN.ANT.CAUD (brain anterior caudate)	50	25.3	5.00E-06
E081 BRN.FET.M (fetal brain male)	48	24	5.00E-06
E105 HRT.VNT.R (right ventricle)	50	25.3	5.00E-06
E113 SPLN (spleen)	57	31	8.00E-06
E043 BLD.CD4.CD25M.TPC (primary T helper cells from peripheral blood)	51	27.1	1.40E-05
E125 BRN.NHA (NH-A astrocytes primary cells)	52	27.8	1.40E-05
E059 SKIN.PEN.FRSK.MEL.01 (foreskin melanocyte primary cells skin01)	45	22.8	1.50E-05
E051 BLD.MOB.CD34.PC.M (primary hematopoietic stem cells G-CSF-mobilized male)	59	33.4	1.90E-05
E079 GI.ESO (esophagus)	41	20.2	1.90E-05
E013 ESDR.CD56.MESO (hESC-derived CD56+ mesoderm cultured cells)	55	30.5	2.20E-05
E110 GI.STMC.MUC (stomach mucosa)	54	29.8	2.20E-05
E047 BLD.CD8.NPC (primary T CD8+ naive cells from peripheral blood)	49	26.3	2.90E-05
E099 PLCNT.AMN (placenta amnion)	45	23.4	3.00E-05
E067 BRN.ANG.GYr (brain angular gyrus)	41	20.6	3.20E-05
E093 THYM.FET (fetal thymus)	59	34.2	3.60E-05
E100 MUS.PSOAS (psoas muscle)	37	18.1	4.30E-05
E094 GI.STMC.GAST (gastric)	37	18.1	4.30E-05
E085 GI.S.INT.FET (fetal intestine small)	57	32.9	4.50E-05
E001 ESC.I3 (ES-I3 cells)	53	30.1	5.50E-05
E123 BLD.K562.CNCR (K562 leukemia cells)	45	24.1	5.50E-05
E077 GI.DUO.MUC (duodenum mucosa)	45	24.6	9.20E-05
E041 BLD.CD4.CD25M.IL17M.PL.TPC (primary T helper cells PMA-I stimulated)	48	26.9	9.30E-05
E022 IPSC.DF.19.11 (iPS DF 19.11 cells)	58	34.8	0.000108
E095 HRT.VENT.L (left ventricle)	53	30.9	0.000108
E109 GI.S.INT (small intestine)	28	12.7	0.000113
E073 BRN.DL.PRFRNTL.CRTX (brain dorsolateral prefrontal cortex)	39	20.4	0.000115
E086 KID.FET (Fetal Kidney)	33	16.3	0.000129
E005 ESDR.H1.BMP4.TROP (H1 BMP4-derived trophoblast cultured cells)	58	35.1	0.000133
E038 BLD.CD4.NPC (primary T helper naive cells from peripheral blood)	40	21.4	0.00015
E050 BLD.MOB.CD34.PC.F (primary hematopoietic stem cells G-CSF-mobilized Female)	59	36.1	0.000155
E065 VAS.AOR (aorta)	25	11.1	0.000182
E030 BLD.CD15.PC (primary neutrophils from peripheral blood)	43	24.1	0.000215
E074 BRN.SUB.NIG (brain substantia nigra)	48	27.9	0.000216
E102 GI.RECT.MUC.31 (rectal mucosa donor 31)	43	24.2	0.000233
E020 IPSC.20B (iPS-20b cells)	39	21.3	0.000252
E096 LNG (lung)	42	23.6	0.000279
E127 SKIN.NHEK (NHEK-epidermal keratinocyte primary cells)	52	31.5	0.00032
E058 SKIN.PEN.FRSK.KER.03 (foreskin keratinocyte primary cells skin03)	55	34.2	0.000383
E069 BRN.CING.GYR (brain cingulate gyrus)	43	25.1	0.000495
E084 GI.L.INT.FET (fetal intestine large)	52	32.3	0.000574
E044 BLD.CD4.CD25.CD127M.TREGPC (primary T regulatory cells from peripheral blood)	33	17.7	0.000578
E103 GI.RECT.SM.MUS (rectal smooth muscle)	35	19.3	0.000614
E034 BLD.CD3.PPC (Primary T cells from peripheral blood)	46	27.7	0.000625
E010 ESDR.H9.NEUR (H9 derived neuron cultured cells)	47	28.6	0.000683
E037 BLD.CD4.MPC (Primary T helper memory cells from peripheral blood 2)	44	26.3	0.000728
E039 BLD.CD4.CD25M.CD45RA.NPC (Primary T helper naive cells from peripheral blood)	40	23.3	0.000735
E029 BLD.CD14.PC (Primary monocytes from peripheral blood)	58	37.6	0.000755
E002 ESC.WA7 (ES-WA7 cells)	24	11.7	0.000927
E071 BRN.HIPP.MID (brain hippocampus middle)	48	29.9	0.000959
E046 BLD.CD56.PC (primary natural killer cells from peripheral blood)	47	29.3	0.001073
E075 GI.CLN.MUC (colonic mucosa)	25	12.6	0.001164
E012 ESDR.CD56.ECTO (hESC derived CD56+ ectoderm cultured cells)	42	25.9	0.00165
E009 ESDR.H9.NEUR.PROG (H9 derived neuronal progenitor cultured cells)	40	24.3	0.001692
E015 ESC.HUES6 (HUES6 cells)	46	29.3	0.001842
E124 BLD.CD14.MONO (monocytes-CD14+ RO01746 primary cells)	39	24.1	0.002452
E036 BLD.CD34.CC (primary hematopoietic stem cells short term culture)	51	34.2	0.003038
E003 ESC.H1 (H1 cells)	43	27.6	0.003075
E101 GI.RECT.MUC.29 (rectal mucosa donor 29)	28	15.9	0.00323
E112 THYM (thymus)	28	16	0.003489
E042 BLD.CD4.CD25M.IL17P.PL.TPC (primary T helper 17 cells PMA-I stimulated)	38	24.2	0.004497
E014 ESC.HUES48 (HUES48 cells)	40	25.9	0.004823
E048 BLD.CD8.MPC (primary T CD8+ memory cells from peripheral blood)	34	21.6	0.007013
E007 ESDR.H1.NEUR.PROG (H1 derived neuronal progenitor cultured cells)	31	19.3	0.00718
E106 GI.CLN.SIG (sigmoid colon)	32	20.1	0.007289
E072 BRN.INF.TMP (brain inferior temporal lobe)	35	22.8	0.008904
E032 BLD.CD19.PPC (primary B cells from peripheral blood)	47	33	0.009892
E062 BLD.PER.MONUC.PC (primary mononuclear cells from peripheral blood)	22	13	0.012835
E070 BRN.GRM.MTRX (brain germinal matrix)	29	18.6	0.013775
E033 BLD.CD3.CPC (primary T cells from cord blood)	30	19.6	0.015527
E011 ESDR.CD184.ENDO (hESC Derived CD184+ endoderm cultured cells)	41	29.3	0.020598
E035 BLD.CD34.PC (primary hematopoietic stem cells)	37	26.3	0.025261
E018 IPSC.15b (iPS-15b cells)	38	27.3	0.027103
E045 BLD.CD4.CD25I.CD127.TMEMPC (primary T cells effector/memory enriched from peripheral blood)	21	13.5	0.033833
E008 ESC.H9 (H9 cells)	23	15.2	0.034026
E082 BRN.FET.F (fetal brain female)	21	13.8	0.040338
E019 IPSC.18 (iPS-18 cells)	36	26.6	0.042991

Significant enrichment for histone modifications associated with enhancer activity occurs at subsets of the growth associated SNPs in a wide range of tissue types. The diversity of tissues that showed significant enrichment reflects the biological processes that interplay to result in a heterogeneous phenotype such as stature.

GWAS3D identified 68 of the 722 SNPs as falling within regions with significant evidence of spatial connections (Supplementary Material, Table S3). Of these SNPs, we observed a growth locus (rs12153391 [lead SNP, P = 4 × 10 ^− 12] and rs26407 [spatial lead SNP, ∼1 kb from lead GWAS SNP]; (9)) located ∼85 kb upstream of FBXW11 (Fig. 3a) with strong evidence of enhancer and DNase hypersensitivity, as well as potential spatial connections. This locus (i.e. 5q35.1) showed at least partial evidence for cis-spatial connections as well as connections with distant loci (2q21.3, 5q13.2, 6p21.1, 14q31.1, 19p13.2 and Xq21.1). However, only the spatial connections between 5q35.1 and the (1) 6p21.1 or (2) 5q13.2 loci had a GWAS3D score > 10 (significant spatial evidence; Supplementary Material, Table S4).

Figure 3.

Spatial connections and trans-eQTL with 5q35.1 form the core of a regulatory network for human growth. (a) Human growth-associated SNP rs12153391, upstream of the closer C5orf50 and FBXW11 genes, instead shows regulatory eQTL evidence with the more distant upstream KCNMB1. (b) The rs12153391 locus is linked to growth disorders and height through GWAS results (hashed red circles), biological pathways (tan fill) or both through spatial associations (dashed gray arrows). This hub represents regulatory and structural associations with the GWAS loci as highlighted by significant eQTLs (green and red arrows for up- and down-regulated, respectively). Arrows pointing directly to genes represent a direct interaction, while arrows pointing to locus markers (e.g. the FBXW11-6p21.1 arrow) highlight spatial interactions within that locus without a specific gene target identified. (c) A DNA region harboring the ‘C’ version of rs12153391 acts as an enhancer. PCR-amplified genomic DNA with the indicated SNP variants were assayed for their ability to drive luciferase expression in HeLa cells. The ‘C’ but not ‘A’ allele of rs12153391 gave greater than 4-fold luminescence, indicating that competence for transcription depends on the alellic version of rs12153391. Error bars represent ±SEM from five biological replicates and significance was determined by one-way analysis of variance (ANOVA) (****P < 0.0001).

The rs12153391 locus was associated with cis-eQTLs within 5q35.1. The spatial connection was associated with the co-regulation of the KCNMB1 gene (P = 8.70 × 10 ^− 5, pituitary tissue) within the 5q35.1 locus (Fig. 3a and b, Table 2). However, despite the fact that FBXW11 is the nearest coding gene to rs12153391, this SNP does not significantly alter the expression of FBXW11 itself in any GTEx tissue type (P value: 0.26–0.95, Supplementary Material, Fig. S1). The standard value of our chosen significance cut-off (P < 1 × 10⁻⁴, see Materials and Methods) was conservative given that the minimum P value observed for the empirically characterized spatial relationship between rs9930506 and IRX3 (Fig. 1a) is 0.11 (P value range 0.11–0.94 and rs9930506-FTO P value range 0.00029–0.85; this analysis does not include cerebellum tissue, which is currently underpowered in GTEx and was the tissue in which the eQTL was identified). Our observations are not surprising, as few of the regulatory linkages to adjacent/nearby genes proposed for non-coding SNPs actually show significant (P < 1x10 ^− 4) evidence of an eQTL (Supplementary Material, Fig, S2) consistent with previous findings (25).

Table 2.

Summary of the trans-eQTL that support the FBXW11 locus

SNP	SNP locus	eQTL gene	Effect size	P value	GTEx tissue	Putative growth phenotype
5q35.1 Cis-eQTLs spatial connections
rs12153391	FBXW11	KCNMB1	−0.6	8.70E-05	Pituitary	Unknown
6p21.1 Trans-eQTLs spatial connections
rs12153391	FBXW11	TBCC	−0.12	8.40E-04	Whole Blood	Unknown
rs12153391	FBXW11	NFYA	0.21	1.10E-03	Esophagus: muscularis	Unknown
rs12153391	FBXW11	CAPN11	−0.3	1.90E-03	Artery: aorta	Unknown
rs12153391	FBXW11	DLK2	−0.26	2.50E-03	Thyroid	Unknown
rs12153391	FBXW11	C6orf226	0.32	3.10E-03	Brain: hypothalamus	Unknown
rs12153391	FBXW11	PPP2R5D	−0.15	3.30E-03	Artery: aorta	Overgrowth
rs12153391	FBXW11	NFKBIE	0.39	4.70E-03	Brain: hypothalamus	Unknown
rs26407	FBXW11	YIPF3	0.26	2.70E-04	Esophagus: muscularis	Unknown
rs26407	FBXW11	TAF8	−0.12	9.30E-04	Thyroid	Unknown
rs26407	FBXW11	NFYA	0.19	1.90E-03	Esophagus: muscularis	Unknown
rs26407	FBXW11	TBCC	−0.098	2.50E-03	Whole Blood	Unknown
rs26407	FBXW11	PGC	0.46	3.50E-03	Pituitary	Unknown
FBXW11 Trans-eQTLs with gh-related genes
rs12153391	FBXW11	GRB2	0.22	3.40E-04	Muscle: skeletal	GH Pathway
rs26407	FBXW11	GRB2	0.2	2.70E-04	Muscle: skeletal	GH Pathway
rs26407	FBXW11	IGFBP3	−0.43	2.00E-03	Pituitary	GH Pathway
rs26407	FBXW11	GH1	−0.25	3.80E-03	Lung	GH Pathway
FBXW11 trans-eQTl with cullin 1
rs12153391	FBXW11	CUL1	0.11	4.30E-03	Muscle: skeletal	Known Cullin1- FBXW11 Protein Interaction
Other FBXW11 SNPs—trans-eQTLs with GH and 3M syndrome genes
rs1368380	FBXW11	GH2	0.31	1.70E-03	Pituitary	GH Pathway
rs153750	FBXW11	AIP	−0.1	4.70E-03	Artery: aorta	Overgrowth
rs153750	FBXW11	GRB2	−0.16	4.80E-03	Muscle: skeletal	GH Pathway
rs33852	FBXW11	CUL9	0.15	5.90E-04	Skin: sun exposed (lower leg)	Cullin-FBXW
rs33852	FBXW11	SSTR2	0.25	1.30E-03	Artery: aorta	Overgrowth
rs33852	FBXW11	AIP	−0.13	1.30E-03	Esophagus: mucosa	Overgrowth
rs33852	FBXW11	CCDC8	−0.13	1.50E-03	Skin: sun exposed (Lower leg)	3M Syndrome
rs4868126	FBXW11	CUL9	−0.15	1.90E-03	Artery: aorta	Cullin-FBXW
rs4868126	FBXW11	TREM1	−0.21	3.30E-03	Nerve: tibial	Co-regulated by GH
rs4868126	FBXW11	SSTR5	0.22	4.40E-03	Esophagus: mucosa	Overgrowth
rs6892884	FBXW11	MEN1	−0.1	2.60E-03	Skin: sun exposed (lower leg)	Co-regulated by GH

The FBXW11 locus, through spatial connections reinforced by mRNA expression, is associated with the regulation of human growth across multiple pathways, including those of short stature and overgrowth. Of note, although not represented in the spatial data at the current level of coverage, the intergenic locus has eQTLs with several growth hormone (GH) pathway genes.

The rs12153391 locus (5q35.1) was captured in spatial proximity with numerous loci in trans but had particularly strong evidence for a spatial association with 6p21.1 (Fig. 3b, Supplementary Material, Table S4). Therefore, we tested for the presence of significant eQTLs between rs12153391 (the lead SNP) and 6p21.1 (Table 2). As above, the choice of significance relies on the results we obtained in testing eQTLs for the empirically characterized spatial relationship between rs9930506 and IRX3. For trans-interactions we selected P < 5 × 10 ^− 3 for significance based on previous studies of cis- (genes < 1 Mb distance from SNP, P < 1 × 10 ^− 4) and trans-eQTLs (longer distance or inter-chromosomal, P < 1 × 10 ^− 3) (25,26). rs12153391 shows significant evidence of eQTLs (P < 5 × 10 ^− 3, Table 2) with the TBCC, NFYA, CAPN11, DLK2, C6orf226, PPP2R5D and NFKBIE genes located within 6p21.1. Consistent with this, PPP2R5D has known roles in regulating growth and has been implicated in overgrowth phenotypes in family-based genetic associations (6).

Our results led us to ask whether the region harboring rs12153391 acts as an enhancer and whether variation at this SNP has potential to influence enhancer activity. Luciferase assays in Hela cells revealed a pronounced allelic effect of rs12153391 on enhancer activity, where the major allele (C) facilitates enhancer activity, while the minor allele (A) did not (Fig. 3c). Therefore, the enhancer activity associated with a region harboring rs12153391 is sensitive to the identity of the base at the SNP position.

Other clusters with significant spatial and eQTL overlaps were identified by our tests for spatial associations with growth associated SNPs (Fig. 4, Table 3, Supplementary Material, Table S3). For example, rs6061231 (GWAS lead SNP, intergenic near RPS21 and CABLES2; (9)) had spatial and trans-eQTL connections with 19q13.32. Thus, both CUL7 (6p21.1) and CCDC8 (19q13.32), which play major roles in 3M syndrome (short stature (27)), are spatially associated with intergenic growth SNPs that incorporate them into networks involved in the modulation of stature. However, the evidence for the spatial connections with both CUL7 and CCDC8 is not supported by a significant eQTL. These findings raise several interesting possibilities. First, cis or trans-eQTLs may exist but averaging across cell types within a particular tissue could hide these functional connections. Second, it is well understood that the number of spatial interactions captured by high-throughput conformation capture (HiC) is a function of sequencing depth and number of restriction sites (28). Third, there may be multiple disconnected regulatory networks working together underlying the growth phenotype. Despite these caveats, our results demonstrate that this method has significant utility for the identification of functional linkages that position intergenic SNPs into biologically relevant pathways.

Figure 4.

Summary of the trans-eQTL that support other GWAS loci. The CABLES2 locus is also linked to growth disorders and height through GWAS results (hashed red circles), biological pathways (tan fill) or both. These loci each show evidence of being linked to these changes via spatial associations (dashed gray arrows) into hubs that are related to the FBXW11 hub through relationships with common growth pathways. These hubs represent regulatory and structural associations with the GWAS loci as highlighted by significant eQTLs (green and red arrows for up- and down-regulated, respectively). Arrows pointing directly to genes represent a direct interaction, while arrows pointing to locus markers highlight spatial interactions within that locus without a specific gene target identified. These results reinforce evidence that these loci can influence the regulation of human growth across multiple pathways, including those of short stature and overgrowth.

Table 3.

Trans-eQTL associated with CABLES2 (20q13.33) and 8q21.11

SNP	SNP locus	eQTL gene	Effect Size	P value	GTEx tissue	Putative growth phenotype
20q13.33-19q13.32 Trans-eQTLs spatial connections
rs1570027	CABLES2	HIF3A	−0.23	5.40E-05	Whole blood	Unknown
rs6061231	CABLES2	PNMAL1	−0.21	1.20E-04	Brain: hypothalamus	Unknown
rs1570027	CABLES2	CBLC	0.24	2.90E-04	Thyroid	Unknown
rs6061231	CABLES2	HIF3A	−0.19	3.40E-04	Whole blood	Unknown
rs6061231	CABLES2	CBLC	0.19	2.00E-03	Thyroid	Unknown
rs1570027	CABLES2	PNMAL1	−0.19	2.60E-03	Brain: hypothalamus	Unknown
rs6061231	CABLES2	GEMIN7	−0.26	3.00E-03	Esophagus: muscularis	Unknown
rs1570027	CABLES2	IGFL2	−0.56	3.20E-03	Brain: hypothalamus	Unknown
rs1570027	CABLES2	CKM	0.52	3.80E-03	Brain: hypothalamus	Unknown
rs6061231	CABLES2	PPM1N	0.19	3.90E-03	Esophagus: mucosa	Unknown
rs6061231	CABLES2	PVRL2	−0.19	4.30E-03	Pituitary	Unknown
8q21.11-11q23.3 Trans-eQTLs spatial connections
rs4735677	8q21.11 Gene desert	POU2F3	−0.21	5.10E-06	Lung	Unknown
rs4735767	8q21.11 Gene desert	POU2F3	−0.21	8.00E-06	Lung	Unknown
rs7846385	8q21.11 Gene desert	POU2F3	−0.2	1.60E-05	Lung	Unknown
rs4735677	8q21.11 Gene desert	TMEM136	0.48	1.30E-03	Pituitary	Unknown
rs7846385	8q21.11 Gene desert	TMEM136	0.48	1.30E-03	Pituitary	Unknown
rs4735767	8q21.11 Gene desert	TMEM136	0.48	1.50E-03	Pituitary	Unknown
rs4735767	8q21.11 Gene desert	UBE4A	0.11	2.90E-03	Lung	Unknown
rs4735767	8q21.11 Gene desert	TRAPPC4	0.21	3.90E-03	Adipose: subcutaneous	Unknown
rs4735767	8q21.11 Gene desert	PVRL1	−0.18	3.90E-03	Esophagus: mucosa	Unknown
rs4735767	8q21.11 Gene desert	CCDC84	0.19	4.60E-03	Stomach	Unknown
rs4735677	8q21.11 Gene desert	TREH	0.29	4.70E-03	Heart: left ventricle	Height GWAS
rs7846385	8q21.11 Gene desert	PVRL1	−0.17	4.80E-03	Esophagus: mucosa	Unknown
rs7846385	8q21.11 Gene desert	UBE4A	0.1	5.00E-03	Lung	Unknown
8q21.11 Gene desert trans-eQTLs spatial connections with stature
rs4735767	8q21.11 Gene desert	GH1	−0.38	7.60E-04	Esophagus: mucosa	Growth hormone pathway

The CABLES2 and 8q21.11 loci have numerous spatial connections with functional (eQTL) support. These results are consistent with a role for CABLES2 and 8q21.11 in the regulation of human growth across multiple pathways. Of note, although not represented in the spatial data at the current level of coverage, the 8q21.11 intergenic locus has an eQTL with growth hormone (GH1).

Collectively, our findings are consistent with the standard GWAS hypothesis that polymorphisms in common single nucleotides (i.e. rs12153391 and rs6061231) have low-effect impacts across multiple pathways (below significance), the sum of which results in measurable, significant SNP-phenotype associations (29). This study serves to illustrate a methodology that enables a better understanding of the impact of polymorphisms on disease and potentially explains intergenic SNP-phenotype associations.

Discussion

By combining spatial connections and RNA expression data, we have assigned function to intergenic regions of DNA that, despite strongly significant phenotypic associations, were not clearly functionally defined by the local genetic landscape. Collectively, our results support the conclusion that rs12153391 (5q35.1) and 6p21.1 form a core component of a gene regulatory hub that underlies growth phenotypes. Our findings suggest that intergenic variants near FBXW11 could be key modifiers of growth and should be further studied to better understand the genetic regulation of growth, with a focus on the diagnosis and treatment of growth disorders. We identified a number of other relationships for significant SNPs but concentrated on rs12153391 because of its potential to be a regulatory hub at the center of the growth phenotype.

Combined, our results show that rs12153391 is associated with altered transcription factor binding sites and significant long-range interactions (supported by eQTLs), implicating this region as a long-range regulator of transcription (i.e. an enhancer). rs12153391 alters a DNA binding motif that is recognized by CCNT2, Maf, SP1 and Zfp281 (Haploreg v4.1 (30)). Moreover, rs12153391 is in linkage disequilibrium (LD) with rs67129974 (r²=0.95), which falls within a NFYA and NFYB binding domain, as well as affecting DNA motifs recognized by Duxl, NF-Y, Nanog and YY1 (Haploreg v4.1 (30)). In combination with the finding of the rs12153391-NFYA eQTL, the fact that rs67129974 alters a NFYA protein binding domain implicates the trans-interaction in self-regulation. The potential effect of this connection on growth is further emphasized by the functional association of NFYA and SP1 in the transcriptional regulation of the ARNTL gene, which contains growth associated SNPs (Supplementary Material, Table S2 (9)). Finally, rs12153391 and rs67129974 alter DNA binding motifs that are recognized by Zfp281 and Nanog, respectively, two repressors that act together in embryonic stem cell differentiation (31). Collectively, our results are consistent with the presence of histone modifications associated with enhancer activity surrounding rs12153391 in seven different tissue types, including fetal lung and muscle (Supplementary Material, Table S2, Haploreg v4.1 (30)).

The spatial information we incorporate into our analysis provides a mechanism whereby we can increase the information content of GWAS studies without resorting to simply trying to increase the sample sizes. This is particularly relevant as Wood et al. (9) noted that SNP-gene attribution on 10⁵ samples had a diminishing marginal effect leading to the identification of variants that represented new members of previously described growth-associated gene sets. In contrast, the study of the three-dimensional organization of genomes identifies spatial associations that contribute to gene regulation, providing an alternative explanation for the source of pathogenicity of the majority of the genetic variants. Notably, out of the 722 SNPs in our analysis, we identified 22 GWAS SNPs as having CCCTC-binding factor (CTCF) binding sites and 10 more as possibly altering CTCF binding motifs (Supplementary Material, Table S2, Haploreg v4.1) (32). Of those, rs6061231 (CABLES2, 20q13.33), rs1380294 (LCORL, 4p15.31) and rs960273 (GNA12, 7p22.2) also showed evidence of spatial regulation (Table 3, Supplementary Material, Table S3). CTCF is recognized as a major contributor to the 3D structure of chromatin (33). Therefore, as higher coverage spatial data becomes available, incorporating 3D structural aspects into GWAS analyses will impact our understanding of many of the most significant SNP-gene relationships in growth. We contend that, all things being equal, similar success rates can be expected for assigning SNPs into regulatory networks that underlie other complex phenotypes.

There are two major caveats to this study. First, the spatial and expression data used to define these networks was limited to relatively low coverage. Specifically, proximity-ligation techniques capture trans-interactions between chromosomes at a reduced level that reflects the balance between the intra-chromosomal folding and the mixing that occurs between chromosome territories. Similarly, the expression (eQTL) data are powered to detect higher-effect cis-acting RNA expression changes. Additionally, the eQTL data only capture a snapshot of expression changes in homogenous cells of a specific lineage and as such does not represent all cells and time-points. Thus, it is likely that our results represent only a sub-network within the complete regulatory network that acts within cells associated with human growth. Therefore, future work will be required to expand and clarify these networks by increasing the resolution and cell-type specificity (34) of the spatial networks (e.g. (35)) and expression data (25).

The second major caveat to this study is an extension of the first; namely the choice of P value at which eQTLs are determined to be statistically significant affects the numbers of eQTLs that are detected. However, the P value that should be used is the subject of some debate for both cis- and trans-interactions (21,25,26). A global eQTL analysis typically requires millions of tests for cis-connections, while a trans-eQTL global analysis requires billions of tests. Thus, minimum P values of 10⁻³–10⁻⁴ and 10⁻⁸–10⁻⁹ are required to validate global cis and trans eQTLs at an FDR < 0.05, respectively (26). Therefore, we took a practical approach to limit the impact of the multiple testing corrections by only checking eQTLs for loci that spatially interact. The biological significance of our findings was supported by our analysis of the minimum P values attained for the well-characterized functionally relevant FTO-IRX3 example of long distance cis-regulation. Thus, it is clear that many of the trans-eQTLs that fall below the level typically used to indicate global significance (10⁻⁸) are biologically significant.

The interconnections between biological pathways are complex and multi-faceted. We have shown that analyses that incorporate spatial data and eQTLs are capable of elucidating the functional relationship between statistically significant inter-genic SNPs and the phenotype of interest. Integrating these diverse (but complementary) datasets provides a critical next step that helps to better explain the biological pathways underlying GWAS results, while limiting the effort that is expended chasing false positives.

Materials and Methods

Identification of GWAS variants

We identified adult height associated SNPs from the GWAS Catalog for SNPs (Genome3D Pipeline; Fig. 2). We also identified growth-associated SNPs from three GWAS meta-analyses for infant (7), pubertal (8) and adult height (9). We confirmed height and/or development associations for each selected locus by searching the published literature (PubMed), the online Mendelian Inheritance in man database and Uniprot for functional classification. In total, 722 significant SNPs, identified in the meta-analyses, were used in the downstream analyses. SNPs within the same haplotype block were combined and named according to their locus. Haplotype blocks were defined as all SNPs sharing LD coefficients of at least 0.8 r² and 0.8 D`.

Spatial and functional analysis

To identify how human variation contributes to growth, we used HiC spatial genomic connectivity (HiC) and functional (eQTL) data to determine how the non-genic GWAS loci contribute to gene regulatory networks. The web-based program GWAS3D (32) was used to identify physical connections (as captured by proximity ligation (19,20) and ENCODE DNase hypersensitivity) between the most significant GWAS meta-analysis SNPs. Cell-type-specific eQTLs were tested using the Genotype-Tissue Expression project database (25).

GTex eQTL data are currently powered to detect cis-eQTLs and not trans-eQTLs, which usually lie far away (>1 Mb) on the same chromosome or on a different chromosome (25). Therefore, we limited false positives in our trans-eQTL results by only testing eQTLs supported by evidence of SNP-gene spatial interactions. Significance thresholds for eQTL were chosen based on prior literature (25,26): cis-eQTL (genes <1 Mb distance from SNP, P < 1 × 10 ^− 4) and trans-eQTL (longer distance or inter-chromosomal, P < 1 × 10 ^− 3).

The GTEx database contains data on many broad categories of heterogeneous tissue types. Since the eQTL data are only being used to further reinforce evidence of enhancer activity, the ability to detect eQTLs across different tissue types provides insight into tissue-specific roles for the connection. Thus, a locus with both positive and negative eQTL effects in differing tissues is potentially showing evidence of cell-specific activity across diverse processes in space and time (e.g. growth and development).

It has previously been shown that human height variants from the GWAS catalog were strongly associated with ‘Blood Vessel’ (three GTEx categories: artery aorta, artery coronary and artery tibial) and weakly associated with nerve (GTEx: nerve tibial), breast (GTEx: Breast Mammary Tissue) and fallopian tube (GTEx: fallopian tube) (33). This is surprising, as short and very tall stature are more often associated with defects in regulation within the pituitary or hypothalamus. Therefore, we tested all of the major tissue types available in GTEx: subcutaneous adipose, aortic artery, tibial artery, heart (left ventricle), tibial nerve, sun-exposed skin (lower leg), skeletal muscle, mucosa and muscularis of the esophagus, thyroid and whole blood. We also tested pituitary tissue in specific cases, where appropriate, but the GTEx sample size was underpowered to provide highly significant results for Pituitary.

Luciferase assays

For determining enhancer activity of the 5q35.1 locus, the regions spanning rs12153391 and rs26407 were polymerase chain reaction (PCR) amplified from genomic DNA, cloned into the Gateway-adapted pGL4.23-GW (Addgene Plasmid no. 60323) (36) and sequenced to confirm genotype. Primers for PCR amplification are listed in Supplementary Material, Table S5. HeLa cells were seeded at 5 × 10³ cells per well in a 96-well plate, grown in DMEM media supplemented with 10% fetal bovine serum (Thermo Scientific, no. 11995-065) 1 day prior to transfection with Lipofectamine 3000 (Thermo Scientific, no. L3000008). Forty-eight hours after transfection, luciferase assays were performed using the Promega Dual Glo Luciferase Assay System (no. PME2920). Luminescence was normalized to Renilla, and results are from five independent biological replicates.

Supplementary Material

Supplementary Material is available at HMG online.

Conflict of Interest Statement. None declared.

Funding

This work was supported by an Auckland University Future Research Development Fund grant (3702119 to J.M.O.S.), a University of Auckland Scholarship (to W.S.), a Royal Society of New Zealand Marsden grant (11-UOO-027, to J.A.H.) and an HRC grant (15-623, to J.M.O.S. and J.A.H.).