GeneChip Expression Profiling Reveals the Alterations of Energy Metabolism Related Genes in Osteocytes under Large Gradient High Magnetic Fields

The diamagnetic levitation as a novel ground-based model for simulating a reduced gravity environment has recently been applied in life science research. In this study a specially designed superconducting magnet with a large gradient high magnetic field (LG-HMF), which can provide three apparent gravity levels (μ-g, 1-g, and 2-g), was used to simulate a space-like gravity environment. Osteocyte, as the most important mechanosensor in bone, takes a pivotal position in mediating the mechano-induced bone remodeling. In this study, the effects of LG-HMF on gene expression profiling of osteocyte-like cell line MLO-Y4 were investigated by Affymetrix DNA microarray. LG-HMF affected osteocyte gene expression profiling. Differentially expressed genes (DEGs) and data mining were further analyzed by using bioinfomatic tools, such as DAVID, iReport. 12 energy metabolism related genes (PFKL, AK4, ALDOC, COX7A1, STC1, ADM, CA9, CA12, P4HA1, APLN, GPR35 and GPR84) were further confirmed by real-time PCR. An integrated gene interaction network of 12 DEGs was constructed. Bio-data mining showed that genes involved in glucose metabolic process and apoptosis changed notablly. Our results demostrated that LG-HMF affected the expression of energy metabolism related genes in osteocyte. The identification of sensitive genes to special environments may provide some potential targets for preventing and treating bone loss or osteoporosis.


Introduction
High magnetic fields (HMFs) are one of the most powerful tools for studying the properties of materials because they couple directly to the electronic charge and magnetic moments of the protons, neutrons, and electrons [1]. Recent technologic innovations have led to the generation of man-made static magnetic fields up to 10 Tesla (T). HMFs produced by a superconducting magnet have been widely used in research and medical applications. HMFs (>10 T) affected MLO-Y4 cells exposed to LG-HMF. Down-regulated genes were much more than up-regulated genes in set 1 (Fig. 1A), while it presented a reversed feature in set 2 (Fig. 1B). Set 3 contained the least number of genes with lower fold change than the other three sets (Fig. 1C). In set 4, there was a significant increase in down-expressed genes compared to other sets (Fig. 1D).
The DEGs were refiltered by setting the cutoff limitations. Besides the initial fold change, cutoff of 2-and 1.5-fold were applied. The number of DEGs (FC > 1.5 and 2) in 4 sets were shown in Table 1. The number of genes decreased with the increase of the fold change. There are a few more down-regulated genes in the set1, set 3 and set 4 relative to set 2, but in set 2 there are much more up-regulated genes compared with other sets (Table 1). In set 1 and set 3, all DGEs (FC > 2) were down-regulated genes. In set 2, the number of up-regulated DGEs was more than that in down-regulated DGEs. Figure 1. Volcano plots of differentially expressed genes in MLO-Y4 cells exposed to LG-HMF. Volcano plots displays unstandardized signal against noise-adjusted/standardized signal. The x-axis represents the fold change cutoff, while y-axis shows the negative logarithmic of P value. A: μ-g v.s. control (set 1), B: 2-g v.s. control (set 2), C: 1-g v.s. control (set 3), D: μ-g v.s. 2-g (set 4).
The relationship of DEGs (FC > 2) among set1, set3 and set 4 was further analyzed (Fig. 2). Most of DEGs in set 1 and set 3 are same except 3 genes (CRCT1, ALDOC and Higd1a). There were one gene (CA12) in set 3 and two genes (MGARP and CRCT1) in set 1 different from set 4.

Molecular function and cellular location of DEGs in different sets
In order to further identify interesting new target genes, we used iReport data analysis system to analyze the molecular function and cellular location of DEGs. Molecular function of DEGs (with a 1.5-fold change) in 4 sets was analyzed by iReport data analysis system. More than 20% DEGs in 4 sets belonged to enzyme ( Table 2). In set 4, there were 9 and 11 DEGs pertained to transciption regulator and transporter, respectively. In the set 2, the molecular function of 6 DEGs was related to cytokines and transporter, respectively (Table 2). Moreover, several DEGs related to G-protein coupled receptor were presented in all the four sets.
The percentages of DEGs loacted in cytoplasm ranked the first in all the four sets, and these genes take up 48.9%, 31.2%, 55.3% and 36.3% in set 1, set 2, set 3 and set 4, respectively ( Table 2). More than 10% DEGs distributed in plasma membrane in 4 sets. Genes located in extra cellular space in set 2 were more than those in any other three sets (Table 3).

