Search tips
Search criteria 


Logo of jbmrJBMR Instructions for AuthorsSubscribe to JBMRAbout JBMRJBMR Editorial BoardJBMR
J Bone Miner Res. 2010 July; 25(7): 1572–1580.
Published online 2010 January 29. doi:  10.1002/jbmr.36
PMCID: PMC3153999

Pathway-Based Genome-Wide Association Analysis Identified the Importance of Regulation-of-Autophagy Pathway for Ultradistal Radius BMD


Wrist fracture is not only one of the most common osteoporotic fractures but also a predictor of future fractures at other sites. Wrist bone mineral density (BMD) is an important determinant of wrist fracture risk, with high heritability. Specific genes underlying wrist BMD variation are largely unknown. Most published genome-wide association studies (GWASs) have focused only on a few top-ranking single-nucleotide polymorphisms (SNPs)/genes and considered each of the identified SNPs/genes independently. To identify biologic pathways important to wrist BMD variation, we used a novel pathway-based analysis approach in our GWAS of wrist ultradistal radius (UD) BMD, examining approximately 500,000 SNPs genome-wide from 984 unrelated whites. A total of 963 biologic pathways/gene sets were analyzed. We identified the regulation-of-autophagy (ROA) pathway that achieved the most significant result (p = .005, qfdr = 0.043, pfwer = 0.016) for association with UD BMD. The ROA pathway also showed significant association with arm BMD in the Framingham Heart Study sample containing 2187 subjects, which further confirmed our findings in the discovery cohort. Earlier studies indicated that during endochondral ossification, autophagy occurs prior to apoptosis of hypertrophic chondrocytes, and it also has been shown that some genes in the ROA pathway (e.g., INFG) may play important roles in osteoblastogenesis or osteoclastogenesis. Our study supports the potential role of the ROA pathway in human wrist BMD variation and osteoporosis. Further functional evaluation of this pathway to determine the mechanism by which it regulates wrist BMD should be pursued to provide new insights into the pathogenesis of wrist osteoporosis. © 2010 American Society for Bone and Mineral Research.

Keywords: osteoporosis, bone mineral density, genome-wide association, regulation of autophagy, whites


Osteoporosis (MIM 166710) is a systemic skeletal disease characterized by low bone mass and microarchitectural deterioration of bone tissue with a consequent increase of bone fragility and susceptibility to fractures.(1) Osteoporosis has become one of the most serious public health problems around the world, causing millions of fractures annually.(2) Fractures at the distal radius of the wrist, also known as Colles' fractures, are the most common fractures in women younger than 75 years of age.(3) Most Colles' fractures are the combined result of a low-energy fall and loss of wrist bone mass owing to osteoporosis. Bone mineral density (BMD) at the distal forearm is the main determinant of susceptibility to Colles' fracture.(4)

BMD, the proxy measure for osteoporosis, has a strong genetic determination. The results of twin and family studies indicate that genetics account for 60% to 90% of the variance in BMD at multiple skeletal sites, including the wrist.(5) BMD variation at different skeletal sites has been shown to share some genetic factors.(6) Since wrist fractures were found to be predictive for fractures at other skeletal sites,(79) the wrist, especially the distal radius, may be a useful site for early assessment of the overall osteoporotic fracture risk.

With the emergence of genome-wide linkage disequilibrium (LD)–based marker panels and the improvements in high-throughput genotyping technology, genome-wide association studies (GWASs) have become feasible. GWASs can systematically survey the whole genome for causal genetic variants for complex traits/diseases and are powerful tools for dissecting the genetic basis for osteoporosis. Several promising candidate genes for osteoporosis or related traits have been identified by GWASs.(1013)

Most published GWASs consider the effects of genetic markers or haplotypes independently and focus only on genes or markers of top-ranking statistical significance. Genes and/or their products often work together, however, interacting in functional groups or pathways to contribute to phenotypic variation or susceptibility to a disease. Therefore, the approach just outlined, which focuses on individual genes, may not be optimally effective for identifying pathophysiologically significant genetic pathways underlying complex traits or diseases such as osteoporosis. Recently, motivated by the availability of gene set enrichment analysis (GSEA) for gene expression data, Wang and colleagues(14) developed a pathway-based genome-wide association (GWA) analysis method that uses up-to-date knowledge of functional interactions among genes in biologic pathways. With this method, multiple genes in a related pathway can be considered jointly, and association with a group of genes as a whole (pathway) can be detected efficiently, even though some individual genes may have only moderate association evidence. The advantage of this approach was illustrated in its application to two Parkinson disease (PD) GWA data sets and one amyotrophic lateral sclerosis (AMD) data set.(14)

To identify novel biologic pathways contributing to wrist osteoporosis and Colles' fractures, we performed a pathway-based GWA analysis on BMD variation at the wrist ultradistal radius (UD) using the method proposed by Wang and colleagues.(14) We identified a significant association, at the pathway-wise level, between UD BMD and the regulation-of-autophagy (ROA) pathway. Our findings support the potential role of the ROA pathway in the pathogenesis of wrist osteoporosis.

Materials and Methods


Discovery GWAS sample: This study was approved by the necessary Institutional Review Boards of all the involved institutions. A total of 984 unrelated US whites of European origin were chosen from our established genetic repertoire currently containing more than 10,000 subjects. Signed informed-consent documents were obtained from all study participants before they entered the study.

