Pathway-Based Association Analyses Identified TRAIL Pathway for Osteoporotic Fractures

Introduction Hip OF carries the highest morbidity and mortality. Previous studies revealed that individual genes/loci in the Tumor Necrosis Factor (TNF) -Related Apoptosis-Inducing Ligand (TRAIL) pathway were associated with bone metabolism. This study aims to verify the potential association between hip OF and TRAIL pathway. Methods Using genome-wide genotype data from Affymetrix 500 K SNP arrays, we performed novel pathway-based association analyses for hip OF in 700 elderly Chinese Han subjects (350 with hip OF and 350 healthy matched controls). Results The TRAIL pathway achieved a significant p value (p = 0.01) for association with hip OF. Among the 38 genes in the TRAIL pathway, seven genes achieved nominally significant association with hip OF (p<0.05); the TNFSF10 (TRAIL) gene obtained the most significant p value (p = 1.70×10−4). SNPs (rs719126, rs6533015, rs9594738, rs1805034, rs11160706) from five genes (CFLAR, NFKB1, TNFSF11, TNFRSF11A, TRAF3) of the pathway had minor alleles that appear to be protective to hip OF. SNPs (rs6445063 and rs4259415) from two genes (TNFSF10 and TNFRSF10B) of the pathway had minor alleles (A) that are associated with an increased risk of hip OF, with the ORs (odds ratios) of 16.51 (95%CI:3.83–71.24) and 1.37 (95%CI:1.08–1.74), respectively. Conclusions Our study supports the potential role of the TRAIL pathway in the pathogenesis of hip OF in Chinese Han population. Further functional study of this pathway will be pursued to determine the mechanism by which it confers risk to hip OF.


Introduction
The prevalence of osteoporotic fractures (OF), especially hip OF, has now become a significant public health burden due to the associated high morbidity, mortality and tremendous health care cost [1,2] The incidence of hip OF will reach 6.27 million worldwide with a resultant cost of .$34 billion by the year 2050 [1]. The major demographic change of hip OF will occur in Asia. In 1990, 26% of hip OF worldwide occurred in Asia, and this rate could rise to 37% by the year 2025 and to 45% by the year 2050 [2]. Genetic factors play an important role in susceptibility to hip OF [3,4]. However, to date, the genetic determination of hip OF is still largely unknown.
A pathway-based association analysis approach is based on gene set enrichment analysis. The method ranks genes/SNPs genome-wide by the significance of association with a disease phenotype to generate a gene list, then uses a statistic enrichment score (ES) to examine whether a particular group of genes/SNPs within a certain functional pathway is enriched at the top (or bottom) of the list more than would be expected by chance [5]. This approach considers critical information about the interaction of a set of functionally related genes and their joint effects, and could potentially improve the power to detect genetic variants working in functional pathways. The method may play a major role in genetic analyses of complex diseases; for a pathway that contributes to a disease's risk, a single genetic variant within that pathway may have only limited contribution to the risk that might not be detectable (in regular genetic association analyses) if not otherwise be detected at the whole pathway level. Recently, through this new powerful approach, we have identified two novel pathways for wrist BMD and femur geometry in US whites [6,7].
Previous studies revealed importance in the pathogenesis of bone metabolism for several individual genes/loci, or their transcripts and proteins in Tumor Necrosis Factor (TNF)-Related Apoptosis-Inducing Ligand (TRAIL) pathway, for example, TRAIL, TRAILR1, TRAILR2, CFLAR, RANK, and OPG [8,9,10]. However, it is unknown if the TRAIL pathway may contribute to osteoporotic fracture risk at the whole pathway level. We therefore undertook a novel, pathway-based association study to assess contribution of the pathway to hip OF risk.