Functional annotation clustering of DEGs in different sets
To reduce the burden of associating similar redundant terms and make the biological interpretation more focused, we utilized DAVID funtional clustering to measure relationships among the annotation terms based on the degrees of their co-association genes. We selected the terms with the smallest P value and a enrichment score more than 2. Totally, 5 subsets of genes were culsterd based on GO in the three sets ( Table 4). The subsets belonged to set 1 and set 4 presented an evident association with the glucose metabolic process (Table 4). These results were enhanced by the subsets in SP-PIR-Keywords and the two groups of genes enriched by KEGG-Pathways (Table 4). Genes of set 2 were clustered into two GO categories, one subset was marked by oxidoreductase activity, and another showed notable location in extracellular region (Table 4). In set 3, the clear clustering genes were not found.
Ingenuity iReport can be used to filter, group, and visualize genes by function, biological process, role in pathway or disease. In order to further analyze DEGs associated to biological process, we chose Ingenuity iReport to filter genes statistically significant associated to biological processes. 20 biological processes with the minimum P value and with the most number of DEGs were showed in Fig. 3 Glycolysis of cells ranked first in set1 and set 4 because of the minimum P value (Fig 3A). The apoptosis process was markedly involved in set1 and set 4 (Fig. Expression Profiling of LG-HMF on Osteocytes 3B). In set 2, multiple biological processes, such as cell viability, cell movement were clustered ( Fig. 3A and B). In set 3, the clustering genes are mainly related to disease.

Verification of DEGs sensetive to distinct apparent gravity levels by qPCR
In order to verify the corrections of microarray data, we selected 12 DEGs from microarray data and real-time PCR was used to verify the effects of LG-HMF on these gene expression at mRNA levels. The selected DEGs could be classified into 4 subgroups according to functional clustering: enzyme related genes (CA9, CA12 and P4H A1), G-protein coupled receptors (GPR35 and GPR84), peptide hormone (STC1, ADM and APLN) and genes related to energy metabolism (PFKL, AK4, ALDOC and COX7a1). After being normalized by internal control genes, the relative gene expression levels in experimental groups were obtained comparing with those of control groups. And then, the differences in gene expression between m-g v.s. control, 1-g v.s. control, and m-g v.s. 2-g were analyzed (Table 5 and Fig. 4). Fold changes of 12 selected genes in m-g v.s. control, m-g v.s. 2-g, and 1-g v.s. control by qPCR and microarray  Expression Profiling of LG-HMF on Osteocytes analysis were showed in Table 5. The tendency of microarray analysis and PCR results were similar. Except for GPR84, PH4HA1 and APLN, the other 9 genes (GPR35, PDK1, AK4, ADM, COX7, STC1, ALDOC, CA9 and CA12) expression significantly decreased in m-g v.s. control, 1-g v.s. control and m-g v.s. 2-g (Fig. 4). The expressions of GPR 84 obviously decreased in Table 4. Functional annotation cluster of DEGs in four sets (FC>1.5). Expression Profiling of LG-HMF on Osteocytes m-g v.s. control and 1-g v.s. control but increased in m-g v.s. 2-g. APLN expression decreased in m-g v.s. 2-g and 1-g v.s. control but slightly increased in m-g v.s. control (Table 5 and Fig. 4).

