- Email: [email protected]
Differentially expressed genes strongly correlated with femur strength in rats
Genomics 94 (2009) 257–262
Contents lists available at ScienceDirect
Genomics j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y g e n o
Differentially expressed genes strongly correlated with femur strength in rats Imranul Alam a, Qiwei Sun a, Daniel L. Koller b, Lixiang Liu b, Yunlong Liu c, Howard J. Edenberg d, Jiliang Li e, Tatiana Foroud b, Charles H. Turner a,⁎ a
Department of Biomedical Engineering, Indiana University Purdue University Indianapolis (IUPUI), 1120 South Drive, Fesler Hall 115, Indianapolis, IN 46202-5251, USA Department of Medical and Molecular Genetics, Indiana University Purdue University Indianapolis (IUPUI), IN, USA Department of Medicine, Indiana University Purdue University Indianapolis (IUPUI), IN, USA d Department of Biochemistry and Molecular Biology, Indiana University Purdue University Indianapolis (IUPUI), IN, USA e Department of Biology, Indiana University Purdue University Indianapolis (IUPUI), IN, USA b c
a r t i c l e
i n f o
Article history: Received 30 January 2009 Accepted 25 May 2009 Available online 29 May 2009 Keywords: Femur strength Gene expression Microarray QTLs Osteoporotic fracture
a b s t r a c t The region of chromosome 1q33–q54 harbors quantitative trait loci (QTL) for femur strength in COP × DA and F344 × LEW F2 rats. The purpose of this study is to identify the genes within this QTL region that contribute to the variation in femur strength. Microarray analysis was performed using RNA extracted from femurs of COP, DA, F344 and LEW rats. Genes differentially expressed in the 1q33–q54 region among these rat strains were then ranked based on the strength of correlation with femur strength in F2 animals derived from these rats. A total of 214 genes in this QTL region were differentially expressed among all rat strains, and 81 genes were found to be strongly correlated (r2 N 0.50) with femur strength. Of these, 12 candidate genes were prioritized for further validation, and 8 of these genes (Ifit3, Ppp2r5b, Irf7, Mpeg1, Bloc1s2, Pycard, Sec23ip, and Hps6) were confirmed by quantitative PCR (qPCR). Ingenuity Pathway Analysis suggested that these genes were involved in interferon alpha, nuclear factor-kappa B (NFkB), extracellular signal-related kinase (ERK), hepatocyte nuclear factor 4 alpha (HNF4A) and tumor necrosis factor (TNF) pathways. © 2009 Published by Elsevier Inc.
Introduction Osteoporosis is a common polygenic disorder characterized by reduced bone mineral density (BMD) and compromised bone quality with microarchitectural deterioration leading to an increased susceptibility to fracture at multiple skeletal sites [1]. BMD, structure and strength are the primary skeletal determinants of osteoporotic fracture risk [2–4] and are under substantial genetic control [5–8]. Most of the studies identifying osteoporosis genes primarily focus on BMD as it is the most commonly available tool for predicting fracture risk. However, BMD is not an entirely adequate measure of bone strength and structure. Bone strength, the phenotype that best reflects skeletal fragility also depends on bone geometry (size and shape), microarchitecture and material properties. Bone strength is under the control of multiple genes, each of which has a modest contribution to the genetic variance [15,20]. Identification of genes underlying bone strength will reveal valuable insight regarding the genetics of osteoporosis and fracture risk. Over many years, several linkage studies in human and experimental animal models have identified quantitative trait loci (QTLs) linked to bone structure and strength [9–21], but the genes underlying these phenotypes are largely unknown. Identification of candidate
⁎ Corresponding author. Fax: +1 317 278 9568. E-mail address: [email protected] (C.H. Turner). 0888-7543/$ – see front matter © 2009 Published by Elsevier Inc. doi:10.1016/j.ygeno.2009.05.008
genes within the broad QTL intervals, which typically harbor hundreds of genes, is extremely challenging and requires multiple complementary approaches. One strategy is the development of a congenic strain, which can narrow the critical QTL to a small chromosomal segment harboring fewer candidate genes. However, due to the extensive breeding that is required, this approach is quite time-consuming and expensive. Other approaches include creation of recombinant inbred lines and haplotype-based analyses. However, these strategies do not produce mapping with high enough resolution to detect the causal genes. An alternative strategy involves an integrative genetic approach combining linkage with genomic expression analysis in multiple inbred lines already used for QTL analysis. Such studies are extremely difficult in humans, but they can be undertaken in animal models. Previously, we showed that skeletal mass, structure and strength vary among inbred strains of rats in a site-specific manner [22]. We demonstrated that despite similar body size, substantial variation exists in bone geometry and biomechanical properties between adult Copenhagen 2331 (COP), Dark Agouti (DA), Fischer 344 (F344) and Lewis (LEW) rats [22]. Subsequently, we generated a large number of second filial (F2) female progeny derived from (COP × DA) and (F344 × LEW) crosses, and detected significant linkage for femur ultimate force (Fu) (a biomechanical measurement that determines the maximum load that the bone can support before failing) within the region of 1q33–q54 in both of these crosses [19,20]. Thus, 1q33–q54 region harbors QTLs for a bone phenotype directly related to femur strength.
258
I. Alam et al. / Genomics 94 (2009) 257–262
The purpose of this study was to narrow the list of candidate genes within the region of 1q33–q54 contributing to the variation in femur strength among COP, DA, F344 and LEW rats. We chose these strains of rats because, at present, these are the only available rat models in which bone QTLs have been identified [18–21]. We used the data from microarray-based gene expression analyses to identify genes differentially expressed among these rats. These genes were then ranked based on their correlation with Fu in F2 animals derived from these rat strains. The expression levels of several strongly correlated (r2 N 0.9) candidate genes were confirmed using quantitative real-time PCR (qPCR). In addition, we applied Ingenuity Pathway Analysis (IPA) for analysis of the related pathways among prioritized candidate genes based on their molecular function, biological process and cellular component of the gene products. IPA is a structured, network-based system to analyze genome-wide gene expression in the context of known functional interrelationships among genes, proteins, and small molecules. This approach has been successfully used for various systems including bone [35,36].
Microarray analyses and correlation with femur phenotypes The Affymetrix microarray analyses showed that a total of 214 genes residing in the region of 1q33–q54 (among the 999 total genes on the microarray probe sets) were differentially expressed among COP, DA, F344 and LEW rats with a FDR less than 10%. Regression analyses were performed for all 214 differentially expressed genes to assess the correlation between gene expression and the femur Fu among F2 animals homozygous for that strain's allele at the markers (D1Rat69 for COP × DA; D1Rat291 for F344 × LEW) underlying the QTL. The 1q33–q54 QTL for Fu was replicated in both (COP × DA) and (F344 × LEW) crosses; therefore, we identified candidate genes that showed strong correlation (r2 N 0.50) for this phenotype in these rats. We considered any r-value of above 0.7, which corresponds to a r2-value of above 0.49, as a strong correlation. A total of 81 genes, including 72 candidate genes and 9 predicted genes (Supplementary Table 1), were detected using this method. Among these 72 genes, 12 were ranked as the highest priority because they were found to explain more than 90% (r2 N 0.90) of the variation of femur Fu (Supplementary Table 1).
Results qPCR analyses Genetic loci for femur Fu on 1q33–q54 Evidence of significant linkage of femur Fu was detected in the region of 1q33–q54 in COP × DA F2 animals with LOD scores of 6.7 (Fig. 1) [19]. In addition, significant evidence of linkage of femur Fu was observed in F344 × LEW F2 rats with a LOD score of 4.8, in a broad region that overlapped the COP × DA peak (Fig. 1) [20]. Femur Fu as a function of genotypes Genotypic means for femur Fu was measured on Chr 1 at D1Rat69 (101.5 cM) in the COP × DA F2 rats and at D1Rat291 (101.4 cM) in the F344 × LEW F2 rats. Since a QTL was detected at marker D1Rat69 in the COP × DA cross and at marker D1Rat291 in the F344 × LEW cross, the mean values for the femur Fu were significantly different between the F2 animals homozygous for the COP (c/c) and the DA (d/d) F2 alleles (120.35 ± 1.81 vs. 129.08 ± 1.87; p = 0.001) and the F2 animals homozygous for the F344 (f/f) and LEW (l/l) F2 alleles (130.47 ± 1.11 vs. 138.02 ± 1.08; p b 0.0001).
qPCR analyses of the 12 prioritized candidate genes in Supplementary Table 1 confirmed the strong correlation between gene expression and femur ultimate force (Table 1). The correlation coefficients for Pold4, Ifit3, Hps6, Mrpl3, Map4k2, Sec23ip, Rras2, Pycard, Ppp2r5b, Irf7, Mpeg1 and Bloc1s2 between microarray and qPCR analyses were 0.65, 0.99, 0.70, 0.70, 0.24, 0.81, 0.22, 0.94, 0.96, 0.99, 0.94 and 0.89, respectively, indicating good concordance between the two methods for all genes except Map4k2 and Rras2. Among these 12 candidate genes, 8 genes (Ifit3, Ppp2r5b, Irf7, Mpeg1, Bloc1s2, Pycard, Sec23ip and Hps6) were confirmed to have similar correlations with femur strength by qPCR (r2-values ranging from 0.59 to 0.97) (Table 1). Pathway analysis The 8 candidate genes that were confirmed by qPCR were mapped to pathways using Ingenuity Pathway Analysis. These genes were directly or indirectly connected to interferon alpha, nuclear factorkappa B (NFkB), extracellular signal-related kinase (ERK), hepatocyte nuclear factor 4 alpha (HNF4A) and tumor necrosis factor (TNF) pathways (Fig. 2). Discussion
Fig. 1. Linkage analysis for femur ultimate force (Fu) for rat chromosome 1 in COP × DA and F344 × LEW crosses in the region of 1q33–q54 (D1Rat52–D1Rat82). LOD scores plotted on the y axis vs. the relative microsatellite markers located along chromosome 1 on the x axis.
Genetic mapping for complex traits usually identifies chromosomal regions that are typically quite broad in nature. The 1q33–q54 QTL region linked to various femur phenotypes also encompasses a broad region of roughly 60 cM of the rat genome, harboring nearly 1000 genes. In addition, many of these genes might have pleiotropic effects on different femur phenotypes. Also, there might be a cluster of genes contributing individually to these phenotypes. To prioritize candidate genes for further analysis, we employed an integrative genetic approach combining linkage, gene expression and physical traits. We took the advantage of microarray-based expression analyses to simultaneously screen thousands of genes for their differential expression in COP, DA, F344 and LEW rats, as these are the only available rat models in which bone QTLs have been identified [18–21]. In addition, we used 4-week-old rats for gene expression analyses since we wanted to assess gene regulation of bone accrual and to identify genes contribute to the accumulation of peak bone mass. By comparing the genes differentially expressed among these rats, we were able to significantly reduce the number of candidate genes in the 1q33–q54 QTL region to 214.
I. Alam et al. / Genomics 94 (2009) 257–262 Table 1 r-square values from qPCR analysis of top 12 candidate genes for femur ultimate force within 1q33–q54 region on Chr 1 in COP, DA, F344 and LEW ratsa. Gene symbol
Gene name
Rat Femur ultimate location force (r2)
Human synteny
Ifit3
Interferon-induced protein with tetratricopeptide repeats 3 Protein phosphatase 2, regulatory subunit B, delta isoform Interferon regulatory factor 7 Macrophage expressed gene 1 Biogenesis of lysosome-related organelles complex-1, subunit 2 PYD and CARD domain containing SEC23 interacting protein Hermansky–Pudlak syndrome 6 Polymerase (DNA-directed), delta 4 Mitochondrial ribosomal protein L43 Related RAS viral (r-ras) oncogene homolog 2 Mitogen activated protein kinase kinase kinase kinase 2
0.97
1q52
10q23.3
0.85
1q43
11q13
0.88 0.72 0.74
1q41 1q43 1q54
11p15.5 11q12.1 10q24.31
0.65 0.61 0.59 0.45 0.34 0.24
1q36 1q37 1q54 1q42 1q54 1q34
16p12 10q26.11 10q24.3 11q13.1 10q24.31 11p15.2
0.20
1q43
11q13
Ppp2r5b Irf7 Mpeg1 Bloc1s2 Pycard Sec23ip Hps6 Pold4 Mrpl43 Rras2 Map4k2 a
Confirmed genes (r2 N 0.50) are indicated in bold face.
Prioritization of microarray-based candidate genes is usually done by ranking the genes for the magnitude of expression differences (fold-differences) between different strains. However, with complex traits like osteoporosis, even subtle changes in gene expression could be important. Thus, stronger evidence for the association between a gene and a trait might be identified by examining the correlation between gene expression and a given physical trait. Using this
259
approach, among the 214 candidate genes, we were able to prioritize 72 genes that were strongly correlated (r2 N 0.5) with femur strength. Furthermore, by investigating genes those were correlated with more than 90% variation (r2 N 0.9) with femur ultimate force, we further narrowed the list of genes to 12 candidate (Supplementary Table 1). Among these 12 candidate genes, 8 genes (Ifit3, Ppp2r5b, Irf7, Mpeg1, Bloc1s2, Pycard, Sec23ip and Hps6) were confirmed to have similar strong correlations (r2-values ranging from 0.59 to 0.97) with the femur strength when expression was measured by qPCR (Table 1). Surprisingly, four genes (Pold4, Mrpl43, Map4k2 and Rras2) that showed high correlations for femur strength by Affymetrix analysis were not confirmed by qPCR; the difference may relate to probe location for the arrays at the 3′ end of the target gene. The 8 candidate genes that were confirmed by qPCR were associated with pathways directly or indirectly involved in controlling cell morphology and development, cellular assembly and organization, cell-to-cell signaling, molecular transport, protein trafficking, and cellular function and maintenance (Fig. 2). Genes in the canonical pathways (networks generated using IPA) were related to interferon signaling and activation of pattern recognition receptors. When these genes were categorized based on location of cellular components Ifit3, Mpeg1, Bloc1s2, Pycard, Sec23ip and Hps6 are located in cytoplasm while Ppp2r5b and Irf7 are located in the nucleus. Recently, signaling networks with interferon and other cytokines have been shown to be important in bone metabolism [30]. Interestingly, the top 2 candidate genes (Ifit3 and Irf7) identified in this study are interferon-related genes and pathway analysis revealed several interferon molecules involved in the same network (Fig. 2). Interferon regulatory factors Irf1 and Irf2 have been associated with osteoblast apoptosis [37] but roles for Ifit3 and Irf7 in bone biology have not yet been studied. However, both
Fig. 2. Network of the 8 candidate genes that were confirmed by qPCR (highlighted in red) for femur strength in Ingenuity Pathway Analysis (IPA). The well-known NFkB, ERK and TNF pathways related to bone metabolism is highlighted in blue.
260
I. Alam et al. / Genomics 94 (2009) 257–262
interferon alpha and beta signaling has been shown to be important for maintaining normal bone mass by regulating osteoclastic bone resorption [38]. Interferon receptor alpha 1 knockout mice exhibit an osteoporosis phenotype due to increased osteoclastogenesis and decreased trabecular bone mass. Interferon-beta inhibits RANK-ligand mediated osteoclastogenesis by inhibiting the expression of c-FOS and mice lacking interferon-beta have severe osteopenia due to enhanced osteoclastogenesis [38]. In addition, Ppp2r5b (Ppp2r2d) has been previously shown to have role in TGF-beta signaling [27] and Pycard is a regulator of NFkB activation pathway [28]. Hps6 is associated with Hermansky–Pudlak syndrome, an autosomal recessive syndrome due to a mutation in clathrin-associated adaptor protein [32]. Because clathrin is important in the WNT signaling pathway in bone [33], Hps6 might be involved in the molecular regulation of the WNT pathway. Bloc1s2 has been shown to regulate apoptosis in cancer [31], and Sec23ip is involved in protein export from the endoplasmic reticulum [34], and both of these genes are indirectly related to Hps5/6 genes. In a recent study, Sec23ip has been shown to be essential for normal craniofacial development [39]. The candidate genes and the interconnected molecules identified in this study are thus related to pathways involving both bone formation and resorption: nuclear factor-kappa B (NFkB) regulates osteoclast mediated bone resorption [40,41]; extracellular signal-related kinase (ERK) has critical role in osteoblast differentiation and bone development [42] as well as osteoclast differentiation and survival [43]; members of tumor necrosis factor (TNF) family of ligands and receptors are critical regulators of osteoclastogenesis [44,45]. However, further functional analysis using gene targeting strategies such as knockout and siRNA will be necessary to fully characterize the roles of these candidate genes in bone metabolism. The chromosomal region of the candidate genes for femur strength on Chr 1 in rats is homologous to regions of mouse chromosomes 7, 19 and human chromosomes 10q, 11p, 11q, and 16p (Table 1). Linkage to mouse Chr 7 was previously reported for femur strength and structure [16–17]. Furthermore, the homologous regions in human 10q21, 10q26 and 11p15 were previously linked to wrist bone size, hip BMD and spine bone size, respectively [12,13,29]. In this study, using an integrative genetic approach, we identified several top candidate genes underlying the QTL region of 1q33–q54 that were differentially expressed and strongly correlated with femur strength in rats. The approach we took here is complementary to conventional strategies for the dissection of complex traits and will facilitate the QTL to gene discovery process. However our approach has several limitations. Gene expression study based on microarray analysis explains only transcriptional regulation of genes and does not capture effects of alternative gene splicing or post-translational modification of proteins. In addition, the whole femur RNA came from a variety of different types of bone cells including bone marrow and thus we cannot be sure that the gene expression changes originated only from bone cells. Further studies involving different cell lines are thus necessary to identify which specific cell type expresses the genes identified in this study. Also, our results could be influenced by the genetic contributions from other QTLs from the same or other chromosomes. Some of these problems could be solved by isolating the Chr 1 QTL region from the effects of other QTLs throughout the genome by using congenic rats. Also, whether the same candidate genes regulate femur strength in male rats remains to be determined, as we only studied female rats in this study. Materials and methods In previous studies, we created 423 and 595 female 2nd filial (F2) offspring from COP × DA and F344 × LEW crosses, respectively. The rats were allowed to grow to 6 months (26 weeks) of age before sacrifice. The lower limbs and spleen were dissected out. The right side of the
lower limbs were immediately stored at −20 °C for biomechanical analysis. The spleens were immediately frozen in liquid nitrogen before transferring them to −80 °C. The right femurs were tested in three-point bending by positioning them on the lower supports of a three-point bending fixture and applying load at the midpoint using a material testing machine (Alliance RT/5, MTS Systems Corp., Eden Prairie, MN, USA) and ultimate force (Fu; N) was obtained from the force vs. displacement curves in TestWorks software, version 4.06 as described previously [19,20]. Genotyping for each animal was accomplished with microsatellite markers using genomic DNA isolated from spleen previously shown to be polymorphic for COP, DA, F344 and LEW rats with an average interval of 20 cM from centromere to telomere. Chromosomal positions, marker order, and map positions were obtained from the Rat Genome Database (RGD) website (http://rgd.mcw. edu/). Genotypic data in the F2 animals for each marker were tested to ensure the expected Mendelian ratios. Marker maps were generated using MAPMARKER/EXP with data from all available genotyped F2 animals. Multipoint quantitative linkage analysis was then performed to identify chromosomal regions that harbor QTL for femur Fu using the programs R/qtl and QTL/Cartographer. Permutation testing was employed to obtain appropriate genome-wide significance levels for the linkage results and ANOVA was used for the most significant marker in each QTL region to further characterize significant genotypic group differences as described previously [19,20]. RNA extraction and microarray analysis Whole femora were harvested from 4-week-old COP, DA, F344 and LEW animals and were immediately frozen in liquid nitrogen and stored at − 80 °C until required. RNA from femurs was extracted from 4 animals per strain using Trizol (Invitrogen, Carlsbad, CA), followed by further purification using an RNeasy Mini Kit (Qiagen Inc., Valencia, CA) as described previously [23]. RNA quality was determined using a 2100 Bioanalyzer (Agilent, Palo Alto, CA) and was quantified using a spectrophotometer (NanoDrop, Wilmington, DE, USA). For microarray analysis, 5 μg of total RNA from each sample was labeled and hybridized to Affymetrix Rat Genome 230 2.0 GeneChips (Affymetrix, Santa Clara, CA) as described previously [23]. Quality control for RNA and Affymetrix data was done as described previously [23]. Microarray data analysis and informatics The images from each array were analyzed using Affymetrix GeneChip Operating system (GCOS) 1.2 software. Mapping of probe sets to chromosomal location was accomplished with data provided by Affymetrix. The region of 1q33–q54, limited by the markers of D1Rat52 (80.1 cM) and D1Rat82 (139.0 cM), harbors a total of 999 genes. Among these genes, the identities of the 214 differentially expressed genes within this QTL region 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 [25]. FDR was calculated by the method of Benjamini and Hochberg [26]. Probe sets were considered differentially expressed if the FDR was less than 10%. We used 10% FDR because we wanted to reduce the number of false negatives when we identified differentially expressed genes. The microarray data set was submitted to the NCBI Gene Expression Omnibus (GEO) Express web portal (GEO accession number GSE 11180).
I. Alam et al. / Genomics 94 (2009) 257–262
Quantitative real-time PCR (qPCR) analysis Five micrograms of total RNA (the same RNA used for Affymetrix analysis) was reverse transcribed using Superscript III reverse transcription reagent for first-strand cDNA synthesis (Invitrogen, CA). Twelve strongly correlated (r2 N 0.90) candidate genes for femur Fu were selected for verification of microarray data by qPCR analysis. All real-time PCR reactions contained the first-strand cDNA corresponding to 25 ng of total RNA. Real-time detection of PCR products was performed using ABI PRISM 7300 sequence detector (Applied Biosystem, CA) as described previously [23]. Expression of mRNA was calculated based on a relative standard curve and normalized to beta-actin. All qPCR analysis used triplicates of each of 4 biological samples. Statistics Differential expression of each gene among all strains was calculated by ANOVA using the package Limma [24]. Regression analysis was performed with the average gene expression level for the strain as the dependent variable, and the phenotypic mean value of femur Fu 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 using the statistical software package StatView (Abacus Concepts, Inc., Berkeley, CA). Pathway analysis The interactions between differentially expressed genes confirmed by qPCR for femur strength and all other genes were investigated using Ingenuity Pathway Analysis (IPA 5.0; Ingenuity Systems, Inc., Mountain View, CA). The differentially expressed genes 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. Acknowledgments This work was supported by the US National Institutes of Health through the following grants: RO1AR047822, P01AG018397 and AA10707. 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.). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.ygeno.2009.05.008. References [1] J.A. Kanis, L.J. Melton, C. Christiansen, C.C. Johnston, N. Khaltaev, The diagnosis of osteoporosis, J. Bone. Miner. Res. 9 (1994) 1137–1141. [2] S.R. Cummings, D.M. Black, M.C. Nevitt, W. Browner, J. Cauley, K. Ensrud, H.K. Genant, L. Palermo, J. Scott, T.M. Vogt, Bone density at various sites for prediction of hip fractures, The Study of Osteoporotic Fractures Research Group, Lancet 341 (1993) 72–75. [3] K.G. Faulkner, S.R. Cummings, D. Black, L. Palermo, C.C. Gluer, H.K. Genant, Simple measurement of femoral geometry predicts hip fracture: the study of osteoporotic fractures, J. Bone. Miner. Res. 8 (1993) 1211–1217.
261
[4] M. Peacock, C.H. Turner, G. Liu, A.K. Manatunga, L. Timmerman, C.C. Johnston, Better discrimination of hip fracture using bone density, geometry and architecture, Osteoporos. Int. 5 (1995) 167–173. [5] N.K. Arden, J. Baker, C. Hogg, K. Bann, T.D. Spector, The heritability of bone mineral density, ultrasound of the calcaneus and hip axis length: a study of postmenopausal twins, J. Bone Miner. Res. 11 (1996) 530–534. [6] P. Garnero, N.K. Arden, G. Griffiths, P.D. Delmas, T.D. Spector, Genetic influence on bone turnover in postmenopausal twins, J. Clin. Endocrinol. Metab. 81 (1996) 140–146. [7] M. Peacock, C.H. Turner, M.J. Econs, T. Foroud, Genetics of osteoporosis, Endocr. Rev. 23 (2002) 378–383. [8] S.H. Ralston, Genetic determinants of osteoporosis, Curr. Opin. Rheumatol. 17 (2005) 475–479. [9] D. Karasik, N.A. Shimabuku, Y. Zhou, Y. Zhang, L.A. Cupples, D.P. Kiel, S. Demissie, A genome wide linkage scan of metacarpal size and geometry in the Framingham Study, Am. J. Hum. Biol. 20 (2008) 663–670. [10] E.A. Streeten, T.J. Beck, J.R. O'Connell, E. Rampersand, D.J. McBride, S.L. Takala, T.I. Pollin, K. Uusi-Rasi, B.D. Mitchell, A.R. Shuldiner, Autosome-wide linkage analysis of hip structural phenotypes in the Old Order Amish, Bone 43 (2008) 607–612. [11] H. Shen, J.R. Long, D.H. Xiong, Y.F. Guo, P. Xiao, Y.Z. Liu, L.J. Zhao, Y.J. Liu, H.Y. Deng, J.L. Li, R.R. Recker, H.W. Deng, A genomewide scan for quantitative trait loci underlying areal bone size variation in 451 Caucasian families, J. Med. Genet. 43 (2006) 873–880. [12] H.W. Deng, F.H. Xu, Q.Y. Huang, H. Shen, H. Deng, T. Conway, Y.J. Liu, Y.Z. Liu, J.L. Li, H.T. Zhang, K.M. Davies, R.R. Recker, A whole-genome linkage scan suggests several genomic regions potentially containing quantitative trait loci for osteoporosis, J. Clin. Endocrinol. Metab. 87 (2002) 5151–5159. [13] H.W. Deng, F.H. Xu, Y.Z. Liu, H. Shen, H. Deng, Q.Y. Huang, Y.J. Liu, T. Conway, J.L. Li, K.M. Davies, R.R. Recker, A whole-genome linkage scan suggests several genomic regions potentially containing QTLs underlying the variation of stature, Am. J. Med. Genet. 113 (2002) 29–39. [14] Q.Y. Huang, F.H. Xu, H. Shen, H.Y. Deng, T. Conway, Y.J. Liu, Y.Z. Liu, J.L. Li, M.X. Li, K.M. Davies, R.R. Recker, H.W. Deng, Genome scan for QTLs underlying bone size variation at 10 refined skeletal sites: genetic heterogeneity and the significance of phenotype refinement, Physiol. Genomics 17 (2004) 326–331. [15] X. Li, G. Masinde, W. Gu, J. Wergedal, S. Mohan, D.J. Baylink, Genetic dissection of femur breaking strength in a large population (MRL/MpJ × SJL/J) of F2 Mice: single QTL effects, epistasis, and pleiotropy, Genomics 79 (2002) 734–740. [16] X. Li, G. Masinde, W. Gu, J. Wergedal, M. Hamilton-Ulland, S. Xu, S. Mohan, D.J. Baylink, Chromosomal regions harboring genes for the work to femur failure in mice, Funct. Integr. Genomics 1 (2002) (2002) 367–374. [17] R.F. Klein, R.J. Turner, L.D. Skinner, K.A. Vartanian, M. Serang, A.S. Carlos, M. Shea, J.K. Belknap, E.S. Orwoll, Mapping quantitative trait loci that influence femoral cross-sectional area in mice, J. Bone Miner. Res. 17 (2002) 1752–1760. [18] I. Alam, Q. Sun, L. Liu, D.L. Koller, L.G. Carr, M.J. Econs, T. Foroud, C.H. Turner, Sex-specific genetic loci for femoral neck bone mass and strength identified in inbred COP and DA rats, J. Bone. Miner. Res. 23 (2008) 850–859. [19] Q. Sun, I. Alam, L. Liu, D.L. Koller, L.G. Carr, M.J. Econs, T. Foroud, C.H. Turner, Genetic loci affecting bone structure and strength in inbred COP and DA rats, Bone 42 (2008) 547–553. [20] I. Alam, Q. Sun, L. Liu, D.L. Koller, T. Fishburn, L.G. Carr, M.J. Econs, T. Foroud, C.H. Turner, Whole-genome scan for linkage to bone strength and structure in inbred Fischer 344 and Lewis rats, J. Bone Miner. Res. 20 (2005) 1589–1596. [21] I. Alam, Q. Sun, L. Liu, D.L. Koller, T. Fishburn, L.G. Carr, M.J. Econs, T. Foroud, C.H. Turner, Identification of a quantitative trait locus on rat chromosome 4 that is strongly linked to femoral neck structure and strength, Bone 39 (2006) 93–99. [22] C.H. Turner, R.K. Roeder, A. Wieczorek, T. Foroud, G. Liu, M. Peacock, Variability in skeletal mass, structure, and biomechanical properties among inbred strains of rats, J. Bone Miner. Res. 16 (2001) 1532–1539. [23] I. Alam, Q. Sun, L. Liu, D.L. Koller, L. Liu, H.J. Edenberg, M.J. Econs, T. Foroud, C.H. Turner, Genomic expression analysis of rat chromosome 4 for skeletal traits at femoral neck, Physiol. Genomics 35 (2008) 191–196. [24] G.K. Smyth, Linear models and empirical Bayes methods for assessing differential expression in microarray experiments, Stat. Appl. Genet. Mol. Biol. 3 (2004) Article3. [25] J.N. McClintick, H.J. Edenberg, Effects of filtering by present call on analysis of microarray experiments, BMC. Bioinformatics 7 (2006) 49. [26] Y. Benjamini, Y. Hochberg, Controlling the false discovery rate: a practical and powerful approach to multiple testing, J. R. Stat. Soc. 57 (1995) 289–290. [27] J. Batut, B. Schmierer, J. Cao, L.A. Raftery, C.S. Hill, M. Howell M, Two highly related regulatory subunits of PP2A exert opposite effects on TGF-beta/Activin/Nodal signaling, Development 135 (2008) 2927–2937. [28] C. Stehlik, L. Fiorentino, A. Dorfleutner, J.M. Bruey, E.M. Ariza, J. Sagara, J.C. Reed, The PAAD/PYRIN-family protein ASC is a dual regulator of a conserved step in nuclear factor kappaB activation pathways, J. Exp. Med. 196 (2002) 1605–1615. [29] D.L. Koller, M.J. Econs, P.A. Morin, J.C. Christian, S.L. Hui, P. Parry, M.E. Curran, L.A. Rodriguez, P.M. Conneally, G. Joslyn, M. Peacock, C.C. Johnston, T. Foroud, Genome screen for QTLs contributing to normal variation in bone mineral density and osteoporosis, J. Clin. Endocrinol. Metab. 85 (2000) 3116–3120. [30] H. Takayanagi, K. Sato, A. Takaoka, T. Taniguchi, Interplay between interferon and other cytokine systems in bone metabolism, Immunol. Rev. 208 (2005) 181–193. [31] G. Gdynia, J. Lehmann-Koch, S. Sieber, K.E. Tagscherer, A. Fassl, H. Zentgraf, S. Matsuzawa, J.C. Reed, W. Roth, BLOC1S2 interacts with the HIPPI protein and sensitizes NCH89 glioblastoma cells to apoptosis, Apoptosis 13 (2008) 437–447.
262
I. Alam et al. / Genomics 94 (2009) 257–262
[32] H. Ohno, Physiological roles of clathrin adaptor AP complexes: lessons from mutant animals, J. Biochem. 139 (2006) 943–948. [33] J.T. Blitzer, R. Nusse, A critical role for endocytosis in Wnt signaling, BMC. Cell. Biol. 7 (2006) 28. [34] T. Mizoguchi, K. Nakajima, K. Hatsuzawa, M. Nagahama, H.P. Hauri, M. Tagaya, K. Tani, Determination of functional regions of p125, a novel mammalian Sec23p-interacting protein, Biochem. Biophys. Res. Commun. 279 (2000) 144–149. [35] S.E. Calvano, W. Xiao, D.R. Richards, R.M. Felciano, H.V. Baker, R.J. Cho, R.O. Chen, B.H. Brownstein, J.P. Cobb, S.K. Tschoeke, C. Miller-Graziano, L.L. Moldawer, M.N. Mindrinos, R.W. Davis, R.G. Tompkins, S.F. Lowry, A networkbased analysis of systemic inflammation in humans, Nature 437 (2005) 1032–1037. [36] X. Li, H. Liu, L. Qin, J. Tamasi, M. Bergenstock, S. Shapses, J.H. Feyen, D.A. Notterman, N.C. Partridge, Determination of dual effects of parathyroid hormone on skeletal gene expression in vivo by microarray and network analysis, J. Biol. Chem. 282 (2007) 33086–33097. [37] M.P. Lynch, C. Capparelli, J.L. Stein, G.S. Stein, J.B. Lian, Apoptosis during bone-like tissue development in vitro, J. Cell. Biochem. 68 (1998) 31–49. [38] H. Takayanagi, S. Kim, K. Matsuo, H. Suzuki, T. Suzuki, K. Sato, T. Yokochi, H. Oda, K. Nakamura, N. Ida, E.F. Wagner, T. Taniguchi, RANKL maintains bone home-
[39]
[40]
[41]
[42]
[43]
[44]
ostasis through c-Fos-dependent induction of interferon-beta, Nature 416 (2002) 744–749. A.K. Townley, Y. Feng, K. Schmidt, D.A. Carter, R. Porter, P. Verkade, D.J. Stephens, Efficient coupling of Sec23–Sec24 to Sec13–Sec31 drives COPII-dependent collagen secretion and is essential for normal craniofacial development, J. Cell. Sci. 121 (2008) 3025–3034. A.E. Kearns, S. Khosla, P.J. Kostenuik, Receptor activator of nuclear factor kappaB ligand and osteoprotegerin regulation of bone remodeling in health and disease, Endocr. Rev. 29 (2008) 155–192. A. Gingery, E. Bradley, A. Shaw, M.J. Oursler, Phosphatidylinositol 3-kinase coordinately activates the MEK/ERK and AKT/NFkappaB pathways to maintain osteoclast survival, J. Cell. Biochem. 89 (2003) 165–179. C. Ge, G. Xiao, D. Jiang, R.T. Franceschi, Critical role of the extracellular signalregulated kinase-MAPK pathway in osteoblast differentiation and skeletal development, J. Cell. Biol. 176 (2007) 709–718. S.E. Lee, K.M. Woo, S.Y. Kim, H.M. Kim, K. Kwack, Z.H. Lee, H.H. Kim, The phosphatidylinositol 3-kinase, p38, and extracellular signal-regulated kinase pathways are involved in osteoclast differentiation, Bone 30 (2002) 71–77. M.C. Horowitz, Y. Xi, K. Wilson, M.A. Kacena, Control of osteoclastogenesis and bone resorption by members of the TNF family of receptors and ligands, Cytokine. Growth. Factor. Rev. 12 (2001) 9–18.
Contents lists available at ScienceDirect
Genomics j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y g e n o
Differentially expressed genes strongly correlated with femur strength in rats Imranul Alam a, Qiwei Sun a, Daniel L. Koller b, Lixiang Liu b, Yunlong Liu c, Howard J. Edenberg d, Jiliang Li e, Tatiana Foroud b, Charles H. Turner a,⁎ a
Department of Biomedical Engineering, Indiana University Purdue University Indianapolis (IUPUI), 1120 South Drive, Fesler Hall 115, Indianapolis, IN 46202-5251, USA Department of Medical and Molecular Genetics, Indiana University Purdue University Indianapolis (IUPUI), IN, USA Department of Medicine, Indiana University Purdue University Indianapolis (IUPUI), IN, USA d Department of Biochemistry and Molecular Biology, Indiana University Purdue University Indianapolis (IUPUI), IN, USA e Department of Biology, Indiana University Purdue University Indianapolis (IUPUI), IN, USA b c
a r t i c l e
i n f o
Article history: Received 30 January 2009 Accepted 25 May 2009 Available online 29 May 2009 Keywords: Femur strength Gene expression Microarray QTLs Osteoporotic fracture
a b s t r a c t The region of chromosome 1q33–q54 harbors quantitative trait loci (QTL) for femur strength in COP × DA and F344 × LEW F2 rats. The purpose of this study is to identify the genes within this QTL region that contribute to the variation in femur strength. Microarray analysis was performed using RNA extracted from femurs of COP, DA, F344 and LEW rats. Genes differentially expressed in the 1q33–q54 region among these rat strains were then ranked based on the strength of correlation with femur strength in F2 animals derived from these rats. A total of 214 genes in this QTL region were differentially expressed among all rat strains, and 81 genes were found to be strongly correlated (r2 N 0.50) with femur strength. Of these, 12 candidate genes were prioritized for further validation, and 8 of these genes (Ifit3, Ppp2r5b, Irf7, Mpeg1, Bloc1s2, Pycard, Sec23ip, and Hps6) were confirmed by quantitative PCR (qPCR). Ingenuity Pathway Analysis suggested that these genes were involved in interferon alpha, nuclear factor-kappa B (NFkB), extracellular signal-related kinase (ERK), hepatocyte nuclear factor 4 alpha (HNF4A) and tumor necrosis factor (TNF) pathways. © 2009 Published by Elsevier Inc.
Introduction Osteoporosis is a common polygenic disorder characterized by reduced bone mineral density (BMD) and compromised bone quality with microarchitectural deterioration leading to an increased susceptibility to fracture at multiple skeletal sites [1]. BMD, structure and strength are the primary skeletal determinants of osteoporotic fracture risk [2–4] and are under substantial genetic control [5–8]. Most of the studies identifying osteoporosis genes primarily focus on BMD as it is the most commonly available tool for predicting fracture risk. However, BMD is not an entirely adequate measure of bone strength and structure. Bone strength, the phenotype that best reflects skeletal fragility also depends on bone geometry (size and shape), microarchitecture and material properties. Bone strength is under the control of multiple genes, each of which has a modest contribution to the genetic variance [15,20]. Identification of genes underlying bone strength will reveal valuable insight regarding the genetics of osteoporosis and fracture risk. Over many years, several linkage studies in human and experimental animal models have identified quantitative trait loci (QTLs) linked to bone structure and strength [9–21], but the genes underlying these phenotypes are largely unknown. Identification of candidate
⁎ Corresponding author. Fax: +1 317 278 9568. E-mail address: [email protected] (C.H. Turner). 0888-7543/$ – see front matter © 2009 Published by Elsevier Inc. doi:10.1016/j.ygeno.2009.05.008
genes within the broad QTL intervals, which typically harbor hundreds of genes, is extremely challenging and requires multiple complementary approaches. One strategy is the development of a congenic strain, which can narrow the critical QTL to a small chromosomal segment harboring fewer candidate genes. However, due to the extensive breeding that is required, this approach is quite time-consuming and expensive. Other approaches include creation of recombinant inbred lines and haplotype-based analyses. However, these strategies do not produce mapping with high enough resolution to detect the causal genes. An alternative strategy involves an integrative genetic approach combining linkage with genomic expression analysis in multiple inbred lines already used for QTL analysis. Such studies are extremely difficult in humans, but they can be undertaken in animal models. Previously, we showed that skeletal mass, structure and strength vary among inbred strains of rats in a site-specific manner [22]. We demonstrated that despite similar body size, substantial variation exists in bone geometry and biomechanical properties between adult Copenhagen 2331 (COP), Dark Agouti (DA), Fischer 344 (F344) and Lewis (LEW) rats [22]. Subsequently, we generated a large number of second filial (F2) female progeny derived from (COP × DA) and (F344 × LEW) crosses, and detected significant linkage for femur ultimate force (Fu) (a biomechanical measurement that determines the maximum load that the bone can support before failing) within the region of 1q33–q54 in both of these crosses [19,20]. Thus, 1q33–q54 region harbors QTLs for a bone phenotype directly related to femur strength.
258
I. Alam et al. / Genomics 94 (2009) 257–262
The purpose of this study was to narrow the list of candidate genes within the region of 1q33–q54 contributing to the variation in femur strength among COP, DA, F344 and LEW rats. We chose these strains of rats because, at present, these are the only available rat models in which bone QTLs have been identified [18–21]. We used the data from microarray-based gene expression analyses to identify genes differentially expressed among these rats. These genes were then ranked based on their correlation with Fu in F2 animals derived from these rat strains. The expression levels of several strongly correlated (r2 N 0.9) candidate genes were confirmed using quantitative real-time PCR (qPCR). In addition, we applied Ingenuity Pathway Analysis (IPA) for analysis of the related pathways among prioritized candidate genes based on their molecular function, biological process and cellular component of the gene products. IPA is a structured, network-based system to analyze genome-wide gene expression in the context of known functional interrelationships among genes, proteins, and small molecules. This approach has been successfully used for various systems including bone [35,36].
Microarray analyses and correlation with femur phenotypes The Affymetrix microarray analyses showed that a total of 214 genes residing in the region of 1q33–q54 (among the 999 total genes on the microarray probe sets) were differentially expressed among COP, DA, F344 and LEW rats with a FDR less than 10%. Regression analyses were performed for all 214 differentially expressed genes to assess the correlation between gene expression and the femur Fu among F2 animals homozygous for that strain's allele at the markers (D1Rat69 for COP × DA; D1Rat291 for F344 × LEW) underlying the QTL. The 1q33–q54 QTL for Fu was replicated in both (COP × DA) and (F344 × LEW) crosses; therefore, we identified candidate genes that showed strong correlation (r2 N 0.50) for this phenotype in these rats. We considered any r-value of above 0.7, which corresponds to a r2-value of above 0.49, as a strong correlation. A total of 81 genes, including 72 candidate genes and 9 predicted genes (Supplementary Table 1), were detected using this method. Among these 72 genes, 12 were ranked as the highest priority because they were found to explain more than 90% (r2 N 0.90) of the variation of femur Fu (Supplementary Table 1).
Results qPCR analyses Genetic loci for femur Fu on 1q33–q54 Evidence of significant linkage of femur Fu was detected in the region of 1q33–q54 in COP × DA F2 animals with LOD scores of 6.7 (Fig. 1) [19]. In addition, significant evidence of linkage of femur Fu was observed in F344 × LEW F2 rats with a LOD score of 4.8, in a broad region that overlapped the COP × DA peak (Fig. 1) [20]. Femur Fu as a function of genotypes Genotypic means for femur Fu was measured on Chr 1 at D1Rat69 (101.5 cM) in the COP × DA F2 rats and at D1Rat291 (101.4 cM) in the F344 × LEW F2 rats. Since a QTL was detected at marker D1Rat69 in the COP × DA cross and at marker D1Rat291 in the F344 × LEW cross, the mean values for the femur Fu were significantly different between the F2 animals homozygous for the COP (c/c) and the DA (d/d) F2 alleles (120.35 ± 1.81 vs. 129.08 ± 1.87; p = 0.001) and the F2 animals homozygous for the F344 (f/f) and LEW (l/l) F2 alleles (130.47 ± 1.11 vs. 138.02 ± 1.08; p b 0.0001).
qPCR analyses of the 12 prioritized candidate genes in Supplementary Table 1 confirmed the strong correlation between gene expression and femur ultimate force (Table 1). The correlation coefficients for Pold4, Ifit3, Hps6, Mrpl3, Map4k2, Sec23ip, Rras2, Pycard, Ppp2r5b, Irf7, Mpeg1 and Bloc1s2 between microarray and qPCR analyses were 0.65, 0.99, 0.70, 0.70, 0.24, 0.81, 0.22, 0.94, 0.96, 0.99, 0.94 and 0.89, respectively, indicating good concordance between the two methods for all genes except Map4k2 and Rras2. Among these 12 candidate genes, 8 genes (Ifit3, Ppp2r5b, Irf7, Mpeg1, Bloc1s2, Pycard, Sec23ip and Hps6) were confirmed to have similar correlations with femur strength by qPCR (r2-values ranging from 0.59 to 0.97) (Table 1). Pathway analysis The 8 candidate genes that were confirmed by qPCR were mapped to pathways using Ingenuity Pathway Analysis. These genes were directly or indirectly connected to interferon alpha, nuclear factorkappa B (NFkB), extracellular signal-related kinase (ERK), hepatocyte nuclear factor 4 alpha (HNF4A) and tumor necrosis factor (TNF) pathways (Fig. 2). Discussion
Fig. 1. Linkage analysis for femur ultimate force (Fu) for rat chromosome 1 in COP × DA and F344 × LEW crosses in the region of 1q33–q54 (D1Rat52–D1Rat82). LOD scores plotted on the y axis vs. the relative microsatellite markers located along chromosome 1 on the x axis.
Genetic mapping for complex traits usually identifies chromosomal regions that are typically quite broad in nature. The 1q33–q54 QTL region linked to various femur phenotypes also encompasses a broad region of roughly 60 cM of the rat genome, harboring nearly 1000 genes. In addition, many of these genes might have pleiotropic effects on different femur phenotypes. Also, there might be a cluster of genes contributing individually to these phenotypes. To prioritize candidate genes for further analysis, we employed an integrative genetic approach combining linkage, gene expression and physical traits. We took the advantage of microarray-based expression analyses to simultaneously screen thousands of genes for their differential expression in COP, DA, F344 and LEW rats, as these are the only available rat models in which bone QTLs have been identified [18–21]. In addition, we used 4-week-old rats for gene expression analyses since we wanted to assess gene regulation of bone accrual and to identify genes contribute to the accumulation of peak bone mass. By comparing the genes differentially expressed among these rats, we were able to significantly reduce the number of candidate genes in the 1q33–q54 QTL region to 214.
I. Alam et al. / Genomics 94 (2009) 257–262 Table 1 r-square values from qPCR analysis of top 12 candidate genes for femur ultimate force within 1q33–q54 region on Chr 1 in COP, DA, F344 and LEW ratsa. Gene symbol
Gene name
Rat Femur ultimate location force (r2)
Human synteny
Ifit3
Interferon-induced protein with tetratricopeptide repeats 3 Protein phosphatase 2, regulatory subunit B, delta isoform Interferon regulatory factor 7 Macrophage expressed gene 1 Biogenesis of lysosome-related organelles complex-1, subunit 2 PYD and CARD domain containing SEC23 interacting protein Hermansky–Pudlak syndrome 6 Polymerase (DNA-directed), delta 4 Mitochondrial ribosomal protein L43 Related RAS viral (r-ras) oncogene homolog 2 Mitogen activated protein kinase kinase kinase kinase 2
0.97
1q52
10q23.3
0.85
1q43
11q13
0.88 0.72 0.74
1q41 1q43 1q54
11p15.5 11q12.1 10q24.31
0.65 0.61 0.59 0.45 0.34 0.24
1q36 1q37 1q54 1q42 1q54 1q34
16p12 10q26.11 10q24.3 11q13.1 10q24.31 11p15.2
0.20
1q43
11q13
Ppp2r5b Irf7 Mpeg1 Bloc1s2 Pycard Sec23ip Hps6 Pold4 Mrpl43 Rras2 Map4k2 a
Confirmed genes (r2 N 0.50) are indicated in bold face.
Prioritization of microarray-based candidate genes is usually done by ranking the genes for the magnitude of expression differences (fold-differences) between different strains. However, with complex traits like osteoporosis, even subtle changes in gene expression could be important. Thus, stronger evidence for the association between a gene and a trait might be identified by examining the correlation between gene expression and a given physical trait. Using this
259
approach, among the 214 candidate genes, we were able to prioritize 72 genes that were strongly correlated (r2 N 0.5) with femur strength. Furthermore, by investigating genes those were correlated with more than 90% variation (r2 N 0.9) with femur ultimate force, we further narrowed the list of genes to 12 candidate (Supplementary Table 1). Among these 12 candidate genes, 8 genes (Ifit3, Ppp2r5b, Irf7, Mpeg1, Bloc1s2, Pycard, Sec23ip and Hps6) were confirmed to have similar strong correlations (r2-values ranging from 0.59 to 0.97) with the femur strength when expression was measured by qPCR (Table 1). Surprisingly, four genes (Pold4, Mrpl43, Map4k2 and Rras2) that showed high correlations for femur strength by Affymetrix analysis were not confirmed by qPCR; the difference may relate to probe location for the arrays at the 3′ end of the target gene. The 8 candidate genes that were confirmed by qPCR were associated with pathways directly or indirectly involved in controlling cell morphology and development, cellular assembly and organization, cell-to-cell signaling, molecular transport, protein trafficking, and cellular function and maintenance (Fig. 2). Genes in the canonical pathways (networks generated using IPA) were related to interferon signaling and activation of pattern recognition receptors. When these genes were categorized based on location of cellular components Ifit3, Mpeg1, Bloc1s2, Pycard, Sec23ip and Hps6 are located in cytoplasm while Ppp2r5b and Irf7 are located in the nucleus. Recently, signaling networks with interferon and other cytokines have been shown to be important in bone metabolism [30]. Interestingly, the top 2 candidate genes (Ifit3 and Irf7) identified in this study are interferon-related genes and pathway analysis revealed several interferon molecules involved in the same network (Fig. 2). Interferon regulatory factors Irf1 and Irf2 have been associated with osteoblast apoptosis [37] but roles for Ifit3 and Irf7 in bone biology have not yet been studied. However, both
Fig. 2. Network of the 8 candidate genes that were confirmed by qPCR (highlighted in red) for femur strength in Ingenuity Pathway Analysis (IPA). The well-known NFkB, ERK and TNF pathways related to bone metabolism is highlighted in blue.
260
I. Alam et al. / Genomics 94 (2009) 257–262
interferon alpha and beta signaling has been shown to be important for maintaining normal bone mass by regulating osteoclastic bone resorption [38]. Interferon receptor alpha 1 knockout mice exhibit an osteoporosis phenotype due to increased osteoclastogenesis and decreased trabecular bone mass. Interferon-beta inhibits RANK-ligand mediated osteoclastogenesis by inhibiting the expression of c-FOS and mice lacking interferon-beta have severe osteopenia due to enhanced osteoclastogenesis [38]. In addition, Ppp2r5b (Ppp2r2d) has been previously shown to have role in TGF-beta signaling [27] and Pycard is a regulator of NFkB activation pathway [28]. Hps6 is associated with Hermansky–Pudlak syndrome, an autosomal recessive syndrome due to a mutation in clathrin-associated adaptor protein [32]. Because clathrin is important in the WNT signaling pathway in bone [33], Hps6 might be involved in the molecular regulation of the WNT pathway. Bloc1s2 has been shown to regulate apoptosis in cancer [31], and Sec23ip is involved in protein export from the endoplasmic reticulum [34], and both of these genes are indirectly related to Hps5/6 genes. In a recent study, Sec23ip has been shown to be essential for normal craniofacial development [39]. The candidate genes and the interconnected molecules identified in this study are thus related to pathways involving both bone formation and resorption: nuclear factor-kappa B (NFkB) regulates osteoclast mediated bone resorption [40,41]; extracellular signal-related kinase (ERK) has critical role in osteoblast differentiation and bone development [42] as well as osteoclast differentiation and survival [43]; members of tumor necrosis factor (TNF) family of ligands and receptors are critical regulators of osteoclastogenesis [44,45]. However, further functional analysis using gene targeting strategies such as knockout and siRNA will be necessary to fully characterize the roles of these candidate genes in bone metabolism. The chromosomal region of the candidate genes for femur strength on Chr 1 in rats is homologous to regions of mouse chromosomes 7, 19 and human chromosomes 10q, 11p, 11q, and 16p (Table 1). Linkage to mouse Chr 7 was previously reported for femur strength and structure [16–17]. Furthermore, the homologous regions in human 10q21, 10q26 and 11p15 were previously linked to wrist bone size, hip BMD and spine bone size, respectively [12,13,29]. In this study, using an integrative genetic approach, we identified several top candidate genes underlying the QTL region of 1q33–q54 that were differentially expressed and strongly correlated with femur strength in rats. The approach we took here is complementary to conventional strategies for the dissection of complex traits and will facilitate the QTL to gene discovery process. However our approach has several limitations. Gene expression study based on microarray analysis explains only transcriptional regulation of genes and does not capture effects of alternative gene splicing or post-translational modification of proteins. In addition, the whole femur RNA came from a variety of different types of bone cells including bone marrow and thus we cannot be sure that the gene expression changes originated only from bone cells. Further studies involving different cell lines are thus necessary to identify which specific cell type expresses the genes identified in this study. Also, our results could be influenced by the genetic contributions from other QTLs from the same or other chromosomes. Some of these problems could be solved by isolating the Chr 1 QTL region from the effects of other QTLs throughout the genome by using congenic rats. Also, whether the same candidate genes regulate femur strength in male rats remains to be determined, as we only studied female rats in this study. Materials and methods In previous studies, we created 423 and 595 female 2nd filial (F2) offspring from COP × DA and F344 × LEW crosses, respectively. The rats were allowed to grow to 6 months (26 weeks) of age before sacrifice. The lower limbs and spleen were dissected out. The right side of the
lower limbs were immediately stored at −20 °C for biomechanical analysis. The spleens were immediately frozen in liquid nitrogen before transferring them to −80 °C. The right femurs were tested in three-point bending by positioning them on the lower supports of a three-point bending fixture and applying load at the midpoint using a material testing machine (Alliance RT/5, MTS Systems Corp., Eden Prairie, MN, USA) and ultimate force (Fu; N) was obtained from the force vs. displacement curves in TestWorks software, version 4.06 as described previously [19,20]. Genotyping for each animal was accomplished with microsatellite markers using genomic DNA isolated from spleen previously shown to be polymorphic for COP, DA, F344 and LEW rats with an average interval of 20 cM from centromere to telomere. Chromosomal positions, marker order, and map positions were obtained from the Rat Genome Database (RGD) website (http://rgd.mcw. edu/). Genotypic data in the F2 animals for each marker were tested to ensure the expected Mendelian ratios. Marker maps were generated using MAPMARKER/EXP with data from all available genotyped F2 animals. Multipoint quantitative linkage analysis was then performed to identify chromosomal regions that harbor QTL for femur Fu using the programs R/qtl and QTL/Cartographer. Permutation testing was employed to obtain appropriate genome-wide significance levels for the linkage results and ANOVA was used for the most significant marker in each QTL region to further characterize significant genotypic group differences as described previously [19,20]. RNA extraction and microarray analysis Whole femora were harvested from 4-week-old COP, DA, F344 and LEW animals and were immediately frozen in liquid nitrogen and stored at − 80 °C until required. RNA from femurs was extracted from 4 animals per strain using Trizol (Invitrogen, Carlsbad, CA), followed by further purification using an RNeasy Mini Kit (Qiagen Inc., Valencia, CA) as described previously [23]. RNA quality was determined using a 2100 Bioanalyzer (Agilent, Palo Alto, CA) and was quantified using a spectrophotometer (NanoDrop, Wilmington, DE, USA). For microarray analysis, 5 μg of total RNA from each sample was labeled and hybridized to Affymetrix Rat Genome 230 2.0 GeneChips (Affymetrix, Santa Clara, CA) as described previously [23]. Quality control for RNA and Affymetrix data was done as described previously [23]. Microarray data analysis and informatics The images from each array were analyzed using Affymetrix GeneChip Operating system (GCOS) 1.2 software. Mapping of probe sets to chromosomal location was accomplished with data provided by Affymetrix. The region of 1q33–q54, limited by the markers of D1Rat52 (80.1 cM) and D1Rat82 (139.0 cM), harbors a total of 999 genes. Among these genes, the identities of the 214 differentially expressed genes within this QTL region 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 [25]. FDR was calculated by the method of Benjamini and Hochberg [26]. Probe sets were considered differentially expressed if the FDR was less than 10%. We used 10% FDR because we wanted to reduce the number of false negatives when we identified differentially expressed genes. The microarray data set was submitted to the NCBI Gene Expression Omnibus (GEO) Express web portal (GEO accession number GSE 11180).
I. Alam et al. / Genomics 94 (2009) 257–262
Quantitative real-time PCR (qPCR) analysis Five micrograms of total RNA (the same RNA used for Affymetrix analysis) was reverse transcribed using Superscript III reverse transcription reagent for first-strand cDNA synthesis (Invitrogen, CA). Twelve strongly correlated (r2 N 0.90) candidate genes for femur Fu were selected for verification of microarray data by qPCR analysis. All real-time PCR reactions contained the first-strand cDNA corresponding to 25 ng of total RNA. Real-time detection of PCR products was performed using ABI PRISM 7300 sequence detector (Applied Biosystem, CA) as described previously [23]. Expression of mRNA was calculated based on a relative standard curve and normalized to beta-actin. All qPCR analysis used triplicates of each of 4 biological samples. Statistics Differential expression of each gene among all strains was calculated by ANOVA using the package Limma [24]. Regression analysis was performed with the average gene expression level for the strain as the dependent variable, and the phenotypic mean value of femur Fu 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 using the statistical software package StatView (Abacus Concepts, Inc., Berkeley, CA). Pathway analysis The interactions between differentially expressed genes confirmed by qPCR for femur strength and all other genes were investigated using Ingenuity Pathway Analysis (IPA 5.0; Ingenuity Systems, Inc., Mountain View, CA). The differentially expressed genes 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. Acknowledgments This work was supported by the US National Institutes of Health through the following grants: RO1AR047822, P01AG018397 and AA10707. 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.). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.ygeno.2009.05.008. References [1] J.A. Kanis, L.J. Melton, C. Christiansen, C.C. Johnston, N. Khaltaev, The diagnosis of osteoporosis, J. Bone. Miner. Res. 9 (1994) 1137–1141. [2] S.R. Cummings, D.M. Black, M.C. Nevitt, W. Browner, J. Cauley, K. Ensrud, H.K. Genant, L. Palermo, J. Scott, T.M. Vogt, Bone density at various sites for prediction of hip fractures, The Study of Osteoporotic Fractures Research Group, Lancet 341 (1993) 72–75. [3] K.G. Faulkner, S.R. Cummings, D. Black, L. Palermo, C.C. Gluer, H.K. Genant, Simple measurement of femoral geometry predicts hip fracture: the study of osteoporotic fractures, J. Bone. Miner. Res. 8 (1993) 1211–1217.
261
[4] M. Peacock, C.H. Turner, G. Liu, A.K. Manatunga, L. Timmerman, C.C. Johnston, Better discrimination of hip fracture using bone density, geometry and architecture, Osteoporos. Int. 5 (1995) 167–173. [5] N.K. Arden, J. Baker, C. Hogg, K. Bann, T.D. Spector, The heritability of bone mineral density, ultrasound of the calcaneus and hip axis length: a study of postmenopausal twins, J. Bone Miner. Res. 11 (1996) 530–534. [6] P. Garnero, N.K. Arden, G. Griffiths, P.D. Delmas, T.D. Spector, Genetic influence on bone turnover in postmenopausal twins, J. Clin. Endocrinol. Metab. 81 (1996) 140–146. [7] M. Peacock, C.H. Turner, M.J. Econs, T. Foroud, Genetics of osteoporosis, Endocr. Rev. 23 (2002) 378–383. [8] S.H. Ralston, Genetic determinants of osteoporosis, Curr. Opin. Rheumatol. 17 (2005) 475–479. [9] D. Karasik, N.A. Shimabuku, Y. Zhou, Y. Zhang, L.A. Cupples, D.P. Kiel, S. Demissie, A genome wide linkage scan of metacarpal size and geometry in the Framingham Study, Am. J. Hum. Biol. 20 (2008) 663–670. [10] E.A. Streeten, T.J. Beck, J.R. O'Connell, E. Rampersand, D.J. McBride, S.L. Takala, T.I. Pollin, K. Uusi-Rasi, B.D. Mitchell, A.R. Shuldiner, Autosome-wide linkage analysis of hip structural phenotypes in the Old Order Amish, Bone 43 (2008) 607–612. [11] H. Shen, J.R. Long, D.H. Xiong, Y.F. Guo, P. Xiao, Y.Z. Liu, L.J. Zhao, Y.J. Liu, H.Y. Deng, J.L. Li, R.R. Recker, H.W. Deng, A genomewide scan for quantitative trait loci underlying areal bone size variation in 451 Caucasian families, J. Med. Genet. 43 (2006) 873–880. [12] H.W. Deng, F.H. Xu, Q.Y. Huang, H. Shen, H. Deng, T. Conway, Y.J. Liu, Y.Z. Liu, J.L. Li, H.T. Zhang, K.M. Davies, R.R. Recker, A whole-genome linkage scan suggests several genomic regions potentially containing quantitative trait loci for osteoporosis, J. Clin. Endocrinol. Metab. 87 (2002) 5151–5159. [13] H.W. Deng, F.H. Xu, Y.Z. Liu, H. Shen, H. Deng, Q.Y. Huang, Y.J. Liu, T. Conway, J.L. Li, K.M. Davies, R.R. Recker, A whole-genome linkage scan suggests several genomic regions potentially containing QTLs underlying the variation of stature, Am. J. Med. Genet. 113 (2002) 29–39. [14] Q.Y. Huang, F.H. Xu, H. Shen, H.Y. Deng, T. Conway, Y.J. Liu, Y.Z. Liu, J.L. Li, M.X. Li, K.M. Davies, R.R. Recker, H.W. Deng, Genome scan for QTLs underlying bone size variation at 10 refined skeletal sites: genetic heterogeneity and the significance of phenotype refinement, Physiol. Genomics 17 (2004) 326–331. [15] X. Li, G. Masinde, W. Gu, J. Wergedal, S. Mohan, D.J. Baylink, Genetic dissection of femur breaking strength in a large population (MRL/MpJ × SJL/J) of F2 Mice: single QTL effects, epistasis, and pleiotropy, Genomics 79 (2002) 734–740. [16] X. Li, G. Masinde, W. Gu, J. Wergedal, M. Hamilton-Ulland, S. Xu, S. Mohan, D.J. Baylink, Chromosomal regions harboring genes for the work to femur failure in mice, Funct. Integr. Genomics 1 (2002) (2002) 367–374. [17] R.F. Klein, R.J. Turner, L.D. Skinner, K.A. Vartanian, M. Serang, A.S. Carlos, M. Shea, J.K. Belknap, E.S. Orwoll, Mapping quantitative trait loci that influence femoral cross-sectional area in mice, J. Bone Miner. Res. 17 (2002) 1752–1760. [18] I. Alam, Q. Sun, L. Liu, D.L. Koller, L.G. Carr, M.J. Econs, T. Foroud, C.H. Turner, Sex-specific genetic loci for femoral neck bone mass and strength identified in inbred COP and DA rats, J. Bone. Miner. Res. 23 (2008) 850–859. [19] Q. Sun, I. Alam, L. Liu, D.L. Koller, L.G. Carr, M.J. Econs, T. Foroud, C.H. Turner, Genetic loci affecting bone structure and strength in inbred COP and DA rats, Bone 42 (2008) 547–553. [20] I. Alam, Q. Sun, L. Liu, D.L. Koller, T. Fishburn, L.G. Carr, M.J. Econs, T. Foroud, C.H. Turner, Whole-genome scan for linkage to bone strength and structure in inbred Fischer 344 and Lewis rats, J. Bone Miner. Res. 20 (2005) 1589–1596. [21] I. Alam, Q. Sun, L. Liu, D.L. Koller, T. Fishburn, L.G. Carr, M.J. Econs, T. Foroud, C.H. Turner, Identification of a quantitative trait locus on rat chromosome 4 that is strongly linked to femoral neck structure and strength, Bone 39 (2006) 93–99. [22] C.H. Turner, R.K. Roeder, A. Wieczorek, T. Foroud, G. Liu, M. Peacock, Variability in skeletal mass, structure, and biomechanical properties among inbred strains of rats, J. Bone Miner. Res. 16 (2001) 1532–1539. [23] I. Alam, Q. Sun, L. Liu, D.L. Koller, L. Liu, H.J. Edenberg, M.J. Econs, T. Foroud, C.H. Turner, Genomic expression analysis of rat chromosome 4 for skeletal traits at femoral neck, Physiol. Genomics 35 (2008) 191–196. [24] G.K. Smyth, Linear models and empirical Bayes methods for assessing differential expression in microarray experiments, Stat. Appl. Genet. Mol. Biol. 3 (2004) Article3. [25] J.N. McClintick, H.J. Edenberg, Effects of filtering by present call on analysis of microarray experiments, BMC. Bioinformatics 7 (2006) 49. [26] Y. Benjamini, Y. Hochberg, Controlling the false discovery rate: a practical and powerful approach to multiple testing, J. R. Stat. Soc. 57 (1995) 289–290. [27] J. Batut, B. Schmierer, J. Cao, L.A. Raftery, C.S. Hill, M. Howell M, Two highly related regulatory subunits of PP2A exert opposite effects on TGF-beta/Activin/Nodal signaling, Development 135 (2008) 2927–2937. [28] C. Stehlik, L. Fiorentino, A. Dorfleutner, J.M. Bruey, E.M. Ariza, J. Sagara, J.C. Reed, The PAAD/PYRIN-family protein ASC is a dual regulator of a conserved step in nuclear factor kappaB activation pathways, J. Exp. Med. 196 (2002) 1605–1615. [29] D.L. Koller, M.J. Econs, P.A. Morin, J.C. Christian, S.L. Hui, P. Parry, M.E. Curran, L.A. Rodriguez, P.M. Conneally, G. Joslyn, M. Peacock, C.C. Johnston, T. Foroud, Genome screen for QTLs contributing to normal variation in bone mineral density and osteoporosis, J. Clin. Endocrinol. Metab. 85 (2000) 3116–3120. [30] H. Takayanagi, K. Sato, A. Takaoka, T. Taniguchi, Interplay between interferon and other cytokine systems in bone metabolism, Immunol. Rev. 208 (2005) 181–193. [31] G. Gdynia, J. Lehmann-Koch, S. Sieber, K.E. Tagscherer, A. Fassl, H. Zentgraf, S. Matsuzawa, J.C. Reed, W. Roth, BLOC1S2 interacts with the HIPPI protein and sensitizes NCH89 glioblastoma cells to apoptosis, Apoptosis 13 (2008) 437–447.
262
I. Alam et al. / Genomics 94 (2009) 257–262
[32] H. Ohno, Physiological roles of clathrin adaptor AP complexes: lessons from mutant animals, J. Biochem. 139 (2006) 943–948. [33] J.T. Blitzer, R. Nusse, A critical role for endocytosis in Wnt signaling, BMC. Cell. Biol. 7 (2006) 28. [34] T. Mizoguchi, K. Nakajima, K. Hatsuzawa, M. Nagahama, H.P. Hauri, M. Tagaya, K. Tani, Determination of functional regions of p125, a novel mammalian Sec23p-interacting protein, Biochem. Biophys. Res. Commun. 279 (2000) 144–149. [35] S.E. Calvano, W. Xiao, D.R. Richards, R.M. Felciano, H.V. Baker, R.J. Cho, R.O. Chen, B.H. Brownstein, J.P. Cobb, S.K. Tschoeke, C. Miller-Graziano, L.L. Moldawer, M.N. Mindrinos, R.W. Davis, R.G. Tompkins, S.F. Lowry, A networkbased analysis of systemic inflammation in humans, Nature 437 (2005) 1032–1037. [36] X. Li, H. Liu, L. Qin, J. Tamasi, M. Bergenstock, S. Shapses, J.H. Feyen, D.A. Notterman, N.C. Partridge, Determination of dual effects of parathyroid hormone on skeletal gene expression in vivo by microarray and network analysis, J. Biol. Chem. 282 (2007) 33086–33097. [37] M.P. Lynch, C. Capparelli, J.L. Stein, G.S. Stein, J.B. Lian, Apoptosis during bone-like tissue development in vitro, J. Cell. Biochem. 68 (1998) 31–49. [38] H. Takayanagi, S. Kim, K. Matsuo, H. Suzuki, T. Suzuki, K. Sato, T. Yokochi, H. Oda, K. Nakamura, N. Ida, E.F. Wagner, T. Taniguchi, RANKL maintains bone home-
[39]
[40]
[41]
[42]
[43]
[44]
ostasis through c-Fos-dependent induction of interferon-beta, Nature 416 (2002) 744–749. A.K. Townley, Y. Feng, K. Schmidt, D.A. Carter, R. Porter, P. Verkade, D.J. Stephens, Efficient coupling of Sec23–Sec24 to Sec13–Sec31 drives COPII-dependent collagen secretion and is essential for normal craniofacial development, J. Cell. Sci. 121 (2008) 3025–3034. A.E. Kearns, S. Khosla, P.J. Kostenuik, Receptor activator of nuclear factor kappaB ligand and osteoprotegerin regulation of bone remodeling in health and disease, Endocr. Rev. 29 (2008) 155–192. A. Gingery, E. Bradley, A. Shaw, M.J. Oursler, Phosphatidylinositol 3-kinase coordinately activates the MEK/ERK and AKT/NFkappaB pathways to maintain osteoclast survival, J. Cell. Biochem. 89 (2003) 165–179. C. Ge, G. Xiao, D. Jiang, R.T. Franceschi, Critical role of the extracellular signalregulated kinase-MAPK pathway in osteoblast differentiation and skeletal development, J. Cell. Biol. 176 (2007) 709–718. S.E. Lee, K.M. Woo, S.Y. Kim, H.M. Kim, K. Kwack, Z.H. Lee, H.H. Kim, The phosphatidylinositol 3-kinase, p38, and extracellular signal-regulated kinase pathways are involved in osteoclast differentiation, Bone 30 (2002) 71–77. M.C. Horowitz, Y. Xi, K. Wilson, M.A. Kacena, Control of osteoclastogenesis and bone resorption by members of the TNF family of receptors and ligands, Cytokine. Growth. Factor. Rev. 12 (2001) 9–18.