A Whole-Genome Linkage Scan Suggests Several Genomic Regions Potentially Containing Quantitative Trait Loci for Osteoporosis

Osteoporosis Research Center (H.-W.D., F.-H.X., Q.-Y.H., H.S., H.D., T.C., Y.-J.L., Y.-Z.L., J.-L.L., H.-T.Z., K.M.D., R.R.R.) and Department of Biomedical Sciences (H.-W.D., F.-H.X., Q.-Y.H., H.S., H.D., Y.-J.L., Y.-Z.L., J.-L.L., H.-T.Z.), Creighton University, Omaha, Nebraska 68131; and Laboratory of Molecular and Statistical Genetics (H.-W.D.), College of Life Sciences, Hunan Normal University, ChangSha, Hunan 410081, People’s Republic of China

L OW BONE MINERAL density (BMD) is an important risk factor for osteoporotic fractures (OF), and osteoporosis is mainly characterized by low BMD (1)(2)(3). In the United States alone in 1995, osteoporosis results in more than 1.3 million OF a year, with an estimated direct cost of about 14 billion dollars (4). Extensive data have established that BMD variation is under strong genetic control with heritability estimates ranging from 0.5-0.9 (5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17). Recently, extensive molecular genetic studies (18 -29) have been launched to search for genes underlying BMD variation. The extensive results from population association studies of candidate genes so far have been largely inconsistent (30 -32). It is known that population admixture may induce false positive (33) or false negative (34) results in association studies of candidate genes.
Compared with the extensive association studies, only a limited few studies of alternative approaches such as linkage studies (20,23,24,35,36) and transmission disequilibrium test (37,38) in humans and quantitative trait locus (QTL) mapping (27-29, 39 -41) in mice have been conducted to identify genomic regions underlying population BMD variation. The 7 pedigrees studied (35) (identified through probands with low BMD) in an autosomal genome scan contain 74 sibling pairs genotyped for 330 markers. A total of 153 sibling pairs (randomly ascertained with respect to BMD) were studied in the genome-wide scan (36) using 367 markers (with an average heterozygosity of ϳ0.77). A total of 595 sibling pairs (randomly ascertained) were used (20), with 270 markers (average heterozygosity being ϳ0.70) genotyped throughout the genome. Sibling pair approach for linkage studies is usually of limited power when sibling pairs are ascertained randomly (42)(43)(44). For example, more than 8000 randomly ascertained sibling pairs are necessary in a whole genome scan to detect a major locus with a heritability as large as 30% (44). The limited power for the small samples of the previous sibling pair linkage studies and the few but the largely inconsistent empirical linkage findings warrant that independent studies with alternative linkage approaches and larger samples be conducted. Generally speaking, the linkage detecting power increases rapidly with the marker heterozygosity, reduced genotyping error rate, and marker density that may reflect expected genetic distance between markers and QTLs (45,46). Sampling through extreme probands may improve the linkage power (47,48).
In this study, for a relatively large sample of human pedigrees that was recruited through probands with low BMD values, we genotyped 380 highly heterozygous markers across the whole genome for each subject. We performed a wholegenome linkage analysis to identify genomic regions potentially important for BMD variation. exclusion criteria to be detailed below) were included in the analyses. All the study subjects were Caucasians of European origin as ascertained by questionnaires administered by research nurses. A total of 53 pedigrees with 630 subjects (248 males and 382 females) from 2-4 generations were analyzed. There are more female than male subjects, largely due to the higher propensity to volunteer in a genetic study of this type among females. Among the study pedigrees, 1 is 2-generation, 41 are 3-generation, and 11 are 4-generation pedigrees. The pedigrees vary in size from 3-99 individuals, with a mean of 11.7 (Ϯse ϭ 2.4). The pedigrees were identified through a proband having BMD Z-scores Ϫ1.28 or less at the hip or spine. Hence, the probands were selected from the bottom 10% of the population BMD variation with the intended purpose of achieving higher statistical power than random sampling (47,48). BMD values expressed as Z-scores adjust for age, gender, and ethnic difference in general referent healthy populations. The exclusion criteria for the study subjects were a history of (including past as well as current disease conditions, unless otherwise specified) the medical conditions that were detailed previously (23,49). The exclusion criteria were assessed by nurse-administered questionnaires and/or medical records and applied most rigorously to potential subjects contacted between ages 25 and 50. About 5.1% of the total people screened were excluded from our recruitment due to their meeting at least one of the exclusion criteria.