Bio-data mining of the verified genes
Bio-data mining was performed on the basis of Ingenuity Knowledge Base. Interactions between DEGs were analyzed. We focused on the 12 verified genes and mapped them by the interaction network (Fig. 5). Notably, genes presented much more complex interactive relation in set 2 and set 4 than those in set 1 (Fig. 5A, 5B). However, in set 3 clear interactive relations among 12 verified genes were not found. Interestingly, in view of m-g v.s. 2-g (Fig. 5C), we got several hub genes linking the genes of interest (the red fond genes in Fig. 5). EPAS1, ADM and AGT (Fig. 5), as well as the TNF (Fig. 5B) played as joints between genes. Moreover, according to literature mining of Ingenuity, we mined those genes biological information deeply. The disease processes and pathways associated with DEGs were filtered. Diseases of bone metabolisms were presented with the corresponding set of genes (Table 6). Particularly, CA9 and CA12 play role in osteoporosis (Table 6).

Discussion
As a novel technology, the diamagnetic levitation technique has caused more and more attention and has been applied in many fields, such as material sciences, biology, and chemistry. In this study, the effects of diamagnetic levitation on gene expression profiling in osteocytes have been investigated for the first time. Our previous results showed that the cellular morphology . Biological processes associated to differentially expressed genes in four sets. Biological processes were mapped by Ingenuity knowledgebase. Line A shows the the most statistically significant biological processes of set 1, set 2 set 3 and set 4. line B shows the biological processes (P < 0.05) involving most differentially expressed genes. Fisher's exact test was used to calculate the P value. Set 1: μ-g v.s. control; Set 2: 2-g v.s. control, Set 3: 1-g v.s. control; Set 4: μ-g v.s. 2-g.
doi:10.1371/journal.pone.0116359.g003 Table 5. Fold change of DEGs tested by qPCR & microarray in μ-g vs. control, μ-g vs. 2-g, and 1-g vs. control by qPCR and microarray assays.  . Microarray results were verified using real-time PCR for genes in response to LG-HMF. Total RNA was extracted and qPCR assay was used to further identify for 12 selected genes. The method of relative quantification was used to estimate the relative expression changes of selected gene expression in MLO-Y4 cells exposed to LG-HMF. The changes in selected gene expression, normalized to 18S under LG-HMF were calculated. The difference between μ-g v.s. control, μ-g v.s. 2-g, and 1-g v.s. control was statistically analyzed by one-way ANOVA. μ -g, 1-g, 2 -g v.s control group: ***P< 0.001; **P< 0.01; *P< 0.05. μ -g v.s 2 -g group: ### P < 0.001; ## P < 0.01.  Expression Profiling of LG-HMF on Osteocytes and cytoskeleton of osteocytes changed dramatically after cultured in LG-HMF for 2 days, and the expression of genes presented a quite different scene [23]. Based on these results, we further investigate the effects of LG-HMF on the gene expression profiling in osteocytes. The novel and most significant finding is that exposure of osteocytes to LG-HMF (m-g, 1-g, and 2-g) distinguishes some genes that are sensitive to low gravity, magnet field, and the combined environment. The results are helpful to improve our understandings of how cells sense altered gravity and the mechanisms bone loss induced by weightlessness at a cellular level. Since a high magnetic field coexists with different gravity levels at all time, four groups were designed in this study, namely 1 g group (normal gravity, 16 T), control group (normal gravity, geomagnetic field), 2 g group (2-fold gravity, 12 T), and m-g group (hypogravity, 12 T). In order to relatively distinguish the effects of magnetic fields and different apparent gravities, we named four sets as set 1(m-g v.s. control), set 2 (2-g v.s. control), set 3 (1-g v.s. control) and set 4 (m-g v.s. 2-g). Set 1, set 2, set3 and set 4, respectively, showed the effects of diamagnetic levitation, hypergravity, magnetic field and hypogravity.
By using iReport and DAVID analysis, DEGs in 4 sets were obtained. In set1, set 2 and set 4, the number of down regulated DEGs was much more than that of up regulated DEGs but it was converse in set 3. The results indicate that hypergravity increases most of DEGs expression while hypogravity or magnetic field mainly decreases their expression. Moreover, 10 of DEGs (FC>2) were sensitive to the combined environment (Fig. 2). iReport data analysis also showed that more than 20% DEGs in 4 sets belonged to enzyme, moreover, more than 30% DEGs located in cytoplasm in all the four sets. Manchester reported that space flight obviously affected eleven enzymes in indibidual fibers of soleus and tibilis anterior muscules [41]. The findings suggest that enzyme related genes may be very sensitive to extreme environment.
In order to decrease the similar redundant terms and make the biological interpretation more focused, we used DAVID funtional clustering to detect relationships among the annotation terms. The results showed that DEGs associated with glucose metabolic process or glycolysis were strikingly presented in set 1, set 2 and set 4 groups ( Table 4). The results suggest that abnormal gravity may affect osteocytes metabolism. Glycolysis is the universal pathway used by all the organisms to extract energy from glucose. Ramirez et al., reported that when mice were subjected to hind limb suspension, the glycolysis was inhibited, while gluconeogenesis was up regulated in liver [42]. The results given by iReport also confirmed the effects of abnormal gravity on glucose metabolism in osteocytes (Fig. 3). Genes related to the apoptosis, necrosis and cell movement processes were also sorted in set 1, set 2 and set 4 by iReport, but genes in set 3 did not present the similar clustering (Fig. 3). These results indicated that abnormal gravity affected osteocyte functions, such as apoptosis, necrosis and cell movement processes but the high magnetic field did not involve in these processed. Furthermore, all these results suggest that osteocytes might respond the mechanical changes through one or more of these processes. Totally 12 DEGs were concerned and verified because of their significant changes and their involvements in biological processes. Both of PCR and microarray analysis showed the expression of 12 DEGs (CA9, CA12, P4HA1, ADM, STC1, APLN, GPR35, PFKL, AK4, ALDOC, COX7A1) was significantly changed in m-g v.s. control, 1-g v.s. control and m-g v.s. 2-g. It suggests that these DEGs are sensitive to both altered gravity and high magnetic field. The carbonic anhydrases (CA) belong to a family of enzymes that catalyze the rapid interconversion of carbon dioxide and water to bicarbonate and protons. CA9 and CA12 are two members of carbonic anhydrase family [43][44][45]. It has been known well that another CA isoform, CA2, takes an active part in the bone resorption of osteoclasts by regulating the osteoclastic PH. Moreover, CA9 and CA12 had been supposed to act as the Synergy factor of CA2 [43]. Alteration of CA2 in both flight and suspended animals after readaptation to Earth gravity [46].These evidences taken together prompted that CA9 and CA12 in osteocytes might do something crucial in mechano-induced bone remodeling. Prolyl 4-hydroxylase subunit alpha-1 (P4HA1), like CA9 and CA12 belongs to enzymes, which plays a key role in collagen synthesis. Decrease in P4HA1 expression demonstrates that LG-HMF may affect collagen synthesis in osteocytes and ultimately impact bone formation.
Adrenomedullin (ADM) is a peptide hormone that in humans is encoded by the ADM gene. ADM was reported to inhibit the osteoclastogenesis [47], while stimulate osteoblast growth and proliferation [48]. Up-regulation of the ADM system occurred early and took part in the adaptative changes occurring during simulated microgravity conditions [49]. The decrease of ADM expression in osteocytes under LG-HMF condition suggests ADM may be involved in response of osteocytes to the extreme environment. Stanniocalcin-1 (STC1) is a glycoprotein hormone involved in calcium/phosphate (Pi) homeostasis. Filvaroff, et al. has reported that STC-1 can affect calcium homeostasis, bone and muscle mass and structure, and angiogenesis through effects on osteoblasts, osteoclasts, myoblasts/myocytes, and endothelial cells [50].The reduction of STC1 expression in osteocytes under LG-HMF conditions suggests STC1 may be involved in abnormal bone remodeling process. Apelin (APLN) is a peptide that is encoded by the APLN gene, and APLN receptor expression is observed at the surface of osteoblasts, the cell progenitors involved in bone formation [51], meanwhile, lack in APLN increased bone mass in mice osteoblast [52]. Our results showed that APLN expression significantly decreased in osteocytes in m-g v.s. 2-g, which indicates that APLN may be sensitive to altered gravity.
G-protein coupled receptors, GPR35 and GPR84 expression dramatically changed in osteocytes under LG-HMF. GPR35 functions as a receptor for the kynurenine pathway intermediate kynurenic acid, which elicits calcium mobilization and inositol phosphate production [53]. GPR84 is highly expressed in the bone marrow, and in splenic T cells and B cells [54] and may link fatty acid metabolism to immunologic regulation [55]. Thesed findings indicate that Gprotein coupled receptors may directly or indirectly participate in regulating bone metabolism.
Gluconeogenesis is required for the living organisms to grow at the expense of carbon as energy source other than carbohydrates and capable of synthesizing glucose from simple starting materials [56]. In this study, DAVID funtional clustering showed that several DEGs associated with glucose metabolic process or glycolysis were obviously presented under LG-HMF. PCR and microarray results showed that the expression of energy metabolism related genes, PFKL (6-phosphofructokinase, liver type), AK4 (adenylate kinase 4), ALDOC (aldolase C, fructosebisphosphate) and COX7A1 (cytochrome c oxidase subunit VIIa), dramatically decreased. These results demonstrate that LG-HMF affect metabolic enzymes in osteocytes. The changes in metabolic enzymes in mice liver or muscle fibers were presented under space flight or simulated microgravity conditions [41,42].
In order to further investigate the function of DEGs, we analyzed the interaction among DEGs by the interaction network. The results showed that there were several hub genes linking the genes of interest in m-g v.s. 2-g, such as ADM, P4HA1. And the gene interaction network also showed that some genes, including EPAS1 (endothelial PAS domain protein 1), TNF (tumor necrosis factor), AGT (angiotensinogen) were involved in regulating DEGs. The gene interaction network provides useful clues for further study in future. Disease processes in which DEGs participated were selected by iReport system. Particularly, CA9 and CA12 play role in osteoporosis. These results further imply the importance of carbonic anhydrases in bone disease.
In summary, the present study used DNA microarray analysis to provide a new and comprehensive cognition to the effects of LG-HMF on gene expression profiles in osteocyte-like cells, and has selected 12 genes (CA9, CA12, P4HA1, ADM, STC1, APLN, GPR35, PFKL, AK4, ALDOC, COX7A1) that may be sensitive to altered gravity or magnetic field. The study shows that LG-HMF affects the expression of several kinds of genes related to enzyme, peptide hormone, G-protein coupled receptors and glucose metabolic process. The identification of mechanosensitive genes will help us to understand the mechanism of bone loss to open a new route for the therapeutic control of bone mass and provide new potential countermeasures.

