Height is a complex trait influenced by various nutritional, environmental, and genetic factors. Previous studies have revealed the important role of genetic factors in determining human stature. The heritability of height is reported to be >75%.^1–3 Linkage and association analyses have identified a number of genomic regions and candidate genes underlying human stature variation.^4–11

In recent years, the vitamin D receptor (VDR) has been studied for its relationship with human height. Linkage studies have found that the genomic region of 12p11.2-q14 harbouring the VDR gene is linked to human height variation.^4,7 Association of VDR polymorphisms with height has also been reported.^8,9,12–15 These observations are expected, given the importance of the vitamin D endocrine system to skeletal growth^8,9,16 and to a number of other bone phenotypes closely correlated with height.^9,17–21 Most of the previous association studies adopted random unrelated population samples. Statistically, such a method is more powerful than linkage or family based association methods, and the samples required are much easier to collect. However, given potential population stratification, population associations could lead to spurious results.²² In comparison, family based approaches, such as the transmission disequilibrium test (TDT), are immune to population stratification and able to test linkage and association simultaneously.

In the present study, using the quantitative transmission disequilibrium test (QTDT) method, we examined four highly polymorphic single nucleotide polymorphisms (SNPs) in the VDR gene for linkage and/or association with height in a large sample of 1873 white subjects from 406 nuclear families. Our results showed the importance of the VDR gene, especially the polymorphisms near its 3′ UTR for height variation in whites. Further tests using stratified unrelated samples randomly chosen from the parents of these 406 nuclear families suggested that such an association might be present in females only and be influenced by menstrual status. New findings about another two marker loci, −4817A/G and FokI C/T, are also reported and discussed in this study.

SUBJECTS AND METHODS

Study population

All the subjects were recruited at the Osteoporosis Research Center of Creighton University for the purpose of searching for genes underlying osteoporosis and obesity. The study was approved by the institutional review board of Creighton University. All the study subjects were white, of European origin, and gave their informed consent before entering the project. The sampling scheme and exclusion criteria, which aimed at reducing potential confounding effects on study phenotypes, have been detailed previously by Deng et al.²³ The white nuclear families were recruited randomly in terms of height variation.

A total of 406 nuclear families with 1873 individuals were analysed in this study, including 740 parents, 389 male children, and 744 female children. Of these, 342 families were composed of both parents and at least one offspring; the number of offspring ranged from 3 to 12, with a mean (SD) of 4.62 (1.78) individuals. The remaining 64 families, with one or no parent, contained two or more children. The percentage of nuclear families with one, two, three and more than three children were 27.2%, 22.7%, 22.7%, and 27.4%, respectively, yielding 1512 sibling pairs (amenable to linkage testing) in our sample.

For population studies, three groups of unrelated subjects were chosen from the 740 parents: (a) men, (b) pre-menopausal and peri-menopausal women, and (c) post-menopausal women.

Subjects were asked for their age by questionnaire. All height measurements without shoes were made using a standard wall mounted statiometer.

SNP genotyping

