Rare and low-frequency coding variants alter human adult height

Abstract

Height is a highly heritable, classic polygenic trait with approximately 700 common associated variants identified through genome-wide association studies so far. Here, we report 83 height-associated coding variants with lower minor-allele frequencies (in the range of 0.1-4.8%) and effects of up to 2 centimetres per allele (such as those in IHH, STC2, AR and CRISPLD2), greater than ten times the average effect of common variants. In functional follow-up studies, rare height-increasing alleles of STC2 (giving an increase of 1-2 centimetres per allele) compromised proteolytic inhibition of PAPP-A and increased cleavage of IGFBP-4 in vitro, resulting in higher bioavailability of insulin-like growth factors. These 83 height-associated variants overlap genes that are mutated in monogenic growth disorders and highlight new biological candidates (such as ADAMTS3, IL11RA and NOX4) and pathways (such as proteoglycan and glycosaminoglycan synthesis) involved in growth. Our results demonstrate that sufficiently large sample sizes can uncover rare and low-frequency variants of moderate-to-large effect associated with polygenic human phenotypes, and that these variants implicate relevant genes and pathways.

Extended Data Figure 1

Flowchart of the GIANT ExomeChip height study design.

Extended Data Figure 2

Height ExomeChip association results. (A) Quantile-quantile plot of ExomeChip variants and their association to adult height under an additive genetic model in individuals of European ancestry. We stratified results based on allele frequency. (B) Manhattan plot of all ExomeChip variants and their association to adult height under an additive genetic model in individuals of European ancestry with a focus on the 553 independent SNPs, of which 469 have MAF>5% (grey), 55 have MAF between 1 and 5% (green), and 29 have MAF<1% (blue). (C) Linkage disequilibrium (LD) score regression analysis for the height association results in European-ancestry studies. In the plot, each point represents an LD Score quantile, where the x-axis of the point is the mean LD Score of variants in that quantile and the y-axis is the mean χ² statistic of variants in that quantile. The LD Score regression slope of the black line is calculated based on Equation 1 in Bulik-Sullivan et al. which is estimated upwards due to the small number of common variants (N=15,848) and the design of the ExomeChip. The LD score regression intercept is 1.4, the λ_GC is 2.7, the mean χ² is 7.0, and the ratio statistic of (intercept -1) / (mean χ² -1) is 0.067 (standard error=0.012). (D) Scatter plot comparison of the effect sizes for all variants that reached significance in the European-ancestry discovery results (N=381,625) and results including only studies with sample sizes >5000 individuals (N=241,453).

Extended Data Figure 3

Height ExomeChip association results in African-ancestry populations. Among the all-ancestry results, we found eight variants for which the genetic association with height is mostly driven by individuals of African ancestry. The minor allele frequency of these variants is <1% (or monomorphic) in all ancestries except African-ancestry individuals. In individuals of African ancestry, the variants had allele frequencies between 9 and 40%.

Extended Data Figure 4

Concordance between direct conditional effect sizes using UK Biobank (x-axis) and conditional analysis performed using a combination of imputation-based methodology and approximate conditional analysis (SSimp, y-axis). The Pearson's correlation coefficient is r=0.85. The dashed line indicates the identity line. The 95% confidence interval is indicated in both directions. Red, SNPs with P_cond>0.05 in the UK Biobank; Green, SNPs with P_cond≤0.05 in the UK Biobank.

Extended Data Figure 5

Heritability estimated for all known height variants in the first release of the UK Biobank dataset. (A) We observed a weak but significant positive trend between minor allele frequency (MAF) and heritability explained (P=0.012). (B) Average heritability explained per variant when stratifying the analyses by allele frequency or genomic annotation. For heritability estimations in UKBB, variants were pruned to r²< 0.2 in the 1000 Genomes Project data set, and the heritability figures are based on h²=80% for height.

Extended Data Figure 6

Comparison of DEPICT gene set enrichment results based on coding variation from ExomeChip (EC) or non-coding variation from genome-wide association study data (GWAS). The x-axis indicates the P-value for enrichment of a given gene set using DEPICT adapted for EC data, where the input to DEPICT is the genes implicated by coding EC variants that are independent of known GWAS signals. The y-axis indicates the P-value for gene set enrichment using DEPICT, using as input the GWAS loci that do not overlap the coding signals. Each point represents a meta-gene set, and the best P-value for any gene set within the meta-gene set is shown. Only significant (false discovery rate < 0.01) gene set enrichment results are plotted. Colors correspond to whether the meta-gene set was significant for EC only (blue), GWAS only (green), both but more significant for EC (purple), or both but more significant for GWAS (orange), and the most significant gene sets within each category are labeled. A line is drawn at x = y for ease of comparison.

Extended Data Figure 7