People with chronic diseases and conditions that potentially might affect bone mass, structure, or metabolism were excluded from the study. These diseases/conditions included chronic disorders involving vital organs (eg, heart, lung, liver, kidney, and brain), serious metabolic diseases (eg, diabetes, hypo- and hyperparathyroidism, hyperthyroidism, etc.), other skeletal diseases (eg, Paget disease, osteogenesis imperfecta, rheumatoid arthritis, etc.), chronic use of drugs affecting bone metabolism (eg, hormone-replacement therapy, corticosteroid therapy, and anticonvulsant drugs), and malnutrition conditions (eg, chronic diarrhea, chronic ulcerative colitis, etc.). In addition, subjects taking anti-bone-resorptive or bone anabolic agents/drugs such as bisphosphonates were excluded from the study.

For all subjects, measurement of anthropometric variables was performed, and a structured questionnaire addressing lifestyles and medical history was administered. Areal BMD (g/cm2) at the UD region of the wrist was measured by dual-energy X-ray absorptiometry (DXA) with Hologic QDR 4500W densitometers (Hologic, Inc., Bedford, MA, USA). The scanners were calibrated daily, and long-term precision was monitored with external phantoms. The coefficient of variation (CV%) of wrist UD BMD measurements obtained on the Hologic scanner was approximately 2.3%.

