Gene Expression Profiling in Monocytes and SNP Association Suggest the Importance of the Gene for Osteoporosis in Both Chinese and Caucasians

Osteoporosis is characterized mainly by low bone mineral density (BMD). Many cytokines and chemokines have been related with bone metabolism. Monocytes in the immune system are important sources of cytokines and chemokines for bone metabolism. However, no study has investigated in vivo expression of a large number of various factors simultaneously in human monocytes underlying osteoporosis. This study explored the in vivo expression pattern of general cytokines, chemokines, and their receptor genes in human monocytes and validated the significant genes by qRT-PCR and genetic association analyses. Expression profilings were performed in monocyte samples from 26 Chinese and 20 Caucasian premenopausal women with discordant BMD. Genome-wide association analysis with BMD variation was conducted in 1000 unrelated Caucasians. We selected 168 cytokines, chemokines, osteoclast-related factors, and their receptor genes for analyses. Significantly, the signal transducer and activator of transcription 1 (STAT1) gene was upregulated in the low versus the high BMD groups in both Chinese and Caucasians. We also revealed a significant association of the STAT1 gene with BMD variation in the 1000 Caucasians. Thus we conclude that the STAT1 gene is important in human circulating monocytes in the etiology of osteoporosis. © 2010 American Society for Bone and Mineral Research.


Introduction
O steoporosis is mainly characterized by low bone mineral density (BMD). Genetic factors have important influences on BMD and osteoporosis. (1)(2)(3) Recent studies have shown that the immune system is strongly related to bone metabolism in terms of osteoimmunology. (4)(5)(6)(7)(8) Pathologic bone resorption was observed in immune system-related diseases such as autoimmune arthritis, periodontitis, Paget's disease, and bone tumors. (9) Monocytes, important cells in immune system, produce a wide variety of factors such as interleukin 1 (IL-1), IL-6, tumor necrosis factor (TNF), transforming growth fator beta (TGF-b), and 1,25-dihydroxyvitamin D 3 [1,25(OH) 2 D 3 ] . (10) These factors are involved in bone metabolism by regulating osteoclastic differentiation. Monocytes are also potential precursors of osteoclasts. (11,12) In vitro studies demonstrated that monocytes can differentiate into osteoclasts with bone resorption function. (13,14) However, it is unknown whether other factors and mechanisms to regulate these factors are important in the ability of monocytes to affect bone metabolism. To address these questions, scientists have screened the differential gene expressions in osteoclastogenic cells using a high-throughput microarray platform. (15,16) Microarray technology has been used successfully for detection of gene expression profiles in diseases such as inflammatory breast cancer and urinary bladder cancer. (17,18) Theoretical studies also supported the reliability of using a microarray platform for the quantitative characterization of gene expression. (19,20) However, differential gene glucocorticoid receptor (GCR) genes in circulating monocytes potentially contributed to bone metabolism. (16) The present study aims to identify significantly differentially expressed genes from 168 selected cytokine, chemokine, and osteoclastogenesis-related genes in circulating monocytes between the high and low BMD groups in Chinese Han females and validate the significant expression in Caucasian women. We also performed single-nucleotide polymorphism (SNP) association analysis with BMD to find further evidence of the identified genes at the DNA level.