Heat map showing entire DEPICT gene set enrichment results (analogous to Fig. 2 in the main text). For any given square, the color indicates how strongly the corresponding gene (shown on the x-axis) is predicted to belong to the reconstituted gene set (y-axis). This value is based on the gene's Z-score for gene set inclusion in DEPICT's reconstituted gene sets, where red indicates a higher Z-score and blue indicates a lower one. The proteoglycan binding pathway was uniquely implicated by coding variants (as opposed to common variants) by both DEPICT and the Pascal method. To visually reduce redundancy and increase clarity, we chose one representative “meta-gene set” for each group of highly correlated gene sets based on affinity propagation clustering (see Methods and Supplementary Information). Heat map intensity and DEPICT p-values correspond to the most significantly enriched gene set within the meta-gene set; meta-gene sets are listed with their database source. Annotations for the genes indicate whether the gene has OMIM annotation as underlying a disorder of skeletal growth (black and grey) and the minor allele frequency of the significant EC variant (shades of blue; if multiple variants, the lowest-frequency variant was kept). Annotations for the gene sets indicate if the gene set was also found significant for EC by the Pascal method (yellow and grey) and if the gene set was found significant by DEPICT for EC only or for both EC and GWAS (purple and green). Abbreviations: GO: Gene Ontology; KEGG: Kyoto encyclopedia of genes and genomes; MP: mouse phenotype in the Mouse Genetics Initiative; PPI: protein-protein interaction in the InWeb database.

Extended Data Figure 8

Heatmaps showing associations of the height variants to other complex traits; −log10(P-values) are oriented with beta effect direction for the alternate allele, white are missing values, yellow are non-significant (P>0.05), green to blue shading for hits with positive beta in the other trait and P-values between 0.05 and <2×10^-7 and, orange to red shading for hits with negative beta in the other trait and P-values between 0.05 to <2×10^-7. Short and tall labels are given for the minor alleles. Clustering is done by the complete linkage method with Euclidean distance measure for the loci. Clusters highlight SNPs that are more significantly associated with the same set of traits. (A) Variants for which the minor allele is the height-decreasing allele. (B) Variants for which the minor allele is the height-increasing allele.

Figure 1

Variants with a larger effect size on height variation tend to be rarer. We observed an inverse relationship between the effect size (from the combined “discovery+validation” analysis, in cm on the y-axis) and the minor allele frequency (MAF) for the height variants (x-axis, from 0 to 50%). We included in this figure the 606 height variants with P<2×10^-7.

Figure 2

Heat map showing subset of DEPICT gene set enrichment results. The full heat map is available as Extended Data Fig. 7. For any given square, the color indicates how strongly the corresponding gene (shown on the x-axis) is predicted to belong to the reconstituted gene set (y-axis). This value is based on the gene's Z-score for gene set inclusion in DEPICT's reconstituted gene sets, where red indicates a higher Z-score and blue indicates a lower one. The proteogly can binding pathway (bold) was uniquely implicated by coding variants by DEPICT and PASCAL. To visually reduce redundancy and increase clarity, we chose one representative “meta-gene set” for each group of highly correlated gene sets based on affinity propagation clustering (Supplementary Information). Heat map intensity and DEPICT P-values correspond to the most significantly enriched gene set within the meta-gene set; meta-gene sets are listed with their database source. Annotations for the genes indicate whether the gene has OMIM annotation as underlying a disorder of skeletal growth (black and grey) and the minor allele frequency of the significant ExomeChip (EC) variant (shades of blue; if multiple variants, the lowest-frequency variant was kept). Annotations for the gene sets indicate if the gene set was also found significant for EC by PASCAL (yellow and grey) and if the gene set was found significant by DEPICT for EC only or for both EC and GWAS (purple and green). Abbreviations: GO: Gene Ontology; MP: mouse phenotype in the Mouse Genetics Initiative; PPI: protein-protein interaction in the In Web database.

Figure 3

STC2 mutants p.Arg44Leu (R44L) and p.Met86Ile (M86I) show compromised proteolytic inhibition of PAPP-A. (A) Schematic representation of the role of STC2 in IGF-1 signaling. Partial inactivation of STC2 by height-associated DNA sequence variation could increase bioactive IGF-1 through reduced inhibition of PAPP-A. (B) Western blot analysis of recombinant STC2 wild-type and variants R44L and M86I. (C) Covalent complex formation between PAPP-A and STC2 wild-type or variants R44L and M86I. Separately synthesized proteins were analyzed by PAPP-A Western blotting following incubation for 8 h. In the absence of STC2 (Mock lane), PAPP-A appears as a single 400 kDa band (*). Following incubation with wild-type STC2, the majority of PAPP-A is present as the approximately 500 kDa covalent PAPP-A:STC2 complex (#), in which PAPP-A is devoid of proteolytic activity towards IGFBP-4. Under similar conditions, incubation with variants R44L or M86I appeared to cause less covalent complex formation with PAPP-A. The gels are representative of at least three independent experiments. (D) PAPP-A proteolytic cleavage of IGFBP-4 following incubation with wild-type STC2 or variants for 1-24 h. Wild-type STC2 causes reduction in PAPP-A activity, with complete inhibition of activity following 24 h incubation. Both STC2 variants show increased IGFBP-4 cleavage (i.e. less inhibition) for all time points analyzed. Mean and standard deviations of three independent experiments are shown. One-way repeated measures analysis of variance followed by Dunnett's post-test showed significant differences between STC2 wild-type and variants R44L (P<0.001) and M86I (P<0.01).