Linkage of osteoporosis to chromosome 20p12 and association to BMP2.

Osteoporotic fractures are a major cause of morbidity and mortality in ageing populations. Osteoporosis, defined as low bone mineral density (BMD) and associated fractures, have significant genetic components that are largely unknown. Linkage analysis in a large number of extended osteoporosis families in Iceland, using a phenotype that combines osteoporotic fractures and BMD measurements, showed linkage to Chromosome 20p12.3 (multipoint allele-sharing LOD, 5.10; p value, 6.3 x 10(-7)), results that are statistically significant after adjusting for the number of phenotypes tested and the genome-wide search. A follow-up association analysis using closely spaced polymorphic markers was performed. Three variants in the bone morphogenetic protein 2 (BMP2) gene, a missense polymorphism and two anonymous single nucleotide polymorphism haplotypes, were determined to be associated with osteoporosis in the Icelandic patients. The association is seen with many definitions of an osteoporotic phenotype, including osteoporotic fractures as well as low BMD, both before and after menopause. A replication study with a Danish cohort of postmenopausal women was conducted to confirm the contribution of the three identified variants. In conclusion, we find that a region on the short arm of Chromosome 20 contains a gene or genes that appear to be a major risk factor for osteoporosis and osteoporotic fractures, and our evidence supports the view that BMP2 is at least one of these genes.

Introduction

Osteoporosis is a common disease characterized by low bone mineral density (BMD) and manifesting clinically with fragility fractures of the hip, spine, and other skeletal sites. It is generalized in its most common form, affecting the elderly, both sexes, and all racial groups, although postmenopausal women are at highest risk (Peacock et al. 2002). There are over 1 million osteoporotic fractures per year in the United States alone, primarily in women, and the direct medical costs exceed US$10 billion annually (Ray et al. 1997). One of the challenges in medicine today is to identify those who are at high risk for osteoporosis before they suffer fractures or lose significant bone mass.

Various factors may distinguish those who develop osteoporotic fractures from those who do not, but the most important appears to be BMD values, both attained in young adults as well as with increasing age. Firstly, peak bone mass (the highest BMD, attained in young adulthood) may represent an important measure of predisposition to osteoporosis. Secondly, the rate of postmenopausal bone loss may be higher in some than in others. There is abundant evidence for a genetic contribution to BMD variation, especially the peak bone mass, but a large genetic contribution to BMD at older ages has also been demonstrated. Environmental factors such as diet, medications, and physical activity may also determine the ultimate BMD. Furthermore, the rate of bone loss, bone size and structure, and propensity to fall are all factors with genetic components and all contribute to the risk of osteoporotic fractures beyond BMD itself (Peacock et al. 2002).

Numerous candidate genes, selected on the basis of current knowledge of bone biology, have been tested for association to BMD and to osteoporotic fractures. A polymorphism in the Sp1 transcription factor-binding site in the first intron of the collagen 1A1 (COL1A1) gene has shown the most consistent association to osteoporosis, although varying between study groups and populations (Efstathiadou et al. 2001; Mann et al. 2001). A few genome-wide linkage scans have been reported by those searching for new genes contributing to BMD variation and skeletal geometry (Devoto et al. 1998; Niu et al. 1999; Koller et al. 2000, 2001; Deng et al. 2002; Karasik et al. 2002; Wilson et al. 2003), and several linkage studies analyzing only specific chromosomal locations have also been conducted (Spotila et al. 1996; Koller et al. 1998; Duncan et al. 1999; Carn et al. 2002). However, none of the loci reported have met criteria for genome-wide significance for linkage, and the osteoporosis genes corresponding to these regions have not yet been isolated. The roles in the common forms of osteoporosis of the genes encoding low-density lipoprotein receptor-related protein 5 (LRP5) and the osteoclast-specific vacuolar proton pump (TCIRG1) on Chromosome 11q12-13, which cause three Mendelian BMD-related diseases/traits (Frattini et al. 2000; Kornak et al. 2000; Gong et al. 2001; Boyden et al. 2002; Little et al. 2002; van Hul et al. 2002), also remain to be shown.