Chinese subjects
The study was approved by the Research Administration Department of Hunan Normal University. Eight hundred and seventy-eight females who were Chinese Hans were recruited from Changsha City. All subjects signed informed-consent documents before entering the project. Healthy female subjects aged of 20 to 45 years were included because BMD reaches its peak and is most stable during this period. For each subject, we collected information on age, sex, medical history, family history, menstrual history, smoking history, physical activity, alcohol use, tea and coffee consumption, diet habits, etc. Female subjects must have regular menses to avoid the effects of menopause on BMD. Subjects with chronic diseases and conditions that potentially may affect bone mass have been excluded from the study. These diseases/conditions included chronic disorders involving vital organs (e.g., heart, lung, liver, kidney, and brain), serious metabolic diseases (e.g., diabetes, hypo-and hyperparathyroidism, and hyperthyroidism, etc.), skeletal diseases (e.g., Paget's disease, osteogenesis imperfecta, and rheumatoid arthritis, etc.), chronic use of drugs affecting bone metabolism (e.g., corticosteroid therapy and anticonvulsant drugs), and malnutrition conditions (e.g., chronic diarrhea, chronic ulcerative colitis, etc.). From the 100 top and 100 bottom hip BMD subjects we recruited all who consented to enter our potential future projects, including 14 high hip BMD (mean AE SD ¼ 1.03 AE 0.05 g/cm 2 ) subjects and 12 low hip BMD (mean AE SD ¼ 0.7 AE 0.06 g/cm 2 ) subjects (Table 1). Thirty milliliters of peripheral blood were drawn for each selected subject. BMD measurement BMD (g/cm 2 ) at the lumbar spine (L1-4, anteroposterior view) and total hip (femoral neck, trochanter, and intertrochanter region) was measured by a Hologic 4500-W dual-energy X-ray absorptiometry (DXA) (Hologic Corp., Waltham, MA, USA). The DXA scanner was calibrated daily, and long-term precision was monitored with external spine and hip phantoms. The coefficient of variation (CV) of measured BMD values was 0.80% at the hip.

Monocyte isolation
A monocyte negative isolation kit (Dynal Biotech, Inc., Lake Success, NY, USA) was used to isolate circulating monocytes from 30 mL of whole blood following the procedures recommended by the manufacturer. The kit contains a mixture of antibodies for CD2, CD7, CD16, CD19, CD56, and CD235a to deplete T cells, B

Statistical analysis
Transciptome-wide expression profiling involves a large number of genes and thus incorporates tremendous multiple tests. Some genes with suggestive significance may be excluded after the multiple-testing correction. If some genes are tightly related in a functionally relevant pathway, the P value for any single gene may not be significant after the multiple-testing correction. However, a focused expression screen on potential functional relevant genes largely can reduce the multiple tests and may increase the statistical power. Considering the multicomparison problem and the function of monocytes, 168 candidate genes were selected for statistical analyses. The 168 genes are all available cytokines, chemokines, osteoclast-related factors, and their receptors selected for our focused expression analyses of the Affymetrix HG133 plus 2.0 gene data set (see Appendix  table). Microarray Suite 5.0 (MAS 5.0, Affymetrix) software was used to generate the array raw data files (CEL files). Then the probe-level data in CEL files were converted into expression measures and normalized by the robust multiarray average algorithm (RMA, www.bioconductor.org). (21) The differential expression analysis between low and high BMD samples was conducted by a nonparameter Wilcoxon signed-rank test. A Benjamini and Hochberg (BH) stepwise procedure was used for multiple-comparison adjustment, (22) and an adjusted P .05 was used as the significant criterion. Fisher's exact test was used in the canonical pathway analysis by Ingenuity Pathways Analysis (IPA) software (Ingenuity Systems, www.ingenuity.com) to test the association between genes within a canonical pathway and BMD variation. According to the similarity of gene expression, the differentially expressed genes were further analyzed for twodimensional hierarchical clustering at both gene and sample levels. (23) Array replication in Caucasians In this independent microarray study, we recruited 20 premenopausal Caucasian women, 10 with high BMD (spine or hip Z-score greater than þ0.84) and 10 with low BMD (spine or hip Z-score less than À0.84; see Table 1) from the vicinity of Creighton University in Omaha, Nebraska, USA, for differential expression analyses in their circulating monocytes. Although the weights in the high BMD group were higher than in the low BMD group, no significant correlation of BMD with weight was detected in the low or high BMD group. This study was approved by the Institutional Review Board, and all the subjects signed informed-consent documents before entering the project. The inclusion and exclusion criteria were the same as in the Chinese population, but the age criterion was limited to the narrow range of 39 to 45 years, within the peak BMD range of 20 to 45 years. For the Caucasian samples, we used the Affymetrix HG-133A instead of the HG-U133 plus 2.0 that we used for the Chinese samples (the HG-133A chip contains fewer genes than the HG-U133 plus 2.0), but all the other experimental procedures and statistical analyses were the same as we described for Chinese samples.

Validation by qRT-PCR in Caucasians
We used two-step qRT-PCR to confirm differentially expressed genes. Reverse-transcription reactions were performed in a 50 mL reaction volume containing 5 mL 10Â PCR Buffer II, 11 mL 25 mM MgCl 2 , 10 mL dNTPs, 1.25 mL MULV reverse transcriptase, 1.0 mL RNase inhibitor, 2.5 mL Oligo d(T), 0.5 mg total RNA, and water to 50 mL. All these reagents were supplied by Applied Biosystems (Foster City, CA, USA). Reaction conditions were as follows: 10 minutes at 258C, 30 minutes at 488C, and 5 minutes at 958C. Realtime quantitative PCR was performed in a 25 mL reaction volume using standard protocols on an Applied Biosystems 7900HT. Briefly, 2.5 mL of cDNA was mixed with 12.5 mL of TaqMan universal PCR master mix (2Â), 1.25 mL of TaqMan gene express assay mix (contains forward and reverse primers and labeled probe), 1.25 mL of human GAPDH probe (20Â), and 7.5 mL of water. The thermocycling conditions were as follows: 2 minutes at 508C, 10 minutes at 958C, and 40 cycles of 15 seconds at 95 8C plus 1 minute at 608C. Based on the relative gene expression 2 ÀDDCt , (16) we performed Student's t test to confirm the differential expression genes. All reactions were run in triplicates for each gene.

Subjects and phenotype
For the association study, 1000 unrelated Caucasian subjects were identified from our established and expanding genetic repertoire, currently containing more than 6000 subjects. All subjects were U.S. Caucasians of European origin. The inclusion and exclusion criteria were the same as in the Chinese population for expression study, but males and postmenopausal women were included. The basic characteristics of all subjects are listed in Table 1. BMD values at spine and hip were measured using the Hologic 4500A DXA (Hologic, Inc., Bedford, MA, USA). The coefficient of variation of the DXA measurement was approximately 1.98% for spine BMD and 1.87% for hip BMD.

Genotyping and statistical analysis
Genomic DNA was extracted from whole human blood using a commercial isolation kit (Gentra Systems, Minneapolis, MN, USA). Genotyping with the Affymetrix Mapping 250K Nsp and 250K Sty arrays was performed. Fluorescence intensities were quantified using an Affymetrix Array Scanner 30007G. Data management and analyses were performed using the Affymetrix GeneChip Operating System. The final average Bayesian Robust Linear Model with Mahalanobis (BRLMM) call rate across the entire sample reached the high level of 99.14%. We tested the association of significant genes identified in the expression studies with BMD in the 1000 Caucasian subjects. Parameters such as age, age 2 , sex, age/age 2 -by-sex interaction, height, and weight were tested for their association with BMD at spine and hip. The significant ( P .05) terms then were included as covariates to adjust the raw BMD values for subsequent analyses. Statistical analyses were performed using genotype and haplotype association software implemented in PLINK-1.03 (http://pngu.mgh.harvard.edu/purcell/plink/). (24) Linkage disequilibrium (LD) patterns were analyzed and plotted with the correlation coefficient between pairs of loci based on the 1000 unrelated Caucasians using the the Haploview program (www.broad.mit.edu/mpg/haploview/), (25) which describes more or less combinations of alleles. The haplotype block was used to show chromosome regions with high LD and low haplotype diversity in haplotpye association studies. MAPPER was used for searching transcript factor binding sites in the JASPAR database (http://mapper.chip.org/). (26)

Results
The basic characteristics of the study subjects are shown in Table 1. Although hip BMD is the major study phenotype for the expression analyses in Chinese, spine BMD also was significantly different between the high and low BMD groups ( P ¼ 1.49 Â 10 À6 ). There were no significant differences in age and height traits between the high and low BMD groups for both Chinese and Caucasians for expression analyses. However, the weight and body mass index (BMI, kg/m 2 ) in the low BMD group were significantly lower than in the high BMD group in Caucasians. We performed a general linear regression analyses for STAT1 expression values and BMD status and incorporated weight and height as covariates in Caucasians. However, weight and height are not significant as covariates for STAT1 expression analysis in the high and low BMD groups ( P ¼ .2529 and .2045, respectively). This implied that the weight and height in current study were not confounding factors for BMD in our expression analyses.
All the nominally significant genes ( P < .05) of differential expression with BMD in Chinese are summarized in Table 2. We submitted our gene expression profiling to the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo/), and the access number was GSE7158. After Benjamini and Hochberg correction for multiple comparisons, the differential expressions of signal transducer and activator of transcription 1 (STAT1) (adjusted P ¼ .02248) and guanylate binding protein 1 (GBP1) (adjusted P ¼ .03372) genes were still significant between the low and the high BMD groups. Expression fold changes of the significant genes were not large in the present study. The main reason might be that the female subjects in this study all were 20 to 45 years of age and with regular menses, in contrast with including both premenopausal and postmenopausal women in our previous study. (16) Actually, a 1.5-fold change has been shown to be significant in other differential gene expression studies. (27)(28)(29)(30) In addition, BMD is a complex trait, and genes regulating its variation are not expected to have large differential expressions. In Fig. 2, a two-dimensional hierarchical dendrogram (based on the nominally significant genes listed in Table 2) shows the results of the hierarchical clustering analyses. Low BMD subjects were mainly clustered to the bottom of the figure.
In the interleukin system, interleukin 1 receptor antagonist (IL1RN) and interleukin 15 (IL15) genes were upregulated in the low BMD group compared with the high BMD group. However, IL6 gene expression was not detected.
Interestingly, the differential expression of the STAT1 gene in circulating monocytes was replicated in our ongoing comparison microarray study between the low and high BMD premenopausal Caucasian women. Consistently, upregulation of the STAT1 gene in the low BMD group also was significant after Benjamini and Hochberg correction for multiple testing ( P ¼ .0028, adjusted P ¼ .048) ( Table 3). qRT-PCR confirmed the significant differential expression of the STAT1 gene in Caucasians ( P ¼ .0046) (see Table 3). For the GBP1 gene, we did not find any significant results in both array and qRT-PCR analyses in Caucasians.
In SNP genotyping analysis using an additive model, two SNPs, rs10199181 ( P ¼ .0028) and rs2030171 ( P ¼ .0264), in the STAT1 gene were associated with spine BMD (Table 4). It was obvious that subjects with the T allele of rs10199181 possessed high spine BMD (Fig. 3). Figure 4 shows the correlation coefficient between pairs of SNPs of the STAT1 gene and reconstructed haplotype blocks. Interestingly, rs10199181 and rs2030171 were located in the same block, and a haplotype composed of SNPs rs16833157-rs2030171-rs10199181 (''G-G-A'') in the block also was demonstrated to be significantly associated with spine BMD ( P ¼ .0029). However, no significant association was detected for 14 SNPs in GBP1 with BMD in Caucasians.

Discussion
In this study we investigated expression of 168 genes related to cytokines, chemokines, osteoclast formation factors, and corresponding receptors in monocytes from Chinese Han women with extremely discordant BMD. Thirteen genes were found to be differentially expressed. A very interesting phenomenon was that among the 13 genes, the STAT1, IFI44L, CXCL10, IFI44, GPB1, and GPB2 genes were expressed in higher levels in the low BMD group than in the high BMD group, which  Note: Hybridization intensity and ''present'' status were given based on MAS5 algorithm. L-BMD intensity means hybridization intensity from the low BMD group; H-BMD intensity means hybridization intensity from the high BMD group. Fold L/H means the ratio of hybridization intensity from the low to that from the high BMD group. Raw P value means P value before multiple testing corrections. Adjusted P value represents P value adjusted with Benjamini/Hochberg method, and asterisks Ã means significant after Benjamini-Hochberg correction considering 281 probes of selected genes.
is very similar to the IFN-induced gene expression pattern (IFN pathway). For instance, immature peripheral blood mononuclear phagocytes stimulated by the type I IFN isoform increased the expression of 44 genes, including STAT1, IFI44L, CXCL10, IFI44, GPB1, and GPB2. (31) Microarray analysis of cells infected with short hairpin RNA vectors pAB319 and pAB322 showed enhanced expression of many IFN target genes, such as STAT1, IFI44L, CXCL10, IFI44, and GPB1. (32) The increased expression of IFN pathway genes also was detected in blood mononuclear cells from patients with systemic lupus erythematosus and juvenile dermatomyositis (33) and in the human fibrosarcoma cell line. (34) The IFN pathway may regulate bone resorption in two ways. First, interferon-g (IFNG) blocks RANKL-induced osteoclast differentiation. (9,35) Second, the IFN pathway in circulating monocytes may stimulate the secretion of cytokines IL-1, IL-6, and TNF to increase bone resorption. (36)(37)(38) In this study, the upregulation of STAT1- Note: Expression value of Affymetrix Microarray was the hybridization intensity based on the MAS5 algorithm; expression value of qRT-PCR was relative quantity based on 2 ÀDDCt ; fold L/H means the ratio of gene expression from the low to that from the high BMD group; Ã indicates the BH-adjusted P values for multiple testing.   mediated IFN pathway genes in the low BMD group suggested the important effect of STAT1 on bone resorption in vivo in humans. For the 13 differentially expressed genes, however, only the STAT1 and GBP1 genes remained significant after correction for multiple testing. In our previous genome-wide array study on the same Chinese samples, we also found significant differential expression of the STAT1 and GBP1 genes in the array data analyses after correcting for multiple testing. (39) However, further qRT-PCR only confirmed the significance of the GBP1 gene but not the STAT1 gene. (39) Thus we tried to replicate significance of the two genes in our Caucasian expression study on circulating monocytes from 20 premenopausal Caucasian women (10 with low BMD and 10 with high BMD) and SNP association study on 1000 unrelated Caucasian subjects. We did not find the significance of the GPBP1 gene in either replication study. Interestingly, however, the significance of the STAT1 gene was found in both replication studies. In particular, significant upregulation of the STAT1 gene was found in both array and qRT-PCR experiments in the Caucasian expression study. The STAT1 gene was not differentially expressed in B cells isolated from peripheral blood between the high and low BMD subjects (data not shown) who were the same Caucasians for our current monocyte study. Therefore, it is likely that the alterations in STAT1 expression only in monocytes, but not in other cells, are responsible for variations in bone mass in humans.
In the IFN signaling pathway, STAT1 is a critical mediator gene. (9,40) In the above-mentioned IFN pathway for regulating bone resorption, STAT1 mediates the effects of IFNG on both inhibition of RANKL-induced osteoclast differentiation (9,35,41) and secretion of IL-1, IL-6, and TNF. (36)(37)(38) In addition, in dexamethasone-treated peripheral blood mononuclear cell (PBMC) cultures, the inhibited IFNG expression suppressed expression of the STAT1 gene. (42) Furthermore, in lupus nephritis patients, basal expression of STAT1 was significantly higher in monocytes. and stimulation of the monocyte cultures with IFNG resulted in phosphorylation of STAT1. (43) In mice, the STAT1 gene plays an important role in bone metabolism in osteoblasts. (44) Recently, STAT1 was reported to be upregulated in femur tissue in osteoporotic mice, (45) and this supports our finding of the high expression of STAT1 in monocytes in human low BMD groups.
Interestingly, our published linkage study of BMD in 4126 human subjects also identified suggestive univariate and significant epistatic linkage signals at 2q32, which harbors the STAT1 gene. (46) Furthermore, our group recently found significant linkage evidence on 2q32 with spine BMD using bivariate linkage analysis. (47) Our current SNP association study also replicated the significance of the STAT1 gene for spine BMD in Caucasian samples. No SNP in the STAT1 gene was associated with hip BMD at the SNP level, perhaps owing to different genetic determinants for spine BMD and hip BMD traits because many studies have shown different heritability and genetic loci underlying the two traits. (47,48) Hence our results tend to reveal the significance of STAT1 on spine BMD. Subjects with the T allele of SNP rs10199181 in the STAT1 gene tended to have a higher spine BMD than those with other alleles (see Fig. 3). According to the transcript factor Jaspar database, the T allele of rs10199181 is likely to bind transcript factor E4BP4, which might be induced by parathyroid hormone (PTH), a well-know hormone for bone growth, in osteoblasts. (49,50) The inducible effect of E4BP4 suggests a negative regulation by glucocorticoids that might decrease BMD. (51) Thus it implies that the T allele of rs10199181 in the STAT1 gene may be involved in bone growth metabolism. Based on the present results and previous knowledge, we developed a novel mechanism for osteoclastogenesis (Fig. 5). In peripheral blood, IFN mediated by STAT1 may stimulate circulating monocytes to produce cytokines such as IL-1, TNF, CXCL10, and IL-15 that increase the bone resorption function of osteoclasts. Another pathway in the bone microenvironment also may be triggered by upregulated STAT1 and IFN, which may inhibit osteoblatogenesis. In this study, no expression of the RANK and TRAF6 genes in circulating monocytes was detected, suggesting that osteoclast formation was completely inhibited in circulating monocytes. The fact that osteoclast differentiation was not initiated in peripheral circulating monocytes in both the high and low BMD groups is possibly because osteoclast formation from monocytes may occur only in a special microenvironment. (52) Our current research studied and found the importance of only the STAT1 gene in the function of monocytes or osteoclasts on bone metabolism, which, however, did not address the reported importance of the STAT1 gene in osteoblasts. (9,44) In summary, our results support the fact that the STAT1 gene in circulating monocytes plays important roles in bone metabolism and also suggests that gene expression of the STAT1-mediated IFN pathway may be important for osteoporosis.

Disclosures
The authors state that they have no conflicts of interest. Note: The data were sorted by adjusted P value. Hybridized intensity and ''present'' status were given based on the MAS5 algorithm. L-BMD intensity means average hybridized intensity in the low BMD group; H-BMD intensity means average hybridized intensity in the high BMD group; Fold L/H means the ration of L-BMD intensity to H-BMD intensity. Raw P value represents P value before multiple testing correction. Adjusted P value represents P value adjusted with Benjamini/Hochberg method; asterisks indicate significant results after Benjamini-Hochberg correction considering 281 probes of selected genes.