SNPs were selected from previously published studies and the public database such as dbSNP (http://www.ncbi.nlm.nih.gov/SNP/). Four SNPs were chosen within the VDR gene based on a comprehensive consideration of the criteria: (a) functional relevance and importance; (b) degree of heterozygosity (that is, minor allele frequencies); (c) relative position in the gene; (d) their use in previous genetic epidemiology studies. They were rs2238136 (−4817A/G) at intron 1, rs2228570 (FokI C/T) at exon 2 start codon, rs1544410 (BsmI A/G) at intron 8, and rs731236 (TaqI T/C) at exon 9. The latter three SNPs have been commonly studied previously, and for convenience are referred to as FokI, BsmI, and TaqI, respectively, throughout this paper. The information for the selected markers and their PCR primers are presented in table 1, and their positions in the VDR gene illustrated in fig 1.

DNA was extracted from whole blood using a commercial isolation kit (Gentra Systems, Minneapolis, MN, USA). The genotyping procedure for all SNPs was similar, involving PCR and Invader assay reaction (Third Wave Technology, Madison, WI, USA). PCR was performed in a 10 μl reaction volume with 35 ng genomic DNA, 0.2 mmol/l each dNTP, 1× PCR buffer and 1.5 mmol/l MgCl₂, 0.4 μmol/l each of the primers, and 0.35 U Taq polymerase (ABI, Applied Biosystems, Foster City, CA, USA). The following procedure was used: 95°C for 5 minutes, 30 cycles of 94°C for 1 minute, 50°C for 1 minute, 72°C for 1 minute, and then 72°C for 5 minutes. After amplification, an Invader reaction was performed in a 7.5 μl reaction volume, with 3.75 μl diluted PCR product (1:20), 1.5 μl probe mix, 1.75 μl Cleavase FRET mix, and 0.5 μl Cleavase enzyme/MgCl₂ solution (Third Wave Technology). The reaction mix was overlaid with 15 μl mineral oil and denatured at 95°C for 5 minutes, then incubated at 63°C for 20 minutes. All the PCR and Invader assay reactions were performed on an ABI 9700 thermal cycler. After the Invader reaction, the fluorescence intensity for both colours (FAM dye and Red dye) was read using a Cytofluor 4000 multiwell plate reader (ABI). The data were analysed with Invader Analyzer software (Third Wave Technology), and the genotype for every sample was classified according to the ratio of the fluorescence intensity of the two dyes. PedCheck software²⁴ was used to verify Mendelian inheritance of the alleles within each family and the alleged familial relationships, and to aid in check of genotyping errors. To correct errors of non-Mendelian inheritance for problematic individuals, several closely related individuals from the same pedigree were re-genotyped along with controls. Genotype data that did not pass PedCheck after 5–7 attempts at re-genotyping were deleted and treated as missing data. At the end of all analyses, the average rate of missing genotype data (here, it equals Mendelian segregation error rate) was ∼0.4%, and the average genotyping error rate estimated through blind duplicating was ∼0.1%. Thus, the overall error rate of genotype data in our study was ∼0.5%.

Statistical analyses

Using the genotype data of all the 1873 subjects, we used a maximum likelihood method implemented in the SOLAR program (available at http://www.sfbr.org/solar) to estimate the allele frequencies of the four testing SNPs. The estimated allele frequencies were then compared with the corresponding data in the literature and dbSNP. Genotype frequencies for each of the four SNPs among unrelated parents were also tested for their concordance with Hardy-Weinberg ratio by χ² test.

To test for linkage and association with adult height variation, the program QTDT (http://csg.sph.umich.edu/pn/index.php?furl=/abecasis/QTDT) was used. This is based on the variance component framework and can test total association, within family association, population stratification, linkage, and linkage while modelling association simultaneously. The orthogonal model of Abecasis et a.²⁵ was adopted in the analyses, where the total association was partitioned into orthogonal within and between family components (β_b and β_w, respectively). The between family component is specific for each nuclear family and could be confounded by population stratification. However, the within family component is significant only in the presence of linkage disequilibrium caused by close linkage. Thus, a QTDT test of the significance of β_w based on allelic transmission is robust to population stratification.

In addition, population stratification is tested according to whether β_b = β_w.²⁶ If there is no evidence for population stratification (β_b = β_w), the test results of total association should be more powerful.²⁵ Linkage tests are based on the standard variance component methods and identity by descent (IBD) sharing among relatives. The program SimWalk2 (http://genetics.ucla.edu/home/software/simwalk2) was used to calculate IBD probabilities. In addition, a test for linkage while modelling association, suggested by Fulker et al,²⁶ was also implemented in QTDT to evaluate whether a candidate marker could be a functional variant underlying the trait. If linkage signal is still significant after modelling association, it means the study marker is not the functional mutation.

Genehunter (version 2.1; http://www.fhcrc.org/labs/kruglyak/downloads/gh-v2.1_r5.tar.gz) was used to analyse linkage disequilibrium (LD) and reconstruct the haplotypes. Genehunter extracts complete multipoint inheritance information to infer maximum likelihood haplotypes for all individuals in nuclear families.²⁷ Pairwise LD between the SNPs was calculated by the normalised measure in the parents,²⁸D′ = D/D_max, where D = π₁₁π₂₂−π₁₂π₂₁ (π is sample estimate of haplotype frequencies); D_max = min(π₁₊π₊₂, π₊₁π₂₊) when D>0, and D_max = min(π₁₊π₊₁, π₊₂π₂₊) when D<0 (π is estimate of SNP frequencies). The statistical significance of the observed LD was examined by the Monte-Carlo approximation of Fisher’s exact test. We excluded −4817A/G and FokI loci from haplotype reconstruction and performed the QTDT using the haplotypes reconstructed by the BsmI and TaqI loci, where the only significant strong LD was found.

For the three samples of unrelated subjects in the parental generation, population based association analyses were performed using the PROC GLM section of the SAS program (version 8.0e; SAS Institute Inc. Cary, NC, USA). Normality tests for height residuals adjusted by age and sex were performed by MINITAB (Minitab Inc., State College, PA, USA).

RESULTS

Allele frequencies and haplotype reconstruction

The basic characteristics of the study subjects are presented in table 2. All the 1873 individuals were aged >19 years at the time of examination and had reached (or nearly so) their final height.^7–9,14 After adjusting for the effects of both age and sex, the residuals of height of all the subjects followed normal distribution (p>0.15). The height residuals adjusted by age in each group of unrelated subjects (men, pre-menopausal and peri-menopausal women, and post-menopausal women) followed normal distribution also (p>0.15).

The results of allele frequencies of the four SNPs are summarised in table 3. The range of minor allele frequencies was 0.28–0.42 and all of the frequencies were consistent with those of white subjects in dbSNP and references (table 3; Dvornyk et al²⁹). We estimated the genotype frequencies of these four SNPs among the parents only (data not shown) and did not find deviation from Hardy-Weinberg equilibrium except for the marker −4817 A/G (table 3).

For pairwise linkage disequilibrium (LD) among these four SNPs, significant strong LD was found only between the BsmI and TaqI markers (|D’| = 0.762, p<0.0001; fig 1, table 4). Therefore, we reconstructed four different haplotypes on the basis of these two markers. The most common two haplotypes, haplotypes 2 and 3 (AC and GT) accounted for 90.3% of the total sample (table 4). Haplotypes 1 and 4 (AT and GC) accounted for 5.2% and 4.5% of the total sample, respectively. Comparing with Uitterlinden’s resolved haplotype structure at the 3′ end of the VDR gene for whites,³⁰ our haplotype 3 (GT) corresponds to “bT”, which is the combination of their haplotype 1 (baT) and 3 (bAT) (A/a is the allele of VDRApaI locus). The frequency of the bT haplotype in our study is consistent with theirs (table 4).

Family based QTDT study in white nuclear families

The heritability of height in our sample was about 78%, which lies in the estimated heritability range from 76% to 90%.⁴ Table 5 summarises the results of the linkage and association analyses via the program QTDT for height variation in the whole group (1873 individuals from 406 white nuclear families), with age and sex included as covariates. For both single SNP and haplotype analyses, we did not find any evidence of population stratification. This finding might rule out, to some extent, the confounding effect of population stratification on the test of total association by QTDT. Single locus analyses yielded the following results. Firstly, linkage signals (based on those nuclear families with two or more children) were found for all of the four testing SNPs, suggesting linkage of the VDR gene with adult height variation. After modelling association, there was still linkage signal at the loci of all four SNPs, suggesting that they are not functional polymorphisms underlying height. Total association with height was observed for the −4817A/G, BsmI, and TaqI markers (p = 0.015, 0.012, and 0.017, respectively). Furthermore, we found evidence for within family association with height at both the BsmI and TaqI loci (p = 0.048 and 0.039).

Based on haplotypes reconstructed by BsmI and TaqI markers, multi-allelic testing detected further evidence for linkage and total association of the VDR gene with height variation (p = 0.05 and 0.001, respectively). The linkage signal was still marginally significant after modelling association (p = 0.077). All of these were consistent with results from single locus analysis. Allele-wise analysis yielded interesting and compatible results (table 5). The most significant and consistent association with stature occurred with haplotype 3 (GT). Both total and within family associations were significant for this haplotype (p = 0.0006 and 0.033). Subjects carrying haplotype 3 were, on average, 1% (1.6 cm) taller than those without it (p = 0.003).

Population based association analyses for the VDR bT haplotype

The importance of the VDR haplotype 3 (GT, or bT haplotype³⁰) to white height variation was further investigated by the population association analyses. The frequencies of the bT haplotype in the three parental groups obeyed Hardy-Weinberg equilibrium. The association between height and bT haplotype was most significant in postmenopausal women (p =  0.002, table 6). Using Cohen’s d value as the measure of effect size, the effect was medium when comparing bT homozygotes with non-carriers in postmenopausal group (d = 0.5). A trend for this association was also observed in pre-menopausal and peri-menopausal women (p = 0.05, table 6). Such findings were consistent with our QTDT data. However, we could not find an association of bT haplotype with height in men (p = 0.34, table 6).

DISCUSSION

In this study, a family based approach was used to test linkage and/or association between the VDR gene and stature. Such a method, although less powerful than population based association studies, is less subject to potential population stratification.^22,31 To our knowledge, this study is the first adopting a TDT approach to test the importance of a prominent candidate gene of stature. Population based association analyses were also performed for the VDR bT haplotype, the one most significantly associated with height in our QTDT analyses.

The sample size in this study was quite large compared with many earlier studies. This sample, comprising 1873 subjects from 406 white nuclear families, renders a relatively high statistical power. Assuming a marker is in strong LD (|D′| = 0.9) with a functional mutation that accounts for about 2% of phenotypic variation, our simulation shows that the present study has >90% power to detect association via the TDT. In addition, the minor allele frequencies of the four tested SNPs are fairly high, ensuring a sufficient number of nuclear families informative for the TDT tests. However, this power calculation was only a rough estimate.

Among the four SNPs studied, the BsmI and TaqI markers and their haplotypes could represent the genetic diversity of the 3′ UTR of the VDR gene, which is a 3.2 kb sequence immediately adjacent to exon 9.³² In addition, recent data show that BsmI and TaqI localise to the so called block B that covers exons 3–9, and can to some extent capture information about SNPs in this block.³³ It may be assumed that the short 3′ UTR is located in block B.^32,33 −4817 A/G and FokI polymorphisms are present in and near the promoter region of the VDR gene. They are not in LD with the 3′ polymorphisms in the way that BsmI and TaqI are (Uitterlinden et al³² and fig 1). More interestingly, we did not find LD between −4817 A/G and FokI, although they are positioned close to each other (∼4.8 kb apart((fig 1). These findings are consistent with the recent opinion that FokI is located in a 1.3 kb LD breaking spot separating the so called blocks B and C.³³ Further, it can be deduced that −4817 A/G is within block C.³³ Thus, they should be considered as different independent markers. Testing these two markers may partially reveal the association situation of the 5′ end of the VDR gene with height variation.

In our sample, subjects carrying the b (G base at BsmI) or T (T base at TaqI) alleles were taller than those carrying the B (A base at BsmI) or t (C base at TaqI) alleles. The carriers of the most common haplotype, bT, (haplotype 3; GT) were also taller than non-carriers. These findings are largely supported by the previous linkage and association studies that had addressed the significance of VDR to adult height variation. The unanimous linkage with stature at all of the studied loci in VDR (table 5) was consistent with the findings of the first two genomewide linkage studies of height in whites, showing that the lod scores of the genomic region centered on the VDR gene were 3.35 and 1.86, respectively.^4,7 A recent genomewide linkage scan for stature by our group also revealed the importance of the VDR harbouring region to height variation in white subjects (lod = 1.04 and 1.33, respectively, in two separate sets of pedigrees that are different from pedigrees used in this study, but sampled from the same US white population).³⁴ Population based association studies in white populations provide direct evidence of support. In an ethnically homogenous group of white children and young adults, the VDR haplotype bAT was associated with increased height.⁹ Another study showed that the BB boys were shorter at birth, had slower growth up to late puberty, and shorter adult stature than their bb counterparts.⁸ Moreover, a recent systematic review of our group based on data from a large number of association studies suggested that bb or TT genotypes were more advantageous in terms of bone metabolism, calcium homeostasis, bone accrual during childhood, and bone retention later in life.³⁵

It was generally believed that BsmI and TaqI loci are in strong LD, such that the b allele is always associated with the T allele while the B allele associated with the t allele.^8,9,30,32,35 This notion was further confirmed in our study (fig 1, table 4). The BsmI and TaqI SNPs were not likely to be the functional polymorphisms per se^36,37 (table 5); however, their LD extends into the 3′ UTR, which harbours a possible functional relevant polymorphism, the poly(A) microsatellite, which might influence either VDR mRNA stability³⁸ or VDR transcriptional activity.³⁹ This could be one of the reasons underlying the association between BsmI/TaqI polymorphisms with white stature. Alternatively, as BsmI and TaqI also localise to the so called block B, covering exons 3–9,³³ they might be tightly linked to other sites that are truly functional in this regard. More research is necessary to find out the causal variant(s).

The association between the VDR bT haplotype with height in white subjects was most significant in postmenopausal women, marginally significant in pre-menopausal and peri-menopausal women and not significant in men. This suggests the association of VDR bT haplotype with white height may be modified by sex and menstrual status, such that the “height retaining” effect of bT haplotype might be specific to females and more obvious after menopause. Our results regarding postmenopausal height seem to be inconsistent with the data of Uitterlinden et al,⁴⁰ which did not show evidence of association between the VDR baT haplotype with height in postmenopausal women (table 2 of that paper⁴⁰). However, our bT haplotype is the combination of baT and bAT haplotypes, thus the association outcome of bT haplotype with height could be different from that of baT. In addition, the height data in the study of Uitterlinden et al were not adjusted by the confounding factor of fracture in their sample, thus the potential “height retaining” effect of baT may be masked by its association with vertebral fracture, which leads to decreased height.

Another commonly studied SNP in the VDR gene, the FokI polymorphism at the first translation initiation site, was an independent marker in the LD breaking region.^32,33 Our result supported this notion by showing that there was no evidence of LD between FokI with any of the other three markers. The absence of association with height at FokI locus seemed to disagree with an earlier study,¹³ which found that FokI was associated with stature in 108 prepubertal and pubertal homozygous β thalassaemic patients.

The −4817A/G polymorphism was found in intron 1 of VDR, and has never previously been studied. For this marker, more heterozygotes than would be predicted by Hardy-Weinberg equilibrium were observed; such deviation from equilibrium may result from negative assortment mating with respect to the genotypes at this locus—that is, the mating between individuals with dissimilar genotypes (AA with GG in this case) occurs at higher than random frequencies. However, other possibilities may also exist, such as genotyping error or population stratification.⁴¹ So far, there is no consensus on the effects of Hardy-Weinberg disequilibrium on population based association results.⁴¹ This study used the approach of family based QTDT, and such disequilibrium will not affect our analyses. Firstly, the essence of the TDT is to compare the frequencies of transmitted and non-transmitted alleles of heterozygous parents within family.⁴² Whether the tested markers are in Hardy-Weinberg equilibrium or not does not affect the validity of the TDT. This is also true in linkage study, which attempts to capture co-inheritance patterns of trait and marker loci in families. Second, Genehunter, the software used to construct haplotypes in our nuclear families, does not assume Hardy-Weinberg equilibrium,²⁷ thus Hardy-Weinberg disequilibrium at the loci should not affect the haplotype reconstruction in our study.

Although there was no evidence of within family association with height for the −4817A/G polymorphism, we did find significant total association at this locus. This may be due to that total association confers more power and can be more sensitive than within family association given the absence of population stratification. As this locus is located between exon 1B and 1C, which are covered by block C,³³ it might be in LD with causal variants leading to alternative splicing of VDR transcripts⁴³ and thus indirectly influence height variation. However, whether this locus is truly associated with adult height awaits further investigation.

Linkage signals were found for all of these four markers, spanning ∼39 kb across the VDR gene. However, except for BsmI and TaqI, the association outcomes with height were different among them (table 5). This is not unexpected because linkage analysis (largely based on sibling pairs in this study) can rarely narrow candidate regions beyond the level of ∼10 cM.⁴⁴ There are generally not enough meioses in pedigrees of a few generations to detect recombination events between closely spaced loci under complex polygenic inheritance. Different associations with height were consistent with the observation that the degrees of LD vary across the SNP loci studied.

In summary, with a large homogenous white sample and the robust QTDT approach, this study demonstrates that VDR gene polymorphisms are significantly linked to and associated with height in white adults. Particularly, it highlights the importance of VDR polymorphisms near the 3′ UTR in final height variation in whites. A population based association test also supported this notion, and suggested that such a genetic effect might be specific to women and influenced by menstrual status. Nevertheless, to confirm our findings, replication studies and meta-analyses may be needed. Moreover, the exact “causal alleles” within VDR that influence final stature are still unknown. Therefore, it is necessary to perform further statistical and molecular genetic studies (including functional studies) in order to better understand the role of VDR in human stature.