Hip fracture is the most devastating osteoporotic fracture type with significant morbidity and mortality. Several studies in humans and animal models identified chromosomal regions linked to hip size and bone mass. Previously, we identified that the region of 4q21-q41 on rat chromosome (Chr) 4 harbors multiple femoral neck quantitative trait loci (QTLs) in inbred Fischer 344 (F344) and Lewis (LEW) rats. The purpose of this study is to identify the candidate genes for femoral neck structure and density by correlating gene expression in the proximal femur with the femoral neck phenotypes linked to the QTLs on Chr 4. RNA was extracted from proximal femora of 4-wk-old rats from F344 and LEW strains, and two other strains, Copenhagen 2331 and Dark Agouti, were used as a negative control. Microarray analysis was performed using Affymetrix Rat Genome 230 2.0 arrays. A total of 99 genes in the 4q21-q41 region were differentially expressed (P < 0.05) among all strains of rats with a false discovery rate <10%. These 99 genes were then ranked based on the strength of correlation between femoral neck phenotypes measured in F2 animals, homozygous for a particular strain's allele at the Chr 4 QTL and the expression level of the gene in that strain. A total of 18 candidate genes were strongly correlated (r2 > 0.50) with femoral neck width and prioritized for further analysis. Quantitative PCR analysis confirmed 14 of 18 of the candidate genes. Ingenuity pathway analysis revealed several direct or indirect relationships among the candidate genes related to angiogenesis (VEGF), bone growth (FGF2), bone formation (IGF2 and IGF2BP3), and resorption (TNF). This study provides a shortened list of genetic determinants of skeletal traits at the hip and may lead to novel approaches for prevention and treatment of hip fracture.
- quantitative trait loci
- gene expression
- bone density
- hip fracture
hip fracture is one of the most devastating osteoporotic fracture types causing significant morbidity and mortality (34). The primary skeletal determinants of hip fracture risk are bone mineral density (BMD), structure, and strength (14, 16, 31). Several studies have demonstrated that there is a substantial genetic component to hip fragility (6, 17). Our goal is to better understand the genetic basis for variation of skeletal traits at the femoral neck as these provide useful information about hip fragility.
Genetic linkage studies in humans have demonstrated that multiple chromosomal regions are linked to hip BMD, axis length and bone size (15, 20, 24, 30, 33). The same approach in rodent models has also detected quantitative trait loci (QTLs) syntenic with skeletal QTLs in the humans (2, 3). Unfortunately, none of these studies has yet identified a gene or genes underlying the variation in these hip phenotypes. The major obstacle to identifying gene/s in a QTL is the expensive and time-consuming process of narrowing a QTL to a few candidate genes that can be tested further. Moreover, multiple genes often contribute to QTLs having large phenotypic effect sizes (9). In some cases, these genes may in turn also contribute to the variation in other correlated traits.
Previously, we showed that skeletal mass, structure, and strength varies among inbred strains of rats in a site-specific manner (41). We also demonstrated that despite similar body size, substantial variation exists in bone geometry and biomechanical properties between adult Fischer 344 (F344) and Lewis (LEW) rats. Subsequently, we generated a large number of second filial (F2) progeny derived from this cross of F344 and LEW inbred rats and detected QTLs on chromosome (Chr) 4 exclusively influencing femoral neck density and structure (3). Evidence of significant linkage (P < 0.01) of femoral neck width, trabecular area, total volumetric BMD (vBMD), total area, and polar moment of inertia (Ip) were observed in the region of 4q21-q41 with logarithm of the odds (LOD) scores of 19.6, 16.4, 11, 8.9 and 8.0, respectively.
The evidence for linkage of multiple femoral neck phenotypes on Chr 4 in F344 and LEW rats suggests either that there are genes having pleiotropic effects on bone structure or there is a cluster of genes contributing individually to these phenotypes. We also completed a genetic linkage study of femoral neck traits in Copenhagen 2331 (COP) and Dark Agouti (DA) rats in which we found no significant QTLs on Chr 4 (2). We infer from these data that the F344/LEW incorporated genetic polymorphisms that contributed to femoral neck phenotypes whereas the COP/DA cross did not (Fig. 1). At present, F344/LEW and COP/DA rat strains are the only available rat models in which bone QTLs have been identified.
The purpose of this study is to identify the candidate genes for femoral neck density and structure by using our previous genetic mapping data for F344/LEW and COP/DA rats. Using Affymetrix microarrays, we identified differentially expressed genes in the F344 and LEW rats, and then compared with gene expression in COP and DA rats. We prioritized for further analysis those genes whose expression was strongly correlated with femoral neck phenotypes measured in F2 animals from the rat crosses.
A structured network-based approach to analyze genome-wide gene expression in the context of known functional interrelationship among genes, proteins, and small molecules has been successfully used for various systems including bone (12, 25). We applied this knowledge-based system for analysis of the related pathways among prioritized candidate genes based on their molecular function, biological process, and cellular component of the gene products.
MATERIALS AND METHODS
In previous studies, we obtained femoral neck phenotypes from 595 female F2 offspring derived from F344 and LEW progenitors and 423 female F2 offspring derived from COP and DA progenitors (2, 3). At 26 wk of age the rats were euthanized, and left femurs were collected for densitometry analysis. vBMD (mg/cm3), cross-sectional area (mm2), and Ip (mm4) of the femoral neck were measured from peripheral quantitative computed tomography slices using the XCT Research SA Plus, software version 5.40 as described previously (3). The femoral neck width of each rat was measured with digital calipers. Genomic DNA was isolated from the spleen of individual F2 rats, and genotyping for each animal was accomplished using PCR with microsatellite markers as described previously (4). Quantitative linkage analysis for femoral neck total vBMD, total area, trabecular area, Ip, and neck width was performed as described previously (4).
RNA extraction and microarray analysis.
Whole femora were harvested from 4-wk-old F344, LEW, COP, and DA animals and were immediately frozen in liquid nitrogen and stored at −80°C until required. RNA from the proximal femurs was extracted (n = 4 per strain) using TRIzol (Invitrogen, Carlsbad, CA) and RNeasy Mini Kit (Qiagen, Valencia, CA). Briefly, frozen femurs were placed into a mortar containing liquid nitrogen and crushed with a pestle. The crushed bones were immediately transferred into TRIzol and homogenized using a Polytron homogenizer (PRO Scientific, Monroe, CT). RNA was then isolated according to manufacturer's instructions and was treated with a DNA-free kit (Ambion, Austin, TX) to remove any residual DNA. RNA quality was determined using a 2100 Bioanalyzer (Agilent, Palo Alto, CA) and was quantified using a spectrophotometer (NanoDrop, Wilmington, DE). For microarray analysis, 5 μg of total RNA from each sample was reverse transcribed and linear amplified using a primer containing oligo (dT) and a T7 RNA polymerase promoter (Invitrogen). After synthesis of the first and second strands of cDNA according to the standard protocols from Affymetrix (Affymetrix: GeneChip® Expression Analysis Technical Manual, Santa Clara, CA), the product was used in an in vitro transcription reaction to generate cRNA in the presence of biotin-labeled dNTPs using an RNA Transcript Labeling Kit (Affymetrix P/N 900449). The labeled cRNA was purified using Affymetrix GeneChip sample cleaning module (Affymetrix P/N 90037) to remove free nucleotides and fragmented nonezymatically. We hybridized 10 mg of fragmented biotin-labeled cRNA of each sample to Rat Genome 230 2.0 Array oligo microarrays (Affymetrix P/N 511056) for 17 h at 45°C with constant rotation. The GeneChip was then washed and stained in the Affymetrix Fluidics Station 400 following standard procedures. Subsequently, each array was scanned by the Agilent GeneChip Scanner 3000 (Agilent Technologies). The Rat 230 GeneChip Microarray has 31,000 probe sets representing 28,700 well-substantiated rat genes.
Quality control for RNA and Affymetrix data.
Quality control (QC) of RNA was done by measuring the ratio between signals from 5′ and 3′ end of the GAPDH and β-actin genes (3′/5′ ratios) and by implementing the RNA degradation plots. Affymetrix data QC was done by determining the percentage of present calls (detection calls) and the scale factors between the arrays. Principal component analysis was conducted to ensure that the overall gene expression profiles in all the samples in each experimental condition are correctly correlated.
Microarray data analysis and informatics.
The images from each array were analyzed using Affymetrix GeneChip Operating system 1.2 software. The .cel files were imported into the statistical programming environment R (38) for further analysis with tools available from the Bioconductor Project (18). Expression data were normalized and log2 transformed using the robust multichip average (RMA) method (11, 21) implemented in the Bioconductor package RMA. Mapping of probe sets to chromosomal location was accomplished with data provided by Affymetrix. The identities of the probe sets were confirmed by comparing the target mRNA sequences on the Affymetrix Rat Genome 230 2.0 GeneChip with the National Center for Biotechnology Information (NCBI) GenBank database (http://www.ncbi.nlm.nih.gov/Genbank/). To increase power and decrease the false discovery rate (FDR), we only analyzed probe sets that were reliably detected on all of the microarrays, based upon the detection call generated by the Affymetrix Microarray Analysis Suite 5.0 algorithm (27).
Quantitative real-time PCR analyses.
Five micrograms of total RNA (same RNA used for Affymetrix analysis) was reverse transcribed using Superscript III reverse transcription reagent for first-strand cDNA synthesis (Invitrogen). The top 18 strongly correlated (r2 > 0.50) candidate genes were selected for verification of microarray data by real-time PCR analysis. All real-time PCR reactions contained the first-strand cDNA corresponding to 25 ng of total RNA. TaqMan predeveloped primers, FAM dye-labeled MGB probes and Universal Master Mix (Applied Biosystems) were used to quantify the relative gene expression. Rat β-actin was used as an endogenous control. Real-time detection of PCR products was performed using ABI PRISM 7300 sequence detector (Applied Biosystems). Relative expression of mRNA was calculated based on a relative standard curve and normalized to β-actin. All quantitative real-time PCR (qPCR) analysis used triplicate independent samples.
The gene-wise P values for microarray analyses between two strains were calculated by Welch's t-test using the package Limma (37), and FDR was calculated by the method of Benjamini and Hochberg (10). Probe sets were considered differentially expressed if the FDR was <10%. The gene-wise P values among all strains were calculated by ANOVA in the same way with the same FDR setting. The microarray data set and MIAME-compliant experimental conditions have been exported to the NCBI Gene Expression Omnibus (GEO) Express web portal (GEO accession number GSE 11180; URL: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE11180). Regression analysis was performed with the average gene expression level for the strain as the dependent variable and the phenotypic mean value in animals homozygous for that strain's allele at the marker underlying the QTL as the independent variable. The proportion of variation (r2 value) in the phenotypic means explained by the variation in gene expression was obtained. This process was repeated for each of the five femoral neck phenotypes of interest. All statistical analyses were performed using the statistical software package StatView (Abacus Concepts, Berkeley, CA).
The interactions between differentially expressed genes for femoral neck phenotypes and all other genes were investigated using Ingenuity Pathway Analysis (IPA 5.0; Ingenuity Systems, Mountain View, CA). Those differentially expressed genes which also explained a significant proportion of the phenotypic variation in the femoral neck phenotypes were uploaded into the application. Each gene identifier was mapped to its corresponding gene in the Ingenuity Pathway Knowledge Base (IPKB). These genes were overlaid onto a global network developed from the information contained in the IPKB. Networks of these genes, defined as the reflection of all interactions of a given gene defined in the literature, were then algorithmically generated based on their connectivity. The interactions indicate physical association, induction/activation, or repression/inactivation of one gene product by the other, directly or through another intermediary molecule.
Femoral neck phenotypes as a function of genotypes.
Genotypic means for femoral neck phenotypes at D4Rat231 (43.1 cM) in the F344 × LEW and at D4Rat103 (41.93 cM) in the COP × DA F2 rats on Chr 4 are summarized in Table 1. As expected due to the QTL observed at marker D4Rat231, all of the mean values for the femoral neck phenotypes were significantly different (P < 0.0001) between the F2 animals homozygous for the F344 (f/f) and the LEW (l/l) F2 alleles. In contrast, for four of the five femoral neck phenotypes at marker D4Rat103, the mean values for the F2 animals homozygous for the COP (c/c) and DA (d/d) F2 alleles were not significantly different [with the exception of femoral neck trabecular area that was significantly higher (P = 0.01) in d/d group compared with c/c group].
A total of 99 out of 458 genes were differentially expressed among the rat strains with a FDR < 10%. Of these 99, 14 were differentially expressed (P < 0.05) between COP and DA rats. Since there was no evidence of femoral neck density or structure QTLs on Chr 4 in the COP × DA F2 rats, these 14 genes were excluded from further analysis and the remaining 85 genes were analyzed by regression. All Affymetrix gene information was collected with the data available through NCBI's GenBank database (http://www.ncbi.nlm.nih.gov/Genbank/).
Genes explaining skeletal phenotypes.
Regression analyses were performed with the 85 prioritized genes to test whether the gene expression in a particular inbred strain explained a substantial proportion of the variation (>50%) in femoral neck width of F2 offspring from both the F344 × LEW and COP × DA crosses. A total of 30 genes, including 18 candidate genes (Table 2) and 12 predicted genes were found to explain >50% of the variation in femoral neck width. The r2 values in Table 2 are sorted based on the femur neck width, which had the highest LOD score QTL obtained from our previous study. Among the top 18 candidate genes determined by regression analysis (Table 2), 14 genes were confirmed to have similar correlations with the femoral neck density and structure phenotypes by qPCR (Table 3). Surprisingly, two of the genes (Pdia4 and Vps24) that showed high correlations (>95%) by Affymetrix analysis were not confirmed by qPCR. Of the 14 genes confirmed by qPCR, only six (Aqp1, Gpnmb, Igf2bp3, Ptprz1, Serbp1, and St3gal5) were highly correlated (>50%) with all femoral neck parameters, whereas two others (Chchd4 and Spr) achieved high correlations with most of the femoral neck parameters.
The eight candidate genes that were highly correlated with femoral neck phenotypes by qPCR were mapped to pathways using Ingenuity. Genes were directly or indirectly connected to insulin-like growth factor 2 (IGF2), fibroblast growth factor 2 (FGF2), vascular endothelial growth factor (VEGF), extracellular signal-related kinase (ERK), nitric oxide synthase (NOS), and tumor necrosis factor pathways (TNF) (Fig. 2).
In this study, we identified candidate genes in the region between 4q21-q41 on rat Chr 4 that has been linked with several distinct femoral neck density and structure phenotypes, i.e., neck width, total vBMD, total area, trabecular area, and Ip. These genes were differentially expressed between F344 and LEW inbred rat strains, and their expression levels explained a significant proportion of the variation in femoral neck density and structure phenotypes in F2 populations. Using qPCR, we confirmed 14 of the 18 top candidate genes. Among the candidate genes, several were previously demonstrated to have a role in bone metabolism, including aquaporin 1 (Aqp1), glycoprotein (transmembrane) nmb (Gpnmb), serpine1 mRNA binding protein 1 (Serbp1), and IGF2-binding protein 3 (Igf2bp3). Pathway analysis identified several cytokines and growth factors known to affect bone metabolism, including IGF2, FGF2, ERK, VEGF, and TNF.
A strong candidate gene for bone traits is Igf2bp3. IGF2, previously known as skeletal growth factor (29), stimulates osteoblasts to form new bone. Binding proteins for IGFs also have been shown to affect bone biology, for instance IGFBP2 has been associated with increased bone turnover and decreased bone mass in both men and women (5). However, IGF2BP3 (or Igf2bp3) is different from previously studied IGF binding proteins because it binds to mRNA, rather than IGF protein, and thus exerts transcriptional regulation of IGF2. IGF2BP3 has not yet been studied with regard to bone biology. Other candidate genes include Gpnmb (involved in osteoclast differentiation) (35), Aqp1 (aquaporins have been shown to regulate osteoclast fusion) (1), and Serpbp1 (serine proteinases are involved in bone resorption) (40).
The pathway analysis revealed several direct or indirect relationships among the candidate genes we obtained from this study with the genes already reported to be related to bone metabolism. Among them IGF2 (discussed above), FGF2 that regulates bone growth (13, 26), and VEGF that promotes angiogenesis and osteoblastogenesis (42). The members of the TNF family of ligands and receptors are critical regulators of osteoclastogenesis (19, 43). Other interesting pathways included pleiotrophin (PTN) and NOS. PTN or heparin binding growth factor 8 has role in chondrogenesis (28, 32). It has been shown that NOS enhances BMD and bone turnover in mice (36), mediates osteoclast activity in vivo and vitro (22), and regulates terminal differentiation of chondrocytes (39). Also, deletion of inducible NOS gene impairs fracture healing in mice (7). Also, the role of ERK/MAP kinase (MAPK) pathway has been shown in all bone cells.
We undertook a novel approach to prioritize candidate genes by analyzing the microarray-based expression of genes underlying the QTLs at rat 4q21-q41 and correlating them with femoral neck phenotypic values. This QTL region on Chr 4 is homologous to regions of mouse chromosome 6 and human chromosomes 1p, 2p, 3p, 3q, 7p, and 7q (Table 3). Linkage to mouse Chr 6 was previously reported for femur BMD and strength (8, 24). Also, the homologous region in human 7p14 was previously linked to femoral neck BMD (33). We used young female rats for gene expression analysis since we wanted to assess gene regulation of bone accrual. Thus, it is likely that we identified genes regulating bone growth that contribute to the accumulation of peak bone mass. The genes and regulatory networks identified in this study will provide insights into the genetic determinants of skeletal trait at hip in humans as QTL genes identified in rat model have shown remarkable relevance to related human disease phenotypes. Of course, our approach has several limitations. Gene expression study based on microarray does not capture effects of alternative gene splicing or posttranslational modification of proteins. Also, we could not account for the genetic contributions from other QTLs, and these may confound our results. This problem would be solved by using congenic or consomic rats in which the Chr 4 QTL is isolated from the effects of other QTLs throughout the genome.
This work was supported by National Institutes of Health Grants: AR-047822, AG-018397, and AA-10707. Microarray analyses were carried out in the Center for Medical Genomics at Indiana University School of Medicine, which is supported in part by the Indiana Genomics Initiative (INGEN, supported in part by The Lilly Endowment Inc.).
Address for reprint requests and other correspondence: C. H. Turner, Dept. of Biomedical Engineering, IUPUI, 1120 South Dr., Fesler Hall 115, Indianapolis, IN 46202 (e-mail:).
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
- Copyright © 2008 the American Physiological Society