Genotyping
For each subject, blood (20 cc) was drawn into lavender cap (EDTA containing) tubes by certified phlebotomists. DNA was extracted by employing a kit (Puregene DNA Isolation Kit, Gentra Systems, Inc., Minneapolis, MN; catalog no. D-5000) following the procedures detailed therein. DNA was genotyped using fluorescently labeled markers commercially available through PE Applied Biosystems (ABI PRISM Linkage Mapping Sets, version 2, Norwalk, CT), as we did before (49,50). A genetic database management system (GenoDB) (51) was employed to manage the genotype data for linkage analyses. PedCheck (52) (available at http://watson.hgen.pitt.edu/register/soft_doc.html) was employed for checking the conformation to Mendelian inheritance pattern at all the marker loci and for checking the relationships of family members within pedigrees. The genotyping error rate, determined by the procedures described earlier (59,51,53), was about 0.3%. A total of 380 markers (including 362 on autosomes) were successfully genotyped. These markers have an average population heterozygosity of approximately 0.79.

Measurement
BMDs of spine, hip, and wrist were measured by a Hologic 1000, 2000ϩ, or 4500 dual energy x-ray absorptiometry scanner (Hologic Corp., Waltham, MA). All machines are calibrated daily, and long-term precision is monitored with external spine and hip phantoms. Shortterm precision in humans as determined by the method of Glü er et al. (54) with repeated measurements on at least 27 subjects for a single type of machine is 0.8% for spine BMD, 0.6% for hip BMD, 1.1% for wrist BMD on Hologic 1000; and 0.8% for spine BMD, 1.0% for hip BMD, 2.6% for wrist BMD on Hologic 2000ϩ; and 0.9% for spine BMD, 1.4% for hip BMD, and 2.3% for wrist BMD on Hologic 4500. We chose to study BMD because BMD is the measure most closely correlated with fracture risk (55). In addition, more than 86% and about 98% of BMD variation are attributable to bone mineral content (BMC) variation at the spine and hip, respectively (56), and the majority of the variation of BMD and BMC at the spine, hip, and wrist can be accounted for by one principle component in factor analyses (11). Our analyses (57,58) show that the genomic regions for BMD and bone size are largely different, whereas those for BMD and BMC are largely the same. For the spine, our quantitative phenotype was combined BMD of L 1-4 . For the hip, it was total BMD of the femoral neck, trochanter, and intertrochanteric region provided by the instruments. For the wrist, it was ultradistal BMD. All dual energy x-ray absorptiometry machines report BMD in g/cm 2 . Weight was measured at the same visit when the BMD measurements were taken. Data obtained from different machines are transformed to a compatible measurement by an algorithm as in our previous studies (59), and members of the same pedigree are usually measured on the same type of machine.

Statistical analyses
A variance component linkage analysis (60 -62) for quantitative traits was performed. The analysis is based on specifying the expected genetic covariances between arbitrary relatives as a function of the identity by descent at a given marker locus. The analysis considers the phenotypic and genetic information from all pedigree members simultaneously. The analysis assumed joint multivariate normality of phenotypic values, additive genetic effects, and no interaction between genes and the residual variation. The common familial environmental effects were assumed to be negligible, which is reasonable and supported by previous studies (7-9, 11, 12). The program employed was SOLAR (Sequential Oligogenic Linkage Analysis Routines) (62), which is available at http://www.sfbr.org/sfbr/public/software/solar/. The ascertainment scheme of pedigrees based on the low BMD values of probands was accounted for in the analyses by identifying to the program the probands and their BMD values for each pedigree. The built-in modules of the SOLAR program can account for the ascertainment scheme by using cut-off BMD values of the probands and conditional likelihoods in LOD score computation.
In linkage analysis, age, sex, and weight were adjusted as covariates for raw BMD values (not the Z-scores), as these generally affect BMD variation (10,63) and tested to be significant in our sample in a statistical screen for important covariates. Z-scores of BMD reported by the Hologic machines adjust for age and sex (but not weight) effects in referent populations provided by the manufacture of the scanners and our BMD phenotype data for linkage analyses are those adjusting the raw BMD data for age, sex, and weight according to their effects in our own sample. Analyses were also performed without adjusting for some or any of these covariates. Adjustment for significant covariates in genetic analyses can generally increase the genetic signal to noise ratio in linkage detection by decreasing the proportion of the residual phenotypic variation attributable to random environmental factors (22). Comparison of the analyses with and without adjustment for a significant covariate may shed light on the genomic regions identified as to their direct importance for the trait per se or indirect importance via the influence on the covariates only. The BMD data were tested by graphical methods (64) and found not to deviate from normal distributions. The variance component analyses implemented in SOLAR are quite robust to reasonable violations of normality of the data (65). Using SOLAR, two-point, and multipoint linkage analyses were performed. When a putative QTL is suggested, the proportion of phenotypic variation attributable to this QTL can be estimated by SOLAR. The estimate is usually an upper bound of the genetic effect due to the locus (66).

Results
The basic characteristics of the study subjects stratified by age and sex are summarized and the data are available upon request. Detailed information about the family structure of the study pedigrees has been summarized earlier (49). Mainly, there are 1249 sibling pairs, 1098 grandparent-grandchild pairs, and 2589 first cousin pairs, etc. The average Zscores (Ϯse) of spine, hip, and wrist BMD are Ϫ0.12 Ϯ 0.05, Ϫ0.22 Ϯ 0.04, and 0.09 Ϯ 0.05, respectively, for all the subjects in the sample, reflecting the sampling scheme of our pedigrees which were ascertained through probands of extremely low BMD values at spine and/or hip.
The results for the chromosomes that had maximum LOD scores (MLS) close to or greater than 2.0 in multipoint linkage analyses at any one of the three skeletal sites were summarized in Fig. 1; the results for the other chromosomes were summarized in Fig. 2. Because the current version of SOLAR does not handle multipoint analyses for X chromosome, we plotted two-point LOD score results for the X chromosome and presented them in Fig. 2 to convey the linkage signal pattern. For the results presented, BMD were adjusted with multiple regression for significant covariates of age, sex, weight, and height. The results (not presented but available upon request) for the linkage analyses for BMD values that were not adjusted for the covariates showed similar pattern in LOD scores and qualitatively the same results, suggesting the importance (if indeed exists) of the genomic regions identified here is direct on BMD rather than indirect through a covariate analyzed. In fact, our whole genome linkage scan results (49) for body mass index (which is weight/height 2 ) revealed genomic regions for body mass index different from those revealed here for BMD.
For the spine BMD, four genomic regions were identified that may contain putative QTLs. Specifically, an MLS of 3.08 at 152 centimorgans (cM) from p terminal (pter) on chromosome 4 in multipoint linkage analyses and a LOD score of 2.12 at the marker D4S413 in two-point linkage analyses were achieved in the genomic region 4q31-32. An MLS of 2.96 at 169 cM from pter on chromosome 12 in multipoint analyses and an MLS of 2.17 at D12S1723 in two-point analyses were achieved in 12q24. An MLS of 2.43 at 103 cM from pter on chromosome 13 in multipoint analyses and a LOD score of 1.77 at D13S285 in two-point analyses were achieved in 13q33-34. An MLS of 1.93 at 0.0 cM from pter on chromosome 7 in multipoint analyses and a LOD score of 2.28 at D7S531 in two-point analyses were achieved in 7p22. In addition, in two-point analyses, a LOD score of 1.6 was achieved at D15S165 in 15p11; however, the MLS is only 0.76 in the genomic region (14cM from pter on chromosome 15). At the spine, for the putative QTL with the strongest evidence of linkage on 4q31, approximately 30.3% of BMD variation (after adjusting for age, sex, weight, and height) may be attributable to this locus.
For the hip BMD, three genomic regions were identified that may contain putative QTLs. An MLS of 2.29 at 170 cM from pter on chromosome 10 in multipoint analyses and a LOD score of 1.97 at D10S1651 in two-point analyses were achieved in 10q26. In addition, in two-point linkage analyses, a LOD score of 1.69 was achieved at D12S368 on 12q13 and a LOD score of 1.58 was achieved at D17S1857 on 17p11, with corresponding multipoint MLSs being 0.99 and 1.2, respectively. At the hip, for the putative QTL with the strongest evidence of linkage on 10q26, approximately 29.2% of BMD variation (after adjusting for age, sex, weight, and height) may be attributable to this locus.
For the wrist BMD, four genomic regions were identified that may contain putative QTLs. An MLS of 2.26 in multipoint analyses at 158 cM from pter on chromosome 4 and a LOD score of 2.53 at D4S413 in two-point analyses were achieved in 4q32. An MLS of 1.87 in multipoint analyses at 16 cM from pter on chromosome 9 and an LOD score of 1.74 in two-point analyses at D9S285 were achieved in 9p22-24. An LOD score of 1.99 was achieved at D17S1852 in 17p13, and a LOD score of 1.82 in two-point analyses was achieved at D3S1297 in 3p26, with corresponding MLSs of 1.48 and 1.25 achieved in multipoint analyses in the respective genomic regions. At the wrist, for the putative QTL with the strongest evidence of linkage on 4q32, approximately 25.5% of BMD variation (after adjusting for age, sex, weight, and height) may be attributable to this locus.

Discussion
Some of the regions we identified are supported by earlier findings or contain important candidate genes (Tables 1 and   2). For the spine BMD, four genomic regions (4q31-32, 7p22, and 13q33-34, and 12q24) were identified. Drake et al. (2001) (35) detected a QTL region for femur BMD with an LOD score greater than 2.3 in mice in regions homologous to human chromosome12q24. IGF1, TBX3, and TBX5 genes that are of importance in BMD metabolism are located in the 12q24 region (Table 1). Several studies have suggested the importance of IGF1 in osteoporosis (67,68). The genomic region 13q33-34 has previously (36) been linked to forearm BMD with a LOD score of 1.67 (Table 2). Potential candidate genes in this region include COL4A1 and COL4A2 (Table 1). Additionally, a candidate gene IL6 was mapped to 7p22. Duncan et al. (24) obtained a two-point LOD score of 1.7 at D12S83 for lumbar spine BMD. Our results for D12S83 do not provide supporting evidence of linkage. However, for the marker D12S368 (in 12q13) that is 8.2 cM from D12S83, a two-point LOD score of 1.69 is obtained for the hip BMD in our analyses. The genomic region near D12S83 and D12S368 contain a number of candidate genes, including the prominent vitamin D receptor gene (Table 1). For the three genomic regions (10q26, 12q13, and 17p11) identified for hip BMD, in mice in a genomic region homologous to human chromosome 10q26, Klein et al. (26) identified a QTL region for whole body BMD (P ϭ 0.0093) and Beamer et al. (39) identified a QTL region for vertebral BMD with a LOD score of 5.01. For the four genomic regions (3p26, 4q32, 9p22-24, and 17p13) identified for wrist BMD, 4q32 coincides with the region identified for spine BMD, and 17p13 coincides with the region identified for hip BMD in our sample. These results suggest that some genetic factors underlying BMD variation at different skeletal sites are located closely in the human genome or are shared. In the region homologous to human chromosome 3p26, Beamer et al. (39) identified in mice a QTL region for femoral BMD with a LOD score of 4.56. Devoto et al. (35) reported a LOD score of 2.34 near D17S261 for hip BMD, which is in the region 17p11-13 identified for hip and wrist BMD in our study. In the region homologous to human chromosome 17p11-12 in mice, Beamer et al. (39) identified a QTL region for femoral and vertebral BMD variation with LOD scores of 6.76 and 2.98, respectively. A QTL for femur peak bone mass in the mouse genomic region homologous to human chromosome 17p11-13 was identified with a LOD score of 10.8 (40).
However, comparing the few earlier whole-genome linkage studies of BMD in humans (20,35,36) and our study here, most of the results cannot be confirmed by each other. Table  2 outlines the comparison of the earlier results and ours. Our results seem to replicate the linkage findings on 4q and 13q ( Table 2) because confirmation of a previously reported link-age in an independent study involves P value criteria different from the whole-genome linkage scan (69). Replication of whole genome linkage studies is difficult due to the potential polygenic nature of inheritance and the limited power associated with the current whole-genome linkage studies (23). A preliminary comparison of mice and human mapping studies for osteoporosis has been made earlier (70). Mitchell et al. (71) studied serum osteocalcin concentrations, which is associated with bone loss and turn over rates that are related to BMD. A whole genome search in 429 subjects from 10 randomly ascertained large multigeneration pedigrees identified linkage of genomic regions on 16q and 20q to serum osteocalcin concentrations (71). These two regions differ from those regions identified for BMD variation in earlier whole genome studies (20,35,36) and our results here.
The power and the robustness of the linkage results critically depend on the study sample sizes employed, among several other things. Generally, the larger the sample, the higher the power and thus the more robust the results. For the same number of subjects involved, pedigree studies may have dramatically increased power over sibling pair studies (72). The LOD scores reported by earlier whole-genome linkage scan studies with smaller samples are sometimes higher than our current study (Table 2). However, this does not  Unless otherwise specified, the candidate genes (if indicated) in the genomic regions, together with their cellular functions, are identified via the skeletal gene database (80).
indicate that those studies are necessarily more powerful and their results are more robust than our current study. Multiple testing and polygenic effects, among other things, may easily yield statistical fluctuation of LOD scores, particularly in small samples (23). A concrete empirical example that demonstrates statistical fluctuation of LOD scores in small study samples is as follows. A sibling pair study with 374 sibships reported a maximum LOD score of 3.50 near D11S987 for linkage to femoral neck BMD variation (73). However, a subsequent linkage study (20) with 595 sibling pairs reported a LOD score less than 2.2 at D11S987. Replication studies with even larger samples for the genomic region near D11S987 yielded no evidence of linkage to this region (23,74,75). Hence, the higher LOD score of the study (73) does not imply that the study (73) is more powerful and the results are more robust than other studies (20,23,74,75) with larger samples but smaller LOD scores. None of our LOD scores meet the criteria (69) for definite linkage. However, the larger sample sizes of our study may render our results statistically more robust than some of the earlier studies with smaller samples. Simulation studies specific for our study pedigrees show that if using a LOD score of 3.0 as a significant criterion, we may achieve statistical powers of 42.5%, 76.5%, and 94.5% for a QTL of a heritability of 15%, 20%, and 25%, respectively.
Age, sex, and weight are well-known major factors that are associated with BMD changes (63). Age (as age and age 2 terms), sex, and weight are adjusted as covariates, which, by the basic principle of statistics and statistical modeling (64), effectively controls for their effects on BMD in our study sample. Age effects in sibling pair studies may be more rigorously controlled experimentally by selecting siblings similar in age and of ages between 20 and 45 (20). However, even for the age range of 20 -45, BMD still changes with age (60), and the age of peak bone mass reached may vary from individual to individual. Hence, statistical testing and adjustment of age effects is still necessary in sibling pair studies (36). Using multigeneration large pedigrees for the linkage study of BMD, the age range covers those of peak bone mass reached and those in the elderly resulting from the subsequent bone loss after the peak bone mass is reached. The natural process of menopause in women and other potential age-related physiology effects such as osteophytes in the lumbar spine may compromise our ability to fully adjust age effects on BMD. Future measures such as obtaining information about other potentially significant factors For those markers that were reported by the other studies but not genotyped in ours, we reported the LOD scores for our markers (the identity given) that are closest to them and give the genetic distance (found through the database "Comprehensive human genetic maps of Marshfield" at http://research.marshfieldclinic.org/genetics/Map_Markers/maps/indexmap.html) between the markers in other studies and the closest markers of ours in parentheses. Markers D18S42 and CD3D were not mapped onto the Genethon and Marshfield human genetic linkage maps (http://www.ncbi.nlm.nih.gov/PMGifs/Genomes/humansearch.html), and hence we are not able to identify our closest markers to them. MLS denotes for maximum LOD scores in the genomic regions.
such as menopause, measuring lumbar spine BMD with lateral scans, detailed recording of therapeutic effects (such as with estrogen in postmenopausal women) may ameliorate potentially confounding effects associated with large age ranges of study subjects and other factors. BMD in the elderly is most relevant to osteoporotic fractures that occur mainly after age 50 (3). Hence, a study that covers a large range of ages (particularly those for the elderly) has a scope that may be able to reveal potentials QTLs not only relevant to peak bone mass but also relevant to bone loss. Our study and those that confine to a more narrow age range aimed at QTL mapping for peak bone mass (20) may be able to complement each other in terms of disentangling genes for peak bone mass and those (if different) for subsequent bone loss. QTL mapping that involves measurement of BMD at different development stages (29) for the same individual organism may provide a way to distinguish peak bone mass QTLs from bone loss QTLs.
The genomic regions identified in this and all previous studies need to be subject to extension studies with more samples and denser markers for confirmation. The consistency and inconsistency of the genomic regions identified may serve as a basis for further endeavors and exploration of some focused genomic regions with larger samples and/or alternative methods (38). Implementation and comparison of studies with various designs and populations and various species may help facilitate confirmation of the results obtained. While QTL mapping in model organisms such as mice is undoubtedly useful, the results found in one species may not always transcend across to other species as shown by several known examples in studies of other phenotypes such as obesity (70). This may not be all unexpected at least because the natural selection challenges and thus evolutionary trajectories of genetic architecture of various species are not all the same during the long time independent evolutionary history. Once linkage to a genomic region is confirmed, subsequent fine mapping, a daunting task, will be pursued. Theoretically, examining the pattern of linkage disequilibrium using a dense panel of markers within positively identified genomic regions can narrow the relatively large genomic regions (ϳ30 cM) identified in linkage analyses to regions less than 1 cM (76,77), a size of the genomic region that is amenable for physical cloning of QTLs. If our ultimate goal is to find genes important for osteoporotic fractures and BMD is just used as a surrogate measure, the genes to be identified important for BMD may eventually need to be tested for their significance for and relevance to osteoporotic fractures (3,78).