Here we present the results of the genealogic approach applied to osteoporosis in Iceland (Gulcher et al. 2001a, 2001b). The aim of the study described here was to find genes contributing to the risk of osteoporotic fractures and its immediate precursor, low BMD. However, instead of simply mapping a quantitative trait locus for BMD itself, we chose to focus on the clinically most relevant manifestation of osteoporosis, low BMD and osteoporotic fractures. Low BMD is the single best predictor of osteoporotic fractures (Marshall et al. 1996) and is most commonly used both as a diagnostic criterion and as the basis for treatment. BMD is, however, only a surrogate marker for osteoporotic fractures and is influenced by many factors. Furthermore, osteoporotic fractures also occur in individuals who do not fall within low BMD range. Therefore, we defined a novel mixed phenotype for osteoporosis for use in this study, combining low BMD values and osteoporotic fractures themselves. By combining age-matched BMD values at both the hip and the spine in the same analysis, we attempted to find the underlying genetic basis for generalized osteoporosis, instead of a skeletal site-specific osteoporosis. Furthermore, although BMD is normally distributed at all ages in the population and is therefore a quantitative trait with presumed polygenic inheritance, it does not necessarily imply that the genes controlling low BMD in families, significant enough to lead to osteoporotic fractures, are the same as those that control the normal distribution of BMD in the population. Therefore, an affected-only approach to linkage analysis, using BMD values below a certain cut-off in combination with osteoporotic fractures, might be an effective approach to the mapping of genes that predispose to osteoporosis. We found highly significant linkage to Chromosome 20p using this approach, discovered variants in the bone morphogenetic protein 2 gene (BMP2) that are associated with osteoporosis, and confirmed the BMP2 association in a cohort of Danish osteoporosis patients.

Results

Linkage Analysis

Linkage analysis, using multipoint allele-sharing methods, was conducted for four phenotype criteria in 207 extended Icelandic osteoporotic families, containing 1,323 study individuals (see Materials and Methods). Descriptions of the phenotypes and pedigree sets along with a summary of the linkage results are provided in Table 1. We first investigated osteoporosis as it is broadly defined (moderate pedigree set). This is the most generalized osteoporosis phenotype, including individuals with a hip and spine BMD approximately one standard deviation (SD) or more below the average and/or those with osteoporotic fractures and/or those receiving bisphosphonate treatment for osteoporosis. The most prominent peak was on Chromosome 20 at D20S905 (19.90 cM) with an allele-sharing logarithm of the odds (LOD) score of 3.39 (p value, 4.2 × 10−5), with four other locations, achieving a LOD score of 1.5 or greater: 16q, two on 18p, and 21q (Figure 1; Table 1).

Table 1

Genome-Wide Scan Results and Fine-Mapping on Chromosome 20p12