FHS replication sample: For replication of our GWAS findings, we used a sample from the Framingham Heart Study (FHS). The FHS is a community-based, multigenerational, longitudinal study of cardiovascular disease and its risk factors, with one substudy focusing on osteoporosis. The phenotype and genotype information of the cohort was downloaded from the Framingham SNP Health Association Resource (SHARe), accessed through the database of Genotypes and Phenotypes (dbGaP) developed by the National Center for Biotechnology Information (NCBI) ( Appropriate procedures have been taken for the use of the data, which include approval from the University of Missouri–Kansas City Institutional Review Board and signatures on the Data Distribution Agreement by all the UMKC investigators who have access to the data.

The FHS comprises the original cohort, offspring, and generation 3 studies. Subjects used as the replication sample in this analysis are those who participated in exam 23 of the first generation and exams 6 and 7 of the second generation with arm BMD as well as other covariates, including age, weight, height, and sex, available. The arm BMD was scanned by DXA (Lunar DPX-L, Lunar Co., Madison, WI, USA). The detailed protocol is available through the NCBI dbGaP. In total, the sample contains 2187 subjects (817 males and 1370 females) from 477 families.


Discovery GWAS sample: Genomic DNA was extracted from whole human blood using a commercial isolation kit (Gentra Systems, Minneapolis, MN, USA) following the protocols detailed in the kit. DNA concentration was assessed by a DU530 UV/VIS spectrophotometer (Beckman Coulter, Inc., Fullerton, CA, USA). Genotyping with the Affymetrix Mapping 250K Nsp and Affymetrix Mapping 250K Sty arrays was performed at the Vanderbilt Microarray Shared Resource at Vanderbilt University Medical Center (Nashville, TN, USA) using the standard protocol recommended by the manufacturer (Affymetrix, Inc., Santa Clara, CA, USA). Briefly, for each array, 250 ng of genomic DNA was digested with either Nsp1 or Sty1 and ligated to adapters that recognize the cohesive 4-base-pair (bp) overhangs. A generic primer that recognizes the ligated adapter sequence was used to amplify the ligation products in a polymerase chain reaction (PCR). The amplified DNA was assayed by agarose gel electrophoresis to verify an average size distribution of 250 to 1500 bp. The amplified DNA was purified per the manufacturer's protocol and quantified by absorbance at 260 and 280 nm. Then 90 µg of purified DNA was fragmented with DNase I and visualized on a 4% agarose gel. Samples with fragment distributions of less than 180 bp were hybridized to the appropriate array (Nsp or Sty). Arrays were stained, washed, and scanned per manufacturer's protocol using immunopure strepavidin (Pierce, Milwaukee, WI, USA), biotinylated antistreptavidin antibody (Vector Labs, Burlingame, CA, USA), and R-phycoerythrin strepavidin (Invitrogen, Carlsbad, CA, USA). Fluorescence intensities were quantitated using an Affymetrix Array Scanner 30007G. Data management and analyses were performed using the Affymetrix GeneChip Operating System. Genotyping calls were determined from the fluorescent intensities using the Dynamic Model-based (DM) algorithm with a 0.33 confidence score setting,(15) as well as the Robust Linear Model with Mahalanobis Distance Classifier with a Bayesian step (B-RLMM) algorithm.(16) According to Affymetrix's guidelines, DM calls were used for quality control. Specifically, subjects with DM call rates of less than 93% were subject to regenotyping. Finally, 99% of all the subjects passed this quality-control criterion, and the average call rate for all the subjects reached greater than 95%. The B-RLMM algorithm was used for the association analyses owing to its improved performance.(16) B-RLMM clustering was performed with 94 samples per cluster.

The final average B-RLMM call rate across the entire sample in this study reached the high level of 98.57%. However, of the initial full set of 500,568 SNPs, we discarded 32,961 SNPs with per-SNP call rates (ie, call rate for each SNP across all samples) of less than 95% and another 36,965 SNPs with allele frequencies deviating from Hardy-Weinberg equilibrium (HWE; p < .001). We further discarded SNPs with a minor allele frequency (MAF) of less than 5% and SNPs that are more than 500 kilobases (kb) away from any gene. The physical distance of 500 kb was set because most enhancers and repressors are less than 500 kb away from genes, and most LD blocks are less than 500 kb. In total, 312,172 SNPs, which covered 14,585 genes, were assayed in the following pathway-based GWA analyses.

FHS replication sample: Genotyping in the FHS sample was carried out with the Affymetrix 500K mapping array plus the Affymetrix 50K supplementary array. For details of the genotyping method, please refer to the Framingham SHARe at the NCBI dbGaP Web site (

Generation of pathway/gene sets collection

We used four different resources to generate a collection of annotated gene sets and pathways to be tested in this study, which are the BioCarta pathway database (, the KEGG pathway database (, the Ambion GeneAssist Pathway Atlas (, and the Gene Ontology (GO) database ( A total of 260, 190, and 380 annotated pathways were retrieved from the BioCarta pathway database, KEGG pathway database, and Ambion GeneAssist Pathway Atlas, respectively. To integrate the GO information into our study, GO annotation files for human genes also were downloaded. We processed the GO annotation files and generated gene sets on the basis of GO level 4 annotations in biologic process and molecular function. Genes whose GO annotations are in level 5 or lower in the GO hierarchy are assigned to their ancestral GO annotations in level 4. In our analysis, we tested only pathways/gene sets in which at least 85% (and in number 10 to 200) of the genes were included in our GWA data set so as to alleviate the multiple-testing problem by avoiding testing too narrowly or too broadly defined functional categories. Overall, 963 pathways/gene sets were analyzed.

Statistical analysis

Analyses in the discovery sample: To detect population stratification that may lead to spurious association results, we used Structure 2.2 software ( and EIGENSOFT 2.0 software ( to investigate the potential substructure/stratification of our sample. The Structure 2.2 program uses a Markov chain Monte Carlo (MCMC) algorithm to cluster individuals into different cryptic subpopulations on the basis of multilocus genotype data.(17) Using the software, we performed independent analyses under three assumed numbers of population strata (k = 2, 3, and 4) using 200 unlinked, randomly selected markers. To confirm the results achieved through Structure 2.2, we further tested population stratification in our sample using EIGENSOFT 2.0 software, which uses both principal-component analysis and a genomic control approach to correct possible statistical bias caused by ancestral differences in whole-genome association studies.(18,19)

Routine GWA statistical analyses on individual SNPs were first conducted by the Wald test implemented in Plink (Version 1.03, to obtain a t statistic for each tested SNP.(20) The parameters, including age, age2, sex, age/age2-by-sex interaction, height, and weight, were tested for their associations with UD BMD. The significant (p ≤ .05) terms (age, age2, sex, and weight) then were included as covariates to adjust the raw UD BMD values for subsequent analyses. The adjusted UD BMD data were normally distributed (p > .1).

Subsequently, pathway-based GWA analysis was performed following the approach developed by Wang and colleagues.(14) Briefly, the main procedures are as follows:

  1. Generation of statistic value of gene-phenotype association: Suppose N genes in the whole genome and each gene with xi (iN) SNPs mapped to it. For a specific gene Gi, the highest SNP-phenotype association statistic value among all xi SNPs was selected to represent the statistic value of the gene, denoted as ri. Mostly, SNPs were mapped to their closest gene. And in rare cases, if an SNP was located within shared regions of two overlapping genes, the SNP was mapped to both genes.
  2. Ranking statistic values of gene-phenotype association: We ranked all genes by sorting their statistic values from the largest to smallest, denoted as gene list L (r1, r2, …, rN).
  3. Enrichment score calculation: For each given pathway/gene set S, composed of NS genes, a statistic enrichment score (ES) was calculated. ES is a weighted Kolmogorov-Smirnov–like running-sum statistic that reflects the overrepresentation of genes within S at the top of the entire ranking list of genes in the genome. The score is calculated by walking down the gene list L, increasing a running-sum statistic when we encounter a gene in S and decreasing it when we encounter genes not in S. ES is the maximum deviation from zero encountered in the random walk.
    equation image
    where An external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu1.jpg, and p(designated as 1 here) is a parameter that gives higher weight to genes with extreme statistic values. Thus, both the rank of each gene in the entire ranked gene list L (NS) and the association statistic value of it (ri) will contribute to the ES value of a pathway.
  4. Permutation and nominal significance assessment: To estimate the significance level of ES, permutation procedures were used to create the null distribution of ES for each pathway/gene set. The phenotype data of sampled subjects were first shuffled. Then the previous three steps were repeated to calculate ES for the pathway/gene set during each permutation. In total, 1000 permutations were done. Finally, a distribution for ES An external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu2.jpg was generated for each pathway/gene set S. The significance of an observed ESS for a pathway/gene set (nominal p value) was estimated as the percentage of permutations whose An external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu3.jpg values were greater than the observed ESS.
  5. Multiple-testing adjustments: To compare the significance among pathways with different numbers of genes, first, a normalizedES, (NES) was constructed based on the observed ES, mean and standard deviation(SD) of An external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu4.jpg.
equation image

Then two measures were used to adjust for multiple-hypothesis testing for more reliable results. False-discovery rate (FDR) is a procedure usually used to control the fraction of expected false-positive findings to stay below a certain threshold. Family-wise error rate (FWER) is another approach to correct for multiple-hypotheses testing that is considered to be a highly conservative procedure seeking to ensure that the list of reported results does not include even a single false-positive gene set.(14,21) For a pathway/gene set, let An external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu5.jpg denote the normalized enrichment score in the observed data. The FRD q value (denoted as An external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu6.jpg) was calculated as the ratio of the fraction of all permutations with NESAn external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu5.jpg to the fraction of observed pathways/gene sets with NESAn external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu5.jpg. FWER p value (denoted as An external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu7.jpg) was calculated as the fraction of pathways/gene sets whose greatest NES among all permutations is greater than An external file that holds a picture, illustration, etc.
Object name is jbmr0025-1572-mu8.jpg. In this study, the significance criteria after multiple-testing correction for a pathway is that both FDR q value and FWER p value are .05 or less.

equation image
equation image

Replication analyses with FHS sample: The statistical analyses for the FHS sample were almost the same as those in the discovery GWAS cohort unless some little modifications regarding that FHS is a family-based study. In brief, first, the raw arm BMD data were adjusted. The parameters, including age, age2, sex, age/age2-by-sex interaction, height, and weight, were tested for their associations with arm BMD. The significant (p ≤ .05) terms (age-by-sex, height, and weight) then were included as covariates to adjust the raw arm BMD values for subsequent analyses. The adjusted UD BMD data were normally distributed (p > .1). Second, the routine GWA statistical analyses on individual SNPs were conducted using FBAT software ( to obtain a Z statistic for each tested SNP. Third, a pathway-based GWA analysis was performed with almost the same method as that in the discovery cohort unless that the Z statistical values for all genes other than the phenotype data were shuffled when doing the permutation and nominal significance assessment so that the family structure of the population could be maintained.(14) A total of 10,000 permutations were done.


Characteristics of study subjects

In the discovery GWAS cohort, a total of 984 individuals were included in the analyses. The basic characteristics of all subjects in this cohort, stratified by sex, are shown in Table 1. The plot of the UD BMD T-scores (mean = −0.53, SD = 2.01) in this population are given in Supplemental Fig. S1. There are 2187 subjects (817 males and 1370 females) from 477 families in the FHS who were studied to replicate the findings from the discovery cohort. The basic characteristics of all subjects in the FHS sample, stratified by sex, are shown in Table 2.

Table 1
Basic Characteristics of the Discovery Cohorta
Table 2
Basic Characteristics of FHS Subjectsa

Pathway-based association analyses in the discovery cohort

No significant population stratification was detected in the discovery cohort (984 subjects) by either Structure 2.2(17) or EIGENSOFT 2.0 software.(18) When using Structure 2.2 software, all subjects were tightly clustered together, suggesting no population stratification (data not shown). The results from EIGENSOFT 2.0 software(18) provided further support for the results from Structure 2.2. Only one principal component was significant in the principal-component analysis (p < .001), suggesting that only one population ancestry existed for our sample. Based on genome-wide SNP information, inflation factor λ, a measure of population stratification, also was estimated by the genomic control(19) approach implemented in EIGENSOFT 2.0 software. Ideally, for a homogeneous population with no stratification, the value of the inflation factor parameter λ should be equal to or near 1. In our study, λ = 1.01, indicating that the sample population is homogeneous.

Among all 963 pathways tested in our analysis, the ROA pathway, which contains 25 genes, showed the strongest significant association with wrist UD BMD. The ROA pathway achieved a high ES of 0.557 with a p value of .005 (Figs. 1 and and2).2). Most of the genes in the ROA pathway ranked at the top in significance in the gene list containing 14,585 genes genome-wide, with 13 genes ranked among the top 2000 genes (Fig. 2). More important, it is the only pathway, among all 963 pathways tested, that achieved statistically significant FDR and FWER values (NES = 2.50, qfdr = 0.043, pfwer = 0.016).

Fig. 1
–Log10 p values of all the 963 pathways with UD BMD. The arrow points to the ROA pathway. It should to be noted that there are three pathways (marked with *) whose nominal p values were less than that of the ROA pathway. However, after multiple-testing ...
Fig. 2
Running-sum plot for the ROA pathway, including the location of the maximum enrichment score (ES) and the leading-edge subset. The x axis is the rank of the 25 genes in the ROA pathway in the whole gene list generated by ranking all the genes by their ...

The ROA pathway is annotated by the KEGG database. It is the only pathway related to autophagy among all the pathways tested. As shown in Table 1, the ROA pathway contains 6 genes (ATG12, ATG3, ATG5, ATG7, BECN1, and GABAPAPL1) encoding essential proteins involved directly in different stages of autophagy, for example, formation and elongation of phagophores and autophagosomes.(23,24) In addition, it contains 14 genes encoding type I and type II interferons (IFNs), which are pleiotropic cytokines that contribute to regulation of autophagy.(25,26) The pathway also contains five genes encoding subunits of kinases, which are autophgy modulators. These kinases include the catalytic subunit and regulatory subunit of lipid kinase class III phosphoinositide 3-kinase (PIK3C3), two catalytic subunits of AMP-activated protein kinase (AMPK), and unc-51-like kinase (ULK2). The individual gene with the most significant p value among all the genes in this pathway is IFNA14 (p = 1.10 × 10−3). Other genes in the ROA pathway that also contribute positively to the ES (ie, the genes that ranked before or at the point of ES, also denoted as leading-edge genes) are ATG12 (p = 2.11 × 10−3), IFNA21 (p = 2.90 × 10−3), IFNA5 (p = 3.16 × 10−3), IFNA17 (p = 3.49 × 10−3), IFNA8 (p = 4.92 × 10−3), IFNA16 (p = 5.46 × 10−3), IFNA4 (p = 6.02 × 10−3), IFNA7 (p = 6.55 × 10−3), ATG7 (p = 8.32 × 10−3), IFNA13 (p = 1.21 × 10−2), PIK3C3 (p = 1.57 × 10−2), and ATG5 (p = 1.62 × 10−2) (Table 3).

Table 3
Genes in Regulation-of-Autophagy (ROA) Pathway

Information for the SNPs whose p values were selected to represent the p values of the 13 leading-edge genes is given in Table 4. All the SNPs are common variants with minor allele frequencies (MAF) of 8% or greater in our sample. We used the FASTSNP program ( to analyze potential functions of these SNPs. The SNP rs7037147 is located in the intron of IFNA14 and is suggested as an intronic enhancer. An A-to-G change at this locus will delete a binding site of the transcription factor SRY. The SNP rs10757219 is located at a transcription factor binding site at the upstream regulatory region of IFNA8. An A-to-G change at this locus will delete a binding site of the transcription factor c-Myc.

Table 4
The Most Significant SNPs Mapped to the Leading-Edge Genes of the ROA Pathway

Replication study in the FHS sample

To further confirm our GWAS findings in the discovery cohort, we carried out a replication analysis with a sample from the FHS population, which contains 2187 subjects. Same as that in the discovery cohort, 963 pathways were analyzed. In total, this analysis covered 15,598 genes genome-wide. All 25 genes of the ROA pathway, as shown in Table 3, were included for association testing. For the ROA pathway, the observed ES is 0.24, observed NES = 13.77; based on 10,000 permutation, nominal p value < 10−4, qfdr = 0.001, and pfwer = 0.007. That is, the ROA pathway showed significant association with arm BMD in this FHS sample even after multiple-testing adjustment. The leading-edge genes of the ROA pathway in the replication FHS cohort are PIK3C3 (p = 1.97 × 10−3), ATG12 (p = 8.16 × 10−3), PRKAA2 (p = 8.28 × 10−3), ATG5 (p = 1.45 × 10−2), GABARAPL1 (p = 3.44 × 10−2), BECN1 (p = 4.45 × 10−2), IFNA13 (p = 4.72 × 10−2), and ATG7 (p = 4.92 × 10−2), five of which (ie, PIK3C3, ATG12, ATG5, IFNA13, and ATG7) are also leading-edge genes in the discovery cohort. Therefore, there is significant overlapping in the genes involved between the discovery and replication cohorts that contribute to the association signals at the ROA pathway level. The information of the most significant SNPs mapped to these leading-edge genes in the replication FHS cohort are presented in Supplemental Table S1.


In this study, a genome-wide association analysis was conducted to identify the pathway(s) involved in BMD determination at a clinically important skeletal site—wrist UD. The ROA pathway was shown to be significantly associated with wrist UD BMD (p = .005) and was the only pathway, among a total of 963 pathways tested, that achieved statistically significant FDR and FWER values. To further confirm our findings in the discovery cohort, 2187 subjects from the FHS with arm BMD data available were studied as a replication study. The ROA pathway showed significant association with arm BMD in the FHS sample even after the multiple-testing adjustment, which further confirmed our findings in the discovery cohort. This is the first time that the autophagy-related biologic process has been implicated as an underlying factor for wrist BMD variation and hence a risk factor for wrist osteoporosis.

Autophagy is a conserved housekeeping function of eukaryotic cells that permits sequestration of cytoplasmic components in membrane-bound vesicles and delivers them to lysosomes for degradation. Autophagy is subject to suppression or further induction in response to different stresses, starvation, specific hormonal regulation, and other stimuli.(24) The core process of autophagy and the functional interactions among the genes in the ROA pathway are depicted in Fig. 3. The autophagy pathway has numerous proposed physiologic functions. Growing evidence suggests that malfunctions of autophagy contribute to many human diseases, for example, myopathies, liver disease, neurodegeneration, heart disease, and cancer.(23,24)

Fig. 3
Autophagy and the functional interactions among the genes in the ROA pathway. Basically, the autophagic process can be divided into three stages: vesicle nucleation, vesicle elongation and completion, and vesicle breakdown and degradation. ATG1, ULK, ...

Although this article is the first to support the ROA pathway's importance to osteoporosis, the pathway's relevance to bone metabolism has been suggested in previous studies. Autophagy has been recognized as an intermediate stage during the terminal differentiation of chondrocytes. Specifically, during endochondral ossification, induction of autophagy occurs prior to apoptotic hypertrophic maturation in terminally differentiated chondrocytes.(27,28) These chondrocytes morphologically exhibit autophagic characteristics and express autophagic proteins.(27,29) The induction of autophagy enhances the survival of chondrocytes in hypoxic microenvironments (eg, low protein, glucose, and O2 concentrations) and facilitates synthesis of the calcified extracellular matrix.(3032)

Some genes in the ROA pathway also have been considered as important modulating factors for bone development or remodeling. For example, IFN-α can suppress the proliferation of osteoprogenitor cells(33,34) and modify the expression of a number of important cytokines that are regulators for human osteoprogenitor cell growth and differentiation.(3537) IFN- α also can inhibit the differentiation of peripheral blood mononuclear cells (PBMCs) to osteoclasts.(38) Interferon-γ (IFN-γ) is considered to be an antagonist of RANKL and is a well-known potent inhibitor of osteoclast function and formation.(3942) IFN-γ also exhibits antiproliferative actions on primary osteoblast cells.(4345) Finally, it has been reported that in cultured osteoblasts, AMPK influences the expression of cyclooxygenase 2 (COX-2), which potentially could impact fracture healing because COX-2 is a crucial mediator in mechanically induced bone formation.(46)

In addition to highlighting the potential contribution of the ROA pathway to human wrist BMD deviation, this study also demonstrates the important advantage of pathway-based GWA analysis. None of the genes in the ROA pathway reached a significant level when considered separately; thus the potential impact of the ROA pathway on wrist BMD would not have been detected by individual gene/SNP GWA analyses. By considering critical information about the interaction of a set of functionally related genes and their joint effects, pathway analysis may be a useful paradigm for revealing the polygenic nature of complex diseases. Recently, some candidate pathway association studies have been done, suggesting that the joint actions of common gene variants within pathways may play a major role in predisposing to complex diseases. Examples include the axon guidance pathways for Parkinson disease (PD)(47) and amyotrophic lateral sclerosis (ALS).(48) Differing from the candidate pathway strategies,(47,48) in the current genome-wide analysis we systematically studied a large number of pathways annotated by public pathway databases with the aim of identifying pathways associated with UD BMD. This hypothesis-free strategy could make full use of GWA data and would be more useful for discovering new disease-related processes and generating novel hypotheses about pathogenic mechanisms of disease.

It would be interesting to clarify the detailed mechanism of the ROA pathway on wrist BMD. Considering the existing evidence linking some ROA pathway genes with osteoblastogenesis or osteoclastogenesis, it can be imagined that the involvement of autophagy in bone may not be limited to chrondrocytes alone. Since autophagy is another mechanism for programmed cell death and shows close interaction with apoptosis, whether antophagy is involved in the survival control of osteoblasts/osteoclasts/osteocytes is an important hypothesis to test. Also, it will be interesting to study whether some lifestyle factors will interact with the ROA pathway or if the ROA pathway would be a candidate target for osteoporosis therapy.

In this study we found that the ROA pathway is associated with human wrist BMD variation. However, we cannot exclude other pathways for their significance in wrist BMD variation, even though these pathways did not pass the significance criteria adopted in this study. Regarding the statistical aspect of this study, it should be noted that some genes inevitably are shared by different pathways. Although the overlap of genes among different pathways will not affect the pathways' relative ranking in terms of NES values, interdependence between pathways will lead to a decreased power by affecting FDR and FWER when the causal genes are shared by multiple pathways (owing to the permutation procedure for multiple-testing adjustment). A larger sample size or more specifically annotated pathways will be helpful to improve the power of this pathway-based GWAS.

It also should be noted that under the significance criteria adopted in this study, the ROA pathway was not detected to be associated with hip or spine BMD in the discovery cohort (data not shown). The possible reasons for the inconsistent association signals detected at these three skeletal sites may include statistical errors (type I error at the wrist or type II error at the hip/spine) and genetic heterogeneity between different skeletal sites. In particular, type II error at the hip/spine may be the major reason. The type II error may be caused by the relatively low statistical power of this newly developed pathway-based method (as discussed earlier). This low statistical power, when complicated by the possible genetic heterogeneity (ie, smaller effects of the ROA pathway on hip and spine than on wrist), may lead to negative association signals at the hip and spine.

In summary, we applied a novel pathway-based GWA analysis method to systematically screen functional pathways/gene sets associated with wrist UD BMD and thus osteoporosis. The significant enrichment of ROA pathway genes among the top-ranking genes associated with wrist UD BMD, together with the pathway's functional relevance to bone metabolism, strongly supports an important role of autophagy in human wrist BMD variation. Further detailed and specific functional studies of the ROA pathway will provide new insights into the pathway's relevance to the physiology of bone and the etiology of wrist osteoporosis.


Investigators of this work were supported in part by grants from the NIH (R01 AR050496, R21 AG027110, R01 AG026564, P50 AR055081, and R21 AA015973). The study also benefited from grants from the National Science Foundation of China, the Huo Ying Dong Education Foundation, HuNan Province, Xi'an Jiaotong University, and the Ministry of Education of China. We also thank all contributors to the Framingham Heart Study. The Framingham Heart Study and the Framingham SHARe project are conducted and supported by the National Heart, Lung, and Blood Institute (NHLBI) in collaboration with Boston University. The Framingham SHARe data used for the analyses described in this article were obtained through dbGaP. This article was not prepared in collaboration with investigators of the Framingham Heart Study and does not necessarily reflect the opinions or views of the Framingham Heart Study, Boston University, or the NHLBI.


Dr. Recker has been paid consultant for Merck, Lilly, Wyeth, Procter&Gamble, Amgen, Roche, Glaxo Smith Kline, Novartis, and NPS Allelix; and has received grant/research support from Merck, Lilly, Wyeth, Procter&Gamble, Amgen, Roche, Glaxo Smith Kline, Novartis, NPS Allelix, and Sanofi-Aventis. All other authors state that they have no conflicts of interest.

Supplementary material


1. Johnell O, Kanis J. Epidemiology of osteoporotic fractures. Osteoporos Int. 2005;16:S3–7. [PubMed]
2. Chao EY, Inoue N, Koo TK, Kim YH. Biomechanical considerations of fracture treatment and bone quality maintenance in elderly patients and patients with osteoporosis. Clin Orthop Relat Res. 2004;425:12–25. [PubMed]
3. Muller ME, Webber CE, Bouxsein ML. Predicting the failure load of the distal radius. Osteoporos Int. 2003;14:345–352. [PubMed]
4. Weichetova M, Stepan JJ, Haas T, Michalska D. The risk of Colles' fracture is associated with the collagen I α1 Sp1 polymorphism and ultrasound transmission velocity in the calcaneus only in heavier postmenopausal women. Calcif Tissue Int. 2005;76:98–106. [PubMed]
5. Tse KY, Macias BR, Meyer RS, Hargens AR. Heritability of bone density: regional and gender differences in monozygotic twins. J Orthop Res. 2009;27:150–154. [PubMed]
6. Duncan EL, Cardon LR, Sinsheimer JS, Wass JA, Brown MA. Site and gender specificity of inheritance of bone mineral density. J Bone Miner Res. 2003;18:1531–1538. [PubMed]
7. Cuddihy MT, Gabriel SE, Crowson CS, O'Fallon WM, Melton LJ., 3rd. Forearm fractures as predictors of subsequent osteoporotic fractures. Osteoporos Int. 1999;9:469–475. [PubMed]
8. Gardsell P, Johnell O, Nilsson BE, Gullberg B. Predicting various fragility fractures in women by forearm bone densitometry: a follow-up study. Calcif Tissue Int. 1993;52:348–353. [PubMed]
9. Mallmin H, Ljunghall S, Persson I, Naessen T, Krusemo UB, Bergstrom R. Fracture of the distal forearm as a forecaster of subsequent hip fracture: a population-based cohort study with 24 years of follow-up. Calcif Tissue Int. 1993;52:269–272. [PubMed]
10. Styrkarsdottir U, Halldorsson BV, Gretarsdottir S, et al. Multiple genetic loci for bone mineral density and fractures. N Engl J Med. 2008;358:2355–2365. [PubMed]
11. Richards JB, Rivadeneira F, Inouye M, et al. Bone mineral density, osteoporosis, and osteoporotic fractures: a genome-wide association study. Lancet. 2008;371:1505–1512. [PMC free article] [PubMed]
12. Liu YZ, Wilson SG, Wang L, et al. Identification of PLCL1 gene for hip bone size variation in females in a genome-wide association study. PLoS ONE. 2008;3:e3160. [PMC free article] [PubMed]
13. Kiel DP, Demissie S, Dupuis J, Lunetta KL, Murabito JM, Karasik D. Genome-wide association with bone mass and geometry in the Framingham Heart Study. BMC Med Genet. 2007;8(Suppl 1):S14. [PMC free article] [PubMed]
14. Wang K, Li M, Bucan M. Pathway-Based Approaches for Analysis of Genomewide Association Studies. Am J Hum Genet. 2007;81:1278–1283. [PubMed]
15. Di X, Matsuzaki H, Webster TA, et al. Dynamic model based algorithms for screening and genotyping over 100K SNPs on oligonucleotide microarrays. Bioinformatics. 2005;21:1958–1963. [PubMed]
16. Rabbee N, Speed TP. A genotype calling algorithm for affymetrix SNP arrays. Bioinformatics. 2006;22:7–12. [PubMed]
17. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–959. [PubMed]
18. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet. 2006;38:904–909. [PubMed]
19. Devlin B, Roeder K. Genomic control for association studies. Biometrics. 1999;55:997–1004. [PubMed]
20. Purcell S, Neale B, Todd-Brown K, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–575. [PubMed]
21. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–15550. [PubMed]
22. Horvath S, Xu X, Laird NM. The family based association test method: strategies for studying general genotype--phenotype associations. Eur J Hum Genet. 2001;9:301–306. [PubMed]
23. Mizushima N, Levine B, Cuervo AM, Klionsky DJ. Autophagy fights disease through cellular self-digestion. Nature. 2008;451:1069–1075. [PMC free article] [PubMed]
24. Espert L, Codogno P, Biard-Piechaczyk M. Involvement of autophagy in viral infections: antiviral function and subversion by viruses. J Mol Med. 2007;85:811–823. [PubMed]
25. Talloczy Z, Jiang W, Virgin HWt, et al. Regulation of starvation- and virus-induced autophagy by the eIF2alpha kinase signaling pathway. Proc Natl Acad Sci U S A. 2002;99:190–195. [PubMed]
26. Pyo JO, Jang MH, Kwon YK, et al. Essential roles of Atg5 and FADD in autophagic cell death: dissection of autophagic cell death into vacuole formation and cell death. J Biol Chem. 2005;280:20722–20729. [PubMed]
27. Shapiro IM, Adams CS, Freeman T, Srinivas V. Fate of the hypertrophic chondrocyte: microenvironmental perspectives on apoptosis and survival in the epiphyseal growth plate. Birth Defects Res C Embryo Today. 2005;75:330–339. [PubMed]
28. Bohensky J, Shapiro IM, Leshinsky S, Watanabe H, Srinivas V. PIM-2 is an independent regulator of chondrocyte survival and autophagy in the epiphyseal growth plate. J Cell Physiol. 2007;213:246–251. [PubMed]
29. Zuscik MJ, Hilton MJ, Zhang X, Chen D, O'Keefe RJ. Regulation of chondrogenesis and chondrocyte differentiation by stress. J Clin Invest. 2008;118:429–438. [PMC free article] [PubMed]
30. Roach HI, Aigner T, Kouri JB. Chondroptosis: a variant of apoptotic cell death in chondrocytes? Apoptosis. 2004;9:265–277. [PubMed]
31. Srinivas V, Shapiro IM. Chondrocytes embedded in the epiphyseal growth plates of long bones undergo autophagy prior to the induction of osteogenesis. Autophagy. 2006;2:215–216. [PubMed]
32. Bohensky J, Shapiro IM, Leshinsky S, Terkhorn SP, Adams CS, Srinivas V. HIF-1 regulation of chondrocyte apoptosis: induction of the autophagic pathway. Autophagy. 2007;3:207–214. [PubMed]
33. Broxmeyer HE, Lu L, Platzer E, Feit C, Juliano L, Rubin BY. Comparative analysis of the influences of human gamma, alpha and beta interferons on human multipotential (CFU-GEMM), erythroid (BFU-E) and granulocyte-macrophage (CFU-GM) progenitor cells. J Immunol. 1983;131:1300–1305. [PubMed]
34. Oreffo RO, Romberg S, Virdi AS, Joyner CJ, Berven S, Triffitt JT. Effects of interferon alpha on human osteoprogenitor cell growth and differentiation in vitro. J Cell Biochem. 1999;74:372–385. [PubMed]
35. Aman MJ, Keller U, Derigs G, Mohamadzadeh M, Huber C, Peschel C. Regulation of cytokine expression by interferon-alpha in human bone marrow stromal cells: inhibition of hematopoietic growth factors and induction of interleukin-1 receptor antagonist. Blood. 1994;84:4142–4150. [PubMed]
36. Aman MJ, Bug G, Aulitzky WE, Huber C, Peschel C. Inhibition of interleukin-11 by interferon-alpha in human bone marrow stromal cells. Exp Hematol. 1996;24:863–867. [PubMed]
37. Gollner G, Aman MJ, Steffens HP, Huber C, Peschel C, Derigs HG. Interferon-alpha (IFN-alpha) inhibits granulocyte-macrophage colony-stimulating factor (GM-CSF) expression at the post-transcriptional level in murine bone marrow stromal cells. Br J Haematol. 1995;91:8–14. [PubMed]
38. Avnet S, Cenni E, Perut F, et al. Interferon-alpha inhibits in vitro osteoclast differentiation and renal cell carcinoma-induced angiogenesis. Int J Oncol. 2007;30:469–476. [PubMed]
39. Nii A, Reynolds DA, Young HA, Ward JM. Osteochondrodysplasia occurring in transgenic mice expressing interferon-gamma. Vet Pathol. 1997;34:431–441. [PubMed]
40. Mann GN, Jacobs TW, Buchinsky FJ, et al. Interferon-gamma causes loss of bone volume in vivo and fails to ameliorate cyclosporin A-induced osteopenia. Endocrinology. 1994;135:1077–1083. [PubMed]
41. Gowen M, Nedwin GE, Mundy GR. Preferential inhibition of cytokine-stimulated bone resorption by recombinant interferon gamma. J Bone Miner Res. 1986;1:469–474. [PubMed]
42. Kurihara N, Roodman GD. Interferons-alpha and -gamma inhibit interleukin-1 beta-stimulated osteoclast-like cell formation in long-term human marrow cultures. J Interferon Res. 1990;10:541–547. [PubMed]
43. Cornish J, Gillespie MT, Callon KE, Horwood NJ, Moseley JM, Reid IR. Interleukin-18 is a novel mitogen of osteogenic and chondrogenic cells. Endocrinology. 2003;144:1194–1201. [PubMed]
44. Nakano T, Otsuka T, Ogo T, et al. Transfection of interferon-gamma gene in a mouse bone marrow stromal preadipocyte cell line causes apoptotic cell death. Exp Hematol. 1993;21:1498–1503. [PubMed]
45. Yoshihara R, Shiozawa S, Imai Y, Fujita T. Tumor necrosis factor alpha and interferon gamma inhibit proliferation and alkaline phosphatase activity of human osteoblastic SaOS-2 cell line. Lymphokine Res. 1990;9:59–66. [PubMed]
46. Hou CH, Tan TW, Tang CH. AMP-activated protein kinase is involved in COX-2 expression in response to ultrasound in cultured osteoblasts. Cell Signal. 2008;20:978–988. [PubMed]
47. Lesnick TG, Papapetropoulos S, Mash DC, et al. A genomic pathway approach to a complex disease: axon guidance and Parkinson disease. PLoS Genet. 2007;3:e98. [PubMed]
48. Lesnick TG, Sorenson EJ, et al. Beyond Parkinson disease: amyotrophic lateral sclerosis and the axon guidance pathway. PLoS ONE. 2008;3:e1449. [PMC free article] [PubMed]

Articles from Journal of Bone and Mineral Research are provided here courtesy of American Society for Bone and Mineral Research