Subjects
The study was approved by the institutional review board of Xi'an Jiaotong University. After signing an informed consent agreement, all subjects were assisted to complete a structured questionnaire including anthropometric variables, lifestyles, and medical history.
The sample consisted of 700 elderly Chinese Han subjects, including 350 with osteoporotic (low trauma) hip fractures (including 124 males and 226 females) and 350 controls (including 173 males and 177 females) (see Table 1 for the basic characteristics of these subjects). All the subjects were unrelated Chinese Han adults living in the city of Xi'an and its neighboring areas. Inclusion criteria for cases were (i) onset age .55 years, all female subjects were postmenopausal women; (ii) age ,80 years to minimize the effect due to age, since previous studies showed that approximately half of females aged 80 years or older have fractures [11]; (iii) minimal or no trauma fractures, usually due to falls from standing height or less; (iv) fracture site at femoral neck or intertrochanter regions; (v) hip fracture was identified/confirmed through diagnosis of orthopedic surgeons/radiologists according to radiological reports and X-rays. Patients with pathological fractures and high-impact fractures (such as due to motor vehicle accidents) were excluded.
Healthy control subjects were selected from our established large database. They were geographically matched to the cases. Inclusion criteria for controls were: (i) age at exam .55 years, without any fracture histories (all female controls were postmenopausal); and (ii) not suffering from chronic diseases and conditions that might potentially affect bone mass, structure, or metabolism. The criteria may ensure that controls are less likely to suffer OF during their remaining life span compared with the general populations. The diseases/conditions mentioned above included chronic disorders involving vital organs (heart, lung, liver, kidney, brain), serious metabolic diseases (diabetes, hypo-and hyperparathyroidism, hyperthyroidism, etc.), other skeletal diseases (Paget disease, osteogenesis imperfecta, rheumatoid arthritis, etc.); chronic use of drugs affecting bone metabolism (e.g., hormone replacement therapy, corticosteroid therapy, anti-convulsant drugs), and malnutrition conditions (such as chronic diarrhea, chronic ulcerative colitis). In addition, subjects taking anti-resorptive or bone anabolic agents/drugs, such as bisphosphonates were also excluded.

Genotyping and quality control
Genomic DNA was extracted from peripheral blood leukocytes using a commercial isolation kit (Gentra systems, Minneapolis, MN, USA) following the standard protocols. SNP genotyping was performed using Affymetrix Human Mapping 500 K array set (Affymetrix, Santa Clara, CA, USA). Genotyping calls were determined from the fluorescent intensities using the dynamic model (DM) algorithm [12] as well as the B-RLMM (Bayesian Robust Linear Model with Mahalanobis Distance Classifier) algorithm [13]. DM calls were used for quality control of the genotyping experiment. BRLMM calls were used for the following pathway-based analyses. SNPs with a sample call rate ,90%, with allele frequencies extremely deviating from Hardy-Weinberg equilibrium (p,10 27 ) and having a minor allele frequency (MAF) ,0.01 in the total sample were discarded. The final SNP number for the analyses is 371,258. This SNP set covers 14,642 genes and yielded an average marker spacing of ,7.9 kb throughout the human genome. Among these SNPs, 350 belong to the genes of the TRAIL pathway.