Materials and Methods
Cell culture MLO-Y4 osteocyte-like cell gifted by Dr. Lynda Bonewald [57] were cultured in α-Modified Eagle's Medium (α-MEM, Gibco, Paisley, UK) containing 5% fetal bovine serum, 5% calf serum (Gibco, Paisley, UK), 1% benzylpenicillin and 1% streptomycin. MLO-Y4 cells grew on culture flask coated with collagen (rat tail collagen type 1, 0.15 mg/ml, BD, USA). Once reaching 80%-85% confluence, cells were digested by trypsin containing 0.03% EDTA, and seeded onto 96-well plates (9102; Corning Costar,Corning, NY, USA) with a density of 30000 per well, then ten wells were placed into a 35-mm tissue-culture plate (Nunc, Inc., Roskilde, Denmark). And then the plate was delivered to the appropriate (m-g, 1-g, and 2-g) in the bore of the superconducting magnet by the object holder to continuously culture for 48 hours at at 37°C with 5% CO 2 . The control group was incubated at 37°C with 5% CO 2 in normal condition.

Superconducting Magnet with Large Gradient High Magnetic Field
Superconducting magnet with LG-HMF was manufactured by Japan Superconductor Technology, Inc. (JASTEC) according to the specific specifications proposed by authors. Specifications of the superconducting magnet were shown in Table 7. The height of the superconducting magnet is 195 centimeters and a F51mm×450mm cylindrical cavity can be used for experiment. The superconducting magnet can generate a magnetic force field (B ·dB/dz) of −1370, 0, and 1370 T 2 /m in a 51-mm diameter room temperature (RT) bore, corresponding to three apparent body force levels (m-g, 1-g, and 2-g) and three magnetic induction intensities (12,16, and 12 T), respectively. The experimental platform for diamagnetic levitation of biological systems has been further developed based on the superconducting magnet by the authors [39,40]. The experimental platform mainly contains four sections: superconducting magnet (JASTEC, Japan) providing large gradient high magnetic gravity environments, temperature control system, object stage, gas control system and observing system. The monitoring device was integrated into the object stage to measure the gravity, temperature, and displacement. The temperature control system includes a water-bath pump and a channel system, and the temperature range for the control system was 37 AE 0.5°C.