The LOD scores above 1.5 and their genomic locations for any of four phenotypes analyzed using the genome-wide scan (GWS) framework marker set (1,100 microsatellite markers) and fine-mapping markers on Chromosome 20p12 are shown. Included are the number of informative families (#fam) for each phenotype and number of genotyped affected members (#aff) and relatives (#rel) in each analysis. The moderate, the hip, and the spine phenotypes include persons in the lower 16th percentile of combined hip and spine BMD, hip BMD, or spine BMD, respectively. The severe phenotype includes persons in the lower 10th percentile of combined hip and spine BMD. All BMD measurements are corrected for age, sex, and weight. The pedigree sets also include as affected those being treated for osteoporosis with bisphosphonates or having had an osteoporotic fracture at the relevant locations: hip fracture in the hip pedigree set, vertebral fracture in the spine pedigree set, or any osteoporotic fracture in the moderate or severe pedigree sets (see Materials and Methods). The number of those additional individuals who did not meet the BMD criteria alone for affection status were 3, 3, 30, and 20 for the hip, spine, moderate, and severe pedigree sets, respectively

Figure 1

Framework Linkage Scan

Framework linkage scan using 1,100 microsatellite markers for the severe (black line), moderate (red line), hip (green line), and spine (blue line) pedigree sets. The LOD score is on the y axis and the distance from the pter in Kosambi cM is on the x axis. Note that the LOD score scales are the same for all chromosomes except Chromosome 20.

We next used a definition of a more severe phenotype, with only individuals in the lower 10th percentile of BMD as affected members, but including fracture patients and patients treated for osteoporosis as in the previous run. The most striking feature of this scan was the increase in the peak on Chromosome 20p (Figure 1) that was again at D20S905, but now with a LOD score of 4.93. Two new LOD peaks were observed on Chromosomes 6p and 17p (Figure 1; Table 1). Compared to the LOD score peaks we had observed when analyzing the moderate pedigree set, for the severe pedigree set both of the peaks on 18p were greatly attenuated, the peak on 16p persisted, and the peak on 21p dropped by about a LOD of 0.5. Two additional analyses were conducted, one considering only hip osteoporosis (hip pedigree set) and the other considering only spine osteoporosis (spine pedigree set). In both cases, the strongest linkage was to Chromosome 20p, but with slightly higher LOD score (3.18) for the hip analysis (Figure 1).

The results at Chromosome 20p were encouraging, so we decided to genotype 30 additional markers in the region in order to increase the information on identity by descent sharing. At this increased fine-mapping density, the information on sharing was over 95% and the LOD score rose to 5.10 (p value, 6.3 × 10−7) at D20S194 (20.35 cM) in the severe phenotype. The LOD score increased for all phenotypes, to 3.99, 3.99, and 3.43 for the moderate, the hip, and the spine pedigree sets, respectively. After applying a Bonferonni adjustment for the four phenotypes tested, the p value for linkage is 1.6 × 10−6, which is genome-wide significant according to the dense marker threshold for multipoint genome scans (Lander and Kruglyak 1995).

Linkage Disequilibrium Mapping, Association Analysis, and BMP2 Variants

We investigated a 6.6 cM region surrounding our linkage peak that corresponds to a 1.7 Mb segment. This region, which includes up to a drop of 1 in LOD score from the linkage peak on the centromeric side and slightly less on the telomeric side, contains six known genes: BMP2, CHGB, LOC51605, C20orf154, C20orf155, and C20orf42 (Figure 2B). The BMP2 gene is a strong candidate gene, based on its role in bone formation and osteoblast differentiation (Wozney et al. 1988; Katagiri et al. 1994; Fujii et al. 1999). However, expression analysis showed that four of the genes in the region are expressed in bone marrow or in an osteoblast cell line (BMP2, C20orf42, C20orf154, and CHGB). Therefore, the expression data did not help us narrow the list of candidate genes.

Figure 2

The Chromosome 20p12 Linkage Region

(A) Chromosome 20 linkage scan for the narrower definition of osteoporosis (severe). The LOD score is on the y axis and the distance from pter in Kosambi cM is on the x axis.

(B) Region under the linkage peak is shown in more detail, including location of microsatellite markers (STRs), location of SNPs, and location of genes in the region and their direction of transcription. The legend bar indicates the distance corresponding to 100 kb.

To decide which gene in the focused region is most likely to contribute to osteoporosis, we performed a case–control association study by typing many additional polymorphic markers. We used the severe-linkage phenotype as affection status and included all available genotyped patients, including additional familial cases and sporadic cases, in the analysis, 705 in total. BMP2 was our leading candidate, and we genotyped two microsatellites and 15 single nucleotide polymorphisms (SNPs) in the 15 kb region containing BMP2. The other five genes clustered in a 220 kb region. For this region of genes, we genotyped nine microsatellites and 26 SNPs. We genotyped additional markers across the entire region in case there were other unidentified genes or nongene regions associated with the disease. In total in the 1.7 Mb region, we genotyped 41 microsatellites, 99 SNPs that are polymorphic in the Icelandic population (some of the public SNPs are exonic), and 20 other SNPs that we identified by screening all exons and flanking intronic sequences of all six genes or expressed sequence tag (EST) matches. Figure 2B displays the location of the microsatellites and SNPs used in this analysis relative to the gene locations. Single-marker association and haplotype analyses showed the strongest association to the region of the BMP2 gene (Figure 3). The BMP2 gene appears to be the only gene within this region, despite our extensive efforts to find additional transcripts (see below). Linkage disequilibrium (LD) analysis across the region, including some SNPs added subsequently, identified about ten blocks of LD, ranging from 10 kb to 250 kb in size. Further blocks of LD are likely to exist in some gaps where no polymorphic markers were found, but these locations tend to be “gene-free.”

Figure 3

Association Analysis in the Focused Region

Association results for one to four consecutive microsatellite markers over the 1.7 Mb region. Only haplotypes with RRs over 1 are plotted. p Values are given on the left and megabase (Mb) locations at bottom.

In order to further understand the haplotype associations and to search for functional variations or SNPs that might capture the risk of osteoporosis, we sequenced the BMP2 gene in 188 patients and 94 controls. The sequenced region covered the 15 kb containing the BMP2 gene, including all exons, introns, and the promoter region, along with another 44 kb of DNA flanking the gene. We genotyped all patients and controls for 63 SNPs from among those identified within the 59 kb region, including SNPs that we had used in this region in the previous analysis. Three of the SNPs change amino acids in BMP2: Ser37Ala, due to a T to G transversion at nucleotide position 116 in exon 2; Ala94Ser, due to a G to T transversion at nucleotide position 287 in exon 2; and Arg189Ser, due to an A to T transversion at nucleotide position 224 in exon 3. The Ala94Ser variant was only found in one family. The Arg189Ser variant is very common and does not show any association to osteoporosis on its own. The Ser37Ala variant, however, changes a conserved amino acid and shows significant association to osteoporosis. It is a relatively rare variant, with allele frequency ranging from 3.0% to 4.9% in patients, depending on the phenotype, versus 0.8% in controls, yielding relative risks (RRs) in the range of 3.8 to 6.3 (Table 2); four of the seven phenotypes tested produce p values smaller than 0.001 after adjusting for the relatedness of the affected (see Materials and Methods). Given that Ser37Ala is only one of three missense variants identified in BMP2, we consider this result statistically significant even with adjustment for multiple comparisons.

Table 2

Association Results for BMP2 SNP Haplotypes and Ser37Ala Missense Variant

Shown are the BMP2 Ser37Ala missense variant, hapB and hapC, number of affected (#aff) and controls (#ctrl) in the analysis, RR (r), frequency (%) of haplotypes (allelic) in affected and controls, and p values. The osteoporosis (OP) phenotypes in the analyses are given to the left (see Materials and Methods for details). The BMD values are corrected for sex, age, and weight and represent the lowest 10th percentile in all cases. The osteoporosis severe phenotype is the severe definition used in the linkage run. Premenopausal and postmenopausal phenotypes consist of low BMD values at the hip or the spine. Only a subset of patients and controls were typed for the Ser37Ala variant in the Icelandic cohort. For Ser37Ala, the p values based on the Fisher's exact test, but not adjusting for relatedness among the affected, are presented (p value) together with p values that have made that adjustment (P adja). The p values for the haplotypes are based on likelihood ratio tests. p Values for the Icelandic data are two-sided, but without adjustment for multiple comparisons (p value) or adjusted for the multiple haplotypes tested using a randomization procedure as described in the text (P adjb). The p values for the Danish replication study are one-sided. NS, not significant

Since Ser37Ala is rare and can only account for a small fraction of the affected, we believed there had to be other at-risk variants of BMP2. Through our effort to search for these other variants, we found three SNP haplotypes covering the BMP2 gene that show the strongest association to osteoporosis (Table 2; Figure 4A). The three haplotypes, hapA, hapB, and hapC, defined by a total of six SNPs, are essentially independent of each other (pairwise R, less than 0.01). The Ser37Ala variant, however, seems to have arisen on the background of hapA so that every chromosome that carries Ser37Ala also carries hapA, but not vice versa. Moreover, it appears that the excess of hapA in the affected can be explained by Ser37Ala; i.e., chromosomes with hapA but not Ser37Ala do not have significant excess in the affected versus controls. The hapB is defined by three SNPs—TSC0271643 (T), P9313 (T), and rs235764 (G)—and spans from within intron 2 of BMP2 and 263 kb upstream of the gene, going across LD blocks (Figure 4). The upstream SNP is in the region of the most significant single-point microsatellite association. HapC is defined by two SNPs, rs1116867 (A) and D35548 (T), downstream of BMP2 and spans 18 kb (Figure 4). HapC is entirely contained within one LD block. There is, however, significant LD between some markers across the hapC LD block and the BMP2 LD block. BMP2 is the only known or predicted gene within the segments spanned by the two haplotypes, as previously noted.

Figure 4

BMP2 Haplotype Region and LD

(A) The location of haplotypes, markers, and the Ser37Ala missense variant is shown in relation to the location of BMP2 exons and to the LD blocks in (B).

(B) Two measures of LD are shown: D′ values on the upper-left side of the plot and p values on the lower-right side. Scales for the LD strength are provided for both measures to the right. This graph shows uniformly distributed SNPs in a 370 kb region in and around the BMP2 gene. The BMP2 block and the next block downstream represent 53 kb.

To explore the role of Ser37Ala and hapB and hapC in relation to various forms of osteoporosis, we used several other osteoporotic phenotypes in the association analysis, including osteoporotic fractures, whole-body BMD, and pre- and postmenopausal BMD, each analyzed separately (see Table 2). The Ser37Ala variant and hapC are associated with all low BMD phenotypes and with osteoporotic fractures. In fact, the RR for low BMD is higher in premenopausal women than postmenopausal women, indicating that BMP2 may influence the attainment of peak bone mass. HapB, on the other hand, seems more associated with osteoporotic fractures than with BMD phenotypes.

Removing patients who carry the Ser37Ala variant or the disease-associated haplotypes from the linkage analysis, leaving 152 affected in families informative for linkage, the LOD score drops from 5.10 to 3.55. Hence, these variants alone cannot account for the entire linkage signal. Even taking into consideration that there is usually a substantial upward bias in estimation of locus-specific effects from genome scans (Goring et al. 2001), we believe there are likely to be other variants, probably rare but possibly having very high penetrance, which have yet to be identified. We also cannot rule out that other genes in the region are conferring osteoporosis risk independently or interacting with BMP2. We have not yet determined the underlying functional variations of disease-associated hapB and hapC. Most importantly, as anonymous haplotypes, proper interpretation of the statistical significance of hapB and hapC has to take multiple comparisons into account. By performing a randomization test (randomizing the patients and controls) for a procedure that considers all haplotypes within the LD block that contains hapC, the smallest adjusted p value is 0.003 for whole-body BMD (see Table 2), 125 times larger than the unadjusted p value of 2.4 × 10−5. Since this adjusted p value still has not taken into account the relatedness of the patients, the multiple phenotypes tested, and the other LD blocks in the region, we consider hapC at best marginally significant with the Icelandic data alone. HapB goes across LD blocks, and the randomization procedure is performed by considering all possible haplotypes within the region spanned by those blocks, leading to an even larger adjustment factor. As a result, none of the adjusted p values for hapB are significant (see Table 2). Hence, instead of haplotypes already established to confer increased risk, we considered hapB and hapC only as candidate at-risk haplotypes, identified based on the Icelandic material that required confirmation from a replication study.

Haplotype Confirmation in Danish Cohort

For replication, we used a cohort of Danish postmenopausal women with persistently low BMD and a group of Danish postmenopausal osteoporotic fracture patients. We only typed those SNPs that comprised hapB and hapC as well as the Ser37Ala variant. The results for the Danish study are also shown in Table 2. For all six tests, two phenotypes times three variants, the patients have a higher frequency of the variant compared to the controls, and the estimated RRs are overall comparable to those observed in the Icelanders. However, the p values are higher in general due to the smaller sample sizes of the Danish cohort. Nonetheless, hapC gives a p value of 0.0038 for low BMD, which is significant even after adjusting for the six tests performed. The Ser37Ala missense variant is nominally significant for osteoporotic fractures, with an estimated RR of 4.2, and hapB is nominally significant for low BMD, with an estimated RR of 2.1.

Another way of handling the problem of testing multiple variants is to consider them as a group. Define a composite haplotype, hapX, as a chromosome that carries at least one of either the Ser37Ala missense variant, hapB, or hapC. For the Danish cohort, hapX is estimated to have a frequency of 16.3% in the low BMD patients versus 7.3% in controls, giving a one-sided p value of 0.004 with an RR of 2.5. The frequency of hapX in the osteoporotic fracture patients is estimated to be 15.3%, giving a one-sided p value of 0.037 with an RR of 2.2.

Screening for Genes Other than BMP2 within the Association Region

We investigated the presence of other genes besides BMP2 within the region covered by the disease-associated haplotypes. All predicted exons and ESTs in the region were considered potential exons in yet undiscovered genes. Reverse transcript polymerase chain reaction (RT–PCR) was performed using RNA from a fetal osteoblastic cell line and bone marrow, where the PCR primers were designed so as to connect adjacent exons of a particular predicted open reading frame or ESTs. Several gene-prediction programs were used, as well as comparisons to the syntenic mouse region, which we sequenced. None of the predicted genes gave us products using these sources of RNA. We therefore conclude that the BMP2 gene is the only gene within the region of the haplotype variants.

Discussion

We have mapped osteoporosis susceptibility to the short arm of Chromosome 20 using linkage analysis in Icelandic families. Case–control association analyses using firstly Icelandic samples and subsequently replicated with an independent Danish cohort support BMP2 as the gene, or one of the genes, in the region contributing to osteoporosis risk. We applied a novel definition of the osteoporosis phenotype, based on corrected BMD values at both the spine and the hip in the same analysis (a more generalized osteoporosis) in addition to including osteoporotic fractures. Hitherto, osteoporotic fractures have not been used as a phenotype in any linkage scan reported, neither on their own nor in combination with BMD values. All previously reported linkage studies were based on the use of BMD at specific skeletal sites as a quantitative trait in the analysis. We want to emphasize that we did not consider our approach to be inherently superior, but rather a reasonable alternative that seems to have worked in this instance. We have also observed that pedigrees with low BMD are not a rarity in the population, as has been suggested (Peacock et al. 2002); in fact, the vast majority of probands in this study had close relatives who, when measured, also turned out to have low BMD. However, the skeletal site of low BMD could vary, not necessarily being the same site among the relatives nor among all members within a particular pedigree. This observation is also consistent with our linkage results, applying a combination of both the hip and the spine BMD measurements in the same analysis. Our skeletal site-specific linkage runs both confirm the results on the Chromosome 20p locus and also show that there are other potential loci that seem to be more site specific; however, the LOD scores at these loci did not reach genome-wide significance.

The linkage region on Chromosome 20 contains several genes. Expression analysis alone was not helpful in narrowing the focus to a single gene. Four of the genes in the region are expressed in an osteoblast cell line or in bone marrow (BMP2, C20orf42, C20orf154, and CHGB). However, initial association analysis did draw attention to the BMP2 gene, which is a priori a strong candidate gene, based on its role in bone formation and osteoblast differentiation (Wozney et al. 1988; Katagiri et al. 1994; Fujii et al. 1999). Our analysis indicates that the influence of BMP2 on the development of osteoporosis is that on the attainment on peak bone mass, rather than increased bone loss. The attainment of peak bone mass is a part of the developmental process. The role of BMP2 in this process is probably through stimulating the differentiation and/or activity of osteoblasts. The number of differentiated osteoblasts may be lower, or their BMP2-induced activity lower, in individuals who attain low peak bone mass. The nature of the underlying haplotype variations in the BMP2 gene is not evident, nor is the effect of the Ser37Ala missense variant. The Ser37Ala variant changes an amino acid in the precursor part of the protein, which might have an effect on the stability of the latent dimer or might play a role in the secretion or activation of the protein. Given the diverse role of BMP2 (Wozney et al. 1990; Zhang and Bradley 1996; Kanzler et al. 2000; Takazawa et al. 2000), variations in this gene may have only subtle effects on its activity and most probably in a spatially and temporally tissue-specific manner, through regulation at the level of transcription, alternative splicing, mRNA transport and stability, or translation.

Given the degree of multiple testing in the Icelandic cohort and the fact that the identified at-risk variants cannot fully account for the linkage signal, we could not be confident in the results without confirmation in an independent cohort. Therefore, we evaluated the results we observed in Iceland by testing the haplotype association and the Ser37Ala variant in an independent cohort of Danish osteoporosis patients. The variants have comparable frequencies and estimated risks in this second European population, supporting the hypothesis that the association of BMP2 with osteoporosis is not unique to the Icelandic population. However, the replication sample is only modest in size. Hence, we hope that the promising results presented here will stimulate other replication studies, necessary to further confirm and to provide a more comprehensive evaluation of the role of BMP2 in osteoporosis. Moreover, even if all of our identified at-risk variants were to be confirmed in replication studies, the fact that they cannot fully explain the linkage signal implies that there may be other genes in the region also contributing to the linkage signal or that there are other at-risk variants of BMP2 that have yet to be identified.

The possibility of exploiting genes in the BMP signalling pathway during osteoblast differentiation for new drug targets in the treatment of osteoporosis was pointed out previously (Mundy et al. 2001). Our results lend support to this idea. Variants in BMP2 may also be useful in addressing the unmet medical need of identifying those at increased risk for osteoporosis and fractures.

Materials and Methods

Study subjects and phenotype

The study group was drawn from a large number of individuals who have had BMD measurements at the National University Hospital of Iceland. Individuals with age-matched BMD (Z score) less than 1 SD below average were invited to participate in the study. Spouses, children (18 years and older), parents, and siblings of participating individuals were also invited to participate.

The BMD of all participants, including relatives, was determined using dual-energy X-ray absorptiometry at the lumbar spine (L2–L4) in posterior–anterior projection and at the total hip (proximal end of femur) and whole body (QDR 4500A, Hologic, Waltham, Massachusetts, United States). Weight and height were measured at the time of BMD measurement. All participants completed a detailed questionnaire regarding their medical history, menstrual period, current and past medications, and history of all fractures and trauma. The BMD values were corrected for sex and age to yield the standard Z score, using reference data on 1,407 BMD measurements on randomly collected individuals in Iceland who had not received therapy influencing BMD, thus creating a population-specific reference for corrections. We further corrected these Z scores for weight, which was found to be a significant variable, as has been reviewed by Reid (2002).

Data on osteoporotic fractures are mostly self-reported, collected via a detailed questionnaire, but not further validated. Validation of fractures through medical records (clinical history, X-rays, hospital records) was, however, done for individuals who reported low-impact fractures but had BMD values higher than 1 SD below average. Osteoporotic fractures were defined as low-impact vertebral compression fractures, low-impact hip fractures, or low-impact fractures at other sites if there were two or more such fractures. The lower age limit for osteoporotic fractures was 40 y, and fractures of fingers, hands, toes, or feet were not included. We excluded as affected subjects those with a history of considerable corticosteroid use (7.5 mg per day for 3 mo or more) or with menopause before the age of 40 y, but the excluded subjects could be used as relatives with unknown status. Also, individuals who were suspected to suffer from diseases or deficiencies that influence BMD were also excluded as affected. If available, an early BMD measurement from before the beginning of hormone replacement therapy or the treatment of bisphosphonates was used.

There were 1,323 individuals with BMD measurements, corrected for age, sex, and weight, as indicated above, and regardless of whether the measurements were high or low, who formed the linkage cohort. Pedigrees were identified using our genealogical database of Icelanders (Gulcher and Stefansson 2000) such that all cohort individuals related at five or fewer meioses (five meioses separate relatives who are between first and second cousins) were in the same pedigree. Families that were too large for the linkage software were split for computational purposes. Linkage analysis was performed by defining an osteoporotic phenotype and considering only those meeting the criteria as affected. Thus, for any particular analysis, only a subset of the initial families were informative for linkage. Four pedigree sets, each corresponding to a slightly different phenotype, were analyzed. The phenotype descriptions, the number of families informative for linkage, the number of affected members, and the number of genotyped relatives in the families for these four pedigree sets are displayed in Table 1. In the moderate pedigree set, there were 79 individuals being treated with bisphosphonates and 78 individuals who had suffered an osteoporotic fracture; in the severe pedigree set, there were 70 being treated and 73 with fractures; in the hip pedigree set, there were 40 treated and five with fractures; and in the spine pedigree set, there were 44 treated and five with fractures.

An additional 478 affected individuals were recruited for association analysis, 109 men (mean age, 53.7 y) and 369 women (mean age, 51.9 y). For the analysis of a particular phenotype, all members of the linkage cohort with that phenotype were included along with those having the phenotype in the additional set. The entire association cohort across all phenotypes consisted of 957 individuals, 781 women and 176 men. The control group was comprised of 710 randomly collected individuals from the general population, 468 women and 242 men, in the age range of 30–85 y, whose medical history was unknown. A subset of these groups, both patients and controls, were typed for the Ser37Ala variant, since this required sequencing for typing instead of our routine fluorescent polarization template-directed dye-terminator incorporation (FP–TDI) assay for SNP typing (see below).

All participants gave informed consent, and both the National Bioethics Committee and the Data Protection Commission of Iceland approved the study. Personal identities associated with blood samples and the Decode Genealogy Database were encrypted by the Data Protection Commission, as previously described (Gulcher et al. 2000).

Danish cohort

The study group was selected from the PERF (Prospective Epidemiological Risk Factors) study in Denmark (Bagger et al. 2001). Two osteoporotic patient groups were analyzed: women with BMD in the lowest 10th percentile of the weight-corrected Z score at either the hip or the spine, both at baseline and at follow-up measurement, and a group of postmenopausal osteoporotic fracture patients. The fracture group included women who had suffered two or more low-impact osteoporotic fractures, meeting the same criteria as described in the Icelandic material. The control group consisted of women who had a weight-corrected Z score around the mean (−0.8 to 0.8 SD), both at baseline and at follow-up, randomly picked.

Linkage analysis

A genome-wide scan was performed as previously described (Gretarsdottir et al. 2002), with a framework scan of microsatellite markers. We used multipoint, affected-only allele-sharing methods (Kong and Cox 1997) to assess the evidence for linkage. All results were obtained using the program Allegro (Gudbjartsson et al. 2000) and our genetic map (Kong et al. 2002). We used the Spairs scoring function (Whittemore and Halpern 1994; Kruglyak et al. 1996) and the exponential allele-sharing model (Kong and Cox 1997) to generate the relevant 1 df (degree of freedom) statistics. When combining the family scores to obtain an overall score, instead of weighting the families equally or weighting the affected pairs equally, we used a weighting scheme that is halfway between the two in the log scale; our family weights are the geometric means of the weights of the two schemes. This weighting scheme tends to give results similar to those of the scheme proposed by Weeks and Lange (1988) as an extension of a weighting scheme of Hodge (1984) designed for sibships.

We computed the p value two different ways and report the less significant one. The first p value was computed based on large sample theory: Zlr = √ is approximately distributed as a standard normal random variable under the null hypothesis of no linkage (Kong and Cox 1997). Because of the concern with small sample behavior, we computed a second p value by comparing the observed LOD score to its complete data-sampling distribution under the null hypothesis (Gudbjartsson et al. 2000). When a dataset consists of more than a handful of families, which is the case here, these two p values tend to be very similar. To ensure that the result was a true reflection of the information contained in the material, for us to consider a linkage result significant, not only was it required that the p value be smaller than 2 × 10−5 (Lander and Kruglyak 1995), but also that the information content in the region was at least 85%. The information measure we used has been defined previously (Nicolae 1999) and implemented in Allegro. This measure is closely related to a classical measure of information (Demster et al. 1977), having the property that it is between 0, if the marker genotypes are completely uninformative, and 1, if the genotypes determine the exact amount of allele sharing by descent among the affected relatives.

Adjustment for relatedness of the affected

For the association tests performed for the Ser37Ala variant, if we do not adjust for the familial relationships among some of the affected, the p values are anticonservative, but the estimated RRs are not biased because one can consider the relatedness as having an effect on the variance but not the mean of a test statistic. To address this issue, we extended a variance adjustment procedure described in Risch and Teng (1998) for sibships so that it can be applied to general familial relationships. Specifically, based on the familial relations known about the patients, we first calculated the inflation factor, relative to having unrelated patients, of the variance of the test statistic, defined as the difference between haplotype frequencies in the patients and controls. The p value obtained from the Fisher's exact test, assuming that the patients are unrelated, was converted to a corresponding 1 df χ2 statistic that gives the same p value. This χ2 statistic was then divided by the inflation factor to obtain the adjusted χ2 statistic and the adjusted p value.

Haplotype analysis

To handle missing genotypes and uncertainty with phase, we utilized NEMO (Gretarsdottir et al. 2003), our haplotype analysis program, in which we implemented a likelihood procedure, using the expectation–maximization algorithm as a computational tool, to estimate haplotype frequencies. Under the null hypothesis, the affected individuals and controls are assumed to have identical haplotype frequencies. Under the alternative hypothesis, the candidate at-risk haplotype is allowed to have a higher frequency in affected individuals than controls, while the ratios of the frequencies of all other haplotypes are assumed to be the same in both groups.

Likelihoods are maximized separately under both hypotheses, and a corresponding 1 df likelihood ratio statistic is used to evaluate statistical significance, a method we have used previously (Stefansson et al. 2002; Gretarsdottir et al. 2003). Even though we only searched for haplotypes that increase risk, all reported p values are two-sided for the Icelandic data, but they have not been adjusted for multiple comparisons or for the relationships of some of the Icelandic patients, except for the Ser37Ala variant. However, the p values for the replication study performed with the Danish cohort are one-sided.

LD mapping

In order to study the LD with microsatellites, an extension to the definitions of D′ (Lewontin 1964) and its statistical significance for bi-allelic markers were utilized. The extended D′ is averaged over all the possible allele combinations of the two markers D′, weighted by the marginal allele probabilities (Hedrick 1987). The corresponding p value is defined as the minimum p value for the pair of markers over the same combinations, provided the joint probability is higher than 0.05. Plotting for all marker combinations, D′ in the upper-left corner and the p value in the lower-right corner suggest the LD structure of the region (see Figure 4B).

SNP identification and genotyping

We screened all exons in the entire region of interest (1.7 Mb), both in known genes (BMP2, CHGB, LOC51605, C20orf154, C20orf155, and C20orf42) and in regions where there was an mRNA (AY007089) or a spliced EST (AI971377, BG822004), by direct sequencing from a PCR template. Subsequently, we screened a 59 kb region containing the BMP2 gene by overlapping PCR fragments. All exonic SNPs discovered were genotyped on the entire material, patients, relatives, and controls, using a method for detecting SNPs with fluorescent polarization template-directed dye-terminator incorporation (the SNP–FP–TDI assay) (Chen et al. 1999). The Ser37Ala variant was genotyped by direct PCR sequencing since a functional FP–TDI assay could not be made. Many of the nonexonic SNPs were also genotyped on all the available material.

Screening for other genes than BMP2

The following ESTs were investigated: AY007089, AI971377, BG822004, BF355539, BE145076, AW852841, and AW852981, as well as several predicted genes/exons. Gene prediction programs included GeneScan, MZef, GRAIL, FgSH, HMMGene, Fgene, and Var3. Reverse transcription was performed using Powerscript Reverse Transcriptase (Clontech, Palo Alto, California, United States) and the ThermoScript RT–PCR system (GIBCO–BRL, Carlsbad, California, United States) according to the manufacturers' protocols. Poly(A)+ RNA from bone marrow (Clontech) and total RNA from hFOB 1.19 (a human fetal osteoblastic cell line from the American Type Culture Collection, Manassas, Virginia, United States) were used for cDNA synthesis.

Supporting Information

Accession Numbers

The LocusLink (http://www.ncbi.nlm.nih.gov/LocusLink/) accession numbers for the genes discussed in this paper are BMP2 (LocusLink ID 650), C20orf42 (LocusLink ID 55612), C20orf154 (LocusLink ID 84515), C20orf155 (LocusLink ID 54675), CHGB (LocusLink ID 1114), COL1A1 (LocusLink ID 1277), LOC51605 (LocusLink ID 51605), LRP5 (LocusLink ID 4041), and TCIRG1 (LocusLink ID 10312).