Statistical Analysis
Statistical analyses on individual SNPs were first conducted by the Wald test implemented in Plink (version 1.03) [14] with age and sex as covariates. The main procedures of pathway-based analysis [15] are briefly summarized as follows: (1) Generation of gene-phenotype association rank and calculation of ES: Among all the SNPs of a given gene G i , the SNP achieving the smallest p value in the single-marker association tests was used to represent the magnitude of association evidence of the gene. We ranked all genes by sorting their statistic values from the largest to smallest, denoted as a vector of the gene list L (r 1 ,r 2 ,…,r N ). ''N'' equals to the total gene number contained in the GeneChipH Human Mapping 500 K set arrays. ''r'' represents the gene-phenotype association statistic value. Then, a weighted Kolmogorov-Smirnovlike running-sum enrichment score (ES S ) was calculated in Equation 1. For the given pathway S composed of N S genes, by walking down the gene list L, we increased the runningsum statistic for the pathway when we encountered a gene in the S and decreased it when we encountered a gene not in the S.
Where N R~P G j Ã [S jr i Ã j and parameter p (designated as 1 here) is designed to give higher weight to genes with extreme statistic values. The magnitude of the increment depends on the correlation of the gene with the phenotype. In short, ES S is the maximum deviation from zero encountered in the random walk. It will be high if the association signal is enriched at the top of list L, reflected by the significance level of observed ES S (i.e. nominal p value). (2) Permutation and nominal significance assessment: The phenotype data was shuffled and permutation (per) was done to compute a ES S per through repeating steps (1),(2) to estimate the nominal p value. A total of 1,000 permutations were done to create a histogram of the corresponding enrichment scores ES S per,null for the given pathway/gene set S. The nominal p value was estimated as the percentage of permutations whose ES S per were greater than the observed ES S .
To detect possible population stratification that may lead to spurious association results, we used Structure 2.2 software (http://pritch.bsd.uchicago.edu/software.html) to investigate the potential population substructure/stratification of the sample. The program uses a Markov chain Monte Carlo (MCMC) algorithm to cluster individuals into different cryptic sub-populations on the basis of multi-locus genotype data [16]. We performed the analysis assuming the number of population strata k = 2 and using 1,000 un-linked SNPs randomly selected genome-wide. To confirm the results achieved through Structure 2.2, we further tested population stratification in our sample using EIGENSOFT 2.0 software that uses both principal component analysis and a genomic control approach to correct possible statistical bias caused by ancestral differences [17,18]. Reported p value was adjusted by the l value estimated by genomic control method.

Characteristics of study subjects
Basic characteristics of the 700 subjects are presented in Table 1. The STRUCTURE program revealed that all subjects in this Chinese sample were clustered together and could not be assigned into any subgroups, indicating that there was no significant population stratification within the sample. We also performed analysis using EIGENSOFT [17] and confirmed the homogeneity of the sample revealed by Structure; there is only one principal component that is significant in the principal component analysis (p,0.001), suggesting only one population ancestry existing for our sample. Finally, we used the genomic control method [18] to estimate the l value, which equals to 1.02, suggesting again there is no population stratification in the sample.

Pathway-based association analysis
The TRAIL pathway, annotated by the Ambion GeneAssist Pathway Atlas (Table 2), achieved a high NES of 1.76 with a p value of 0.01. Among the 38 genes in the TRAIL pathway, the gene with the most significant p value is TNFSF10 (TRAIL, p = 1.70610 24 ) ( Table 3). Another six genes in the TRAIL pathway, which also contribute positively to the ES (i.e., the genes that ranked before or at the point of ES, also denoted as ''leading edge genes''), are: CFLAR (c-FLIP, FLIP, p = 3.17610 23 ), TNFSF11 (RANKL, p = 6.06610 23 ), TNFRSF11A (RANK, p = 7.63610 23 ), TNFRSF10B (TRAILR2, p = 1.07610 22 ), TRAF3 (p = 1.67610 22 ), and NFKB1 (NF-kB1, p = 2.04610 22 ). Among the 38 genes in the TRAIL pathway, 4 genes, NFKBIE, IKBKG, RELB, DIABLO, were not covered by Affymetrix 500 k Array Set, and the SNPs of the XIAP and APAF1 genes were discarded for failure of passing the quality control criteria (Table 3).
For the top SNP of each gene, the effect size of the minor allele is expressed by odds ratio (OR) derived from logistic regression analyses. The OR value is interpreted in the standard manner. A value of 1.0 indicates no effect. A value greater than 1.0 indicates that the minor allele is associated with an increased hip OF risk, and a value less than 1.0 implies that the minor allele is associated with a decreased risk, hence may be protective. Among the 7 SNPs representing the 7 leading edge genes, 5 SNPs had the protective  Table 3).

Discussion
In this study, the TRAIL pathway was shown to be significantly associated with hip OF (p = 0.01) according to the pathway-based association analysis. This is the first time that the TRAIL pathway is implicated as an underlying factor for hip OF in Chinese. The pathway may contribute to bone mass variation via regulating osteoclast metabolism; previous findings have shown that the TRAIL pathway play roles in modulating the differentiation, function, survival and/or apoptosis of osteoclasts [19,20], which may in turn accelerate bone resorption and consequently the susceptibility to hip OF. The core process of TRAIL and the functional interactions among the genes in the TRAIL pathway are depicted in Figure 1.
TRAIL, also known as Apo2L, is a widely recognized member of the tumor necrosis factor (TNF) ligand family [21,22]. TRAIL is a typical type II membrane protein, binds to aggregates type I transmembrane receptors with cytoplasmic death domains (DD), which ultimately activate a protease cascade leading to apoptosis [23]. Interacting with five different receptors of the TNF receptor superfamily, TRAILR1, TRAILR2, TRAILR3, TRAILR4 and OPG, the ligand is expressed constitutively in a wide range of tissues [22,24,25], including bone-related cells, such as osteoclast or osteoclast precursors. The binding of TRAIL to TRAILR2 (DR5), which contains a conserved death domain (DD), may induce the apoptosis of osteoclast [19]. This may play an important role in the etiology of hip OF. TRAIL induces apoptosis in human differentiated osteoclasts by means of TRAILR2, and activates an intracellular pathway involving caspase-8 and Bid cleavage in the active forms [19]. However, in this study, SNP rs6445063 from TNFSF10 (TRAIL) gene and SNP rs4259415 from TNFRSF10B (TRAILR2) gene of the pathway had minor alleles (A) that are associated with an increased risk of hip OF, with the ORs (odds ratios) of 16.51 (95%CI:3.83-71.24) and 1.37 (95%CI:1.08-1.74), respectively. SNP rs6445063 is located on chromosome 3q26.31 and has never been studied previously. SNP rs4259415 is located on chromosome 8p21.3 and functions as intronic enhancer. According to the FASTSNP program (http://fastsnp. ibms.sinica.edu.tw), a change of ''GRA'' at rs4259415 may lead to alter the binding sites for transcription factor from GATA-1 to CdxA and TATA, which may change the function of apoptosis related TRAILR2 gene and accordingly alter the apoptosis of osteoclast. Furthermore, the apoptosis induced by death receptors can be modulated at several levels. Intracellular antiapoptotic molecules can block the apoptotic signaling pathway. Soluble receptors inhibit TRAIL-mediated apoptosis in part by increasing baseline levels of CFLAR (c-FLIP), which competes with caspase-8 for binding to FADD, and inhibits active caspases [26]. Increased c-FLIP levels following exposure to soluble factors may decrease the apoptosis of osteoclast. Intracellular anti-apoptotic molecules can also divert the apoptosis signaling into alternative responses, i.e., the proliferation/survival of osteoclast via activating NF-kB [20]. Although lack of bone phenotype in mice with gene deletion of TRAIL was observed by Sedger et al [27], Yen et al [28] showed that TRAIL can induce osteoclast formation via direct engagement with the TRAIL death receptor. Yen et al [28,29] suggested that TRAIL-induced osteoclastogenesisis was dependent on activation of NF-kB, ERK, and p38 MAP kinase and TRAF6 dependent. However TRAF6 was not found to contribute positively to the ES as TRAF3 in this study, suggests that there may be other factors that also can function in this capacity. As the soluble decoy receptor of TRAIL, OPG was originally identified as a decoy receptor for RANKL [30] later found to be able to bind TRAIL [31]. OPG acts as antagonist for osteoclast apoptosis induced by TRAIL. But when OPG binds to the RANKL on the surface of osteoblast/bone matrix, it prevents RANKL from binding to its receptor, RANK, which inhibits the formation, activation and survival of multinucleated osteoclasts. In this process, the RANKL/OPG ratio is an important determinant of bone mass and skeletal integrity [32]. Up-regulation of the RANKL gene increases the RANKL/OPG ratio and enhances bone loss  [32]. Our study found a significant association of SNPs (rs9594738 and rs1805034) of RANKL and RANK with hip OF. NF-kB1 (p50) is a subunit of NF-kB transcriptional factor complex. The balance of survival (anti-apoptotic) and death (apoptotic) signals through NF-kB activation cascades results in normal bone homeostasis and healthy remodeling [33]. NF-kB activation through RANK/RANKL signal pathway is closely related with osteoporosis and the bone resorbing activity of osteoclasts [34]. After RANKL binds to RANK, a key step in downstream signaling is binding of TRAFs to specific sites of the cytoplasmic domain of RANK [20,35,36]. Stimulation of RANK results in strong NF-kB activation. NF-kB1(p50) is an important signal for osteoclast development and osteoclast function. Moreover, NF-kB can be activated in many signal cascades related to bone metabolism, e.g., NF-kB activation through the Fas/FasL system leads to enhanced osteoclastogenesis and reduced apoptosis, and it may also cause apoptotic cell death in osteoblasts [37].
Population stratification is an important source of spurious association in genetic association studies [38,39]. However, these factors are unlikely to exist in our sample to interfere with our pathway-based association results. Our cohort came from an apparently homogenous Chinese north-west Han ethnicity population, living in Xi'an city and its neighboring areas. More importantly, in the analyses using Structure 2.2 [40], all subjects used in our study consistently clustered together as a single group, suggesting no significant population substructure. In the analysis using EIGENSTRAT [17], the measure for population stratification, l, was very close to 1, which also suggests no stratification in our cohort. For the above reasons, our association results are unlikely to be plagued by spurious associations due to population admixture/stratification.
In summary, we applied a pathway-based analysis method to identify functional gene sets associated with hip OF. The significant enrichment of the TRAIL pathway genes among the top ranking genes associated with hip OF, together with the pathway's functional relevance to bone metabolism, strongly supports the important role of TRAIL in human hip OF risk. Further detailed and specific functional studies of the TRAIL pathway will be pursued to provide new insights into the etiology of hip OF.