Gene expression profiling by DNA microarray
Total RNA was isolated from MLO-Y4 cells exposed to LG-HMF and controls for 48 h using Trizol method as recommended by the manufacturer's protocol (Invitrogen, Carlsbad, CA, USA). Gene expressions patterns were examined by Affymetrix Mouse Gene 1.0 ST arrays. Total RNA was extracted by using Trizol reagent (Life technologies, Carlsbad, CA, US) with the standard operating steps given by the manufacturer. The integrity of RNA samples were checked by an Agilent Bioanalyzer 2100 (Agilent technologies, Santa Clara, CA, US), which performed as a RIN number.
Then, qualified total RNA was further purified by RNeasy micro kit (QIAGEN, GmBH, Germany) and RNase-Free DNase Set (QIAGEN, GmBH, Germany). Purified total RNA were amplified, labeled and purified by using Ambion WT Expression Kit (Ambion, US) and GeneChip WT Terminal Labeling Kit (Affymetrix, Santa Clara, CA, US). After that, array hybridization was in process through GeneChip Hybridization, Wash and Stain Kit (Affymetrix, Santa Clara, CA, US) in Hybridization Oven 645 (Affymetrix, Santa Clara, CA, US). Next was washing arrays in the Fluidics Station 450 (Affymetrix, Santa Clara, CA, US). All of these steps above were followed by their special instructions.
In the end, array slides were scanned by GeneChip Scanner 3000 (Affymetrix, Santa Clara, CA, US). At the same, Quantity control of microarray was tested by Command Console Software 3.1 (Affymetrix, Santa Clara, CA, US) with default settings. Raw data was normalized by Robust Multi-Chip Average (RMA) algorithm. All data have been deposited in NCBI's Gene Expression Omnibus (Qian et al., 2014) and are accessible through the GEO Series accession number GSE62128 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62128).

Quantitative Real-Time PCR
RNA extraction was performed all the same as the steps in DNA microarray test. cDNA was obtained by reversing transcription of purified RNA samples by using PrimeScript RT reagent A relative quantitative analysis method was used to calculate the fold change of differential expression between experimental treatment and control, as well as that between m-g and 2-g. Messenger RNA-specific oligonucleotide primers were designed by primer premier 5 or NCBI primer pick tools, and their sequences were available in table 8, together with their annealing temperatures.

Bioinformatics analysis
We got four groups of comparisons: m-g v.s. control, 2-g v.s. control, 1-g v.s. control, and m-g v.s. 2-g. iReport online software (Ingenuity Systems, USA) was used to identify the differentially expressed genes (DEGs) of each comparison. The analysis technique for filtering genes was LIMMA [37]. The filtering standard was a fold change cutoff of 1.5, with statistical significance of P < 0.05. Cellular locations and molecular functions of genes were mapped to ingenuity Expression Profiling of LG-HMF on Osteocytes knowledgebase through iReport. Bio-data mining processes were also executed by iReport based on ingenuity knowledgebase, which consisted of biological processes, pathways, diseases and interactions. The likelihood of the association between genes and given pathway, biological process, or disease was measured by Fisher's exact test with the statistical significance P < 0.05. DAVID online resource was used to cluster the DEGs. Genes were firstly mapped to three different bio-data categories, Gene ontology, SP-PIR-KEYWORDS and KEGG-PATHWAY. Then genes were clustered according to the corresponding category terms by DAVID.

Statistical Analysis
Statistically significant differences were determined by Prism statistical software (GraphPad Software Inc., LaJolla, CA, USA). A value of P < 0.05 was considered significant in all cases. All data averages or means are accompanied by SDs to indicate the amount of variability in the data.