Metabonomic-Transcriptome Integration Analysis on Osteoarthritis and Rheumatoid Arthritis

Purpose This study is aimed at exploring the potential metabolite/gene biomarkers, as well as the differences between the molecular mechanisms, of osteoarthritis (OA) and rheumatoid arthritis (RA). Methods Transcriptome dataset GSE100786 was downloaded to explore the differentially expressed genes (DEGs) between OA samples and RA samples. Meanwhile, metabolomic dataset MTBLS564 was downloaded and preprocessed to obtain metabolites. Then, the principal component analysis (PCA) and linear models were used to reveal DEG-metabolite relations. Finally, metabolic pathway enrichment analysis was performed to investigate the differences between the molecular mechanisms of OA and RA. Results A total of 976 DEGs and 171 metabolites were explored between OA samples and RA samples. The PCA and linear module analysis investigated 186 DEG-metabolite interactions including Glycogenin 1- (GYG1-) asparagine_54, hedgehog acyltransferase- (HHAT-) glucose_70, and TNF receptor-associated factor 3- (TRAF3-) acetoacetate_35. Finally, the KEGG pathway analysis showed that these metabolites were mainly enriched in pathways like gap junction, phagosome, NF-kappa B, and IL-17 pathway. Conclusions Genes such as HHAT, GYG1, and TRAF3, as well as metabolites including glucose, asparagine, and acetoacetate, might be implicated in the pathogenesis of OA and RA. Metabolites like ethanol and tyrosine might participate differentially in OA and RA progression via the gap junction pathway and phagosome pathway, respectively. TRAF3-acetoacetate interaction may be involved in regulating inflammation in OA and RA by the NF-kappa B and IL-17 pathway.


Introduction
Osteoarthritis (OA) and rheumatoid arthritis (RA) are all inflammatory joint diseases [1]. Globally, approximately 250 million people (3.6% of the population) have OA [2]. Meanwhile, RA will affect about 24.5 million people by the year 2015 [3,4]. Although some factors including cytokines and chemokines help distinguish OA from RA [5], the unclear differences in the molecular mechanisms underlying these two diseases impede the choice of the optimal clinical treatment strategy.
Metabolites have various functions including signaling, defense, and interactions with other organisms [6]. Numerous studies have shown that metabolites are associated with the pathological process of OA and RA [7][8][9]. A previous study shows that the expression levels of reactive oxygen metabolites are upregulated in patients with knee OA [10]. Zhang et al. indicated that a total of 14 metabolites extracted could be potentially used as biomarkers for OA [11]. Moreover, as a biomarker of in vivo mast cell activation, the activity of prostaglandin D2 metabolites is closely associated with the progression of RA [12]. In an animal model, Jahreis et al. showed that mold metabolites drive RA in mice via the promotion of T cells [13]. However, the interpretation of metabonomic data is difficult due to the obstacle in data extraction and disease correlation analysis [14]. A previous study indicates that the expression of interleukin 1 can regulate the progression of both OA and RA through direct stimulation of synoviocytes and augmentation of matrix degradation [15]. Actually, the detection and interpretation of metabolite-transcript coresponses using combined profiling can yield important information on the complex biological regulation mechanism of the disease [16]. Thus, integrative analysis of transcriptome and metabonomic data may contribute to a further understanding of OA and RA progression.
In the present study, a metabonomic-transcriptome integration analysis was performed. To investigate the genes and metabolites that differentially expressed between OA and RA, the differentially expressed genes (DEGs) and metabolites were revealed from microarray data and metabolite expression profile data, respectively. Then, principal component analysis (PCA) and the linear model were used to explore the DEG-metabolite interactions. Finally, based on these interactions, pathway analysis was performed on DEGassociated metabolites to reveal the pathways that participate in the process of OA and RA. This study was expected to investigate the role of key metabolites and genes as well as their interactions in the pathogenesis of OA and RA, and further understand the complex biological regulatory mechanisms of metabolites in these two diseases.

Data
Preprocessing. The normalization for transcriptome data was performed using the Robust Multichip Average (RMA) [17] method in the Affy package (version: 1.56.0) [18] of R (version: 3.4.3) software. The normalization process in this study included background adjustment, quantile normalization, and finally summarization and log base 2 scale. If different probes mapped to the same miRNA (miRNA symbol), the mean value of different probes was considered as the final expression value of this miRNA. Meanwhile, the metabolomic data could be directly read from the processed metabolite data file using R software.
2.3. The Investigation for DEGs. The P value between OA samples and RA samples in transcriptome data was calculated by the Linear Models for Microarray Data (version: 3.34.9, limma) package [19] in R software. Then, P < 0:05 was selected as the threshold for the identification of DEGs. Then, based on Euclidean distance, the bidirectional hierarchical clustering for DEMs was performed by pheatmap software (version: 1.0.8) [20]. The results were visualized using a heat map.

Principal Component
Analysis. In the current study, the average value of each gene in the transcriptome data was calculated and ranked from high to low (deleting the last 10% of the genes). Meanwhile, the proportion of the deletion value in the expression value of each metabolite was counted, and the metabolites with more than 80% deletion value were deleted. Then, the principal component analysis (PCA) was performed on the data in two groups.   [21]. The computation formula is as follows: where "m" and "g" in the formula represent metabolite abundance and gene expression level, respectively; "p " in the formula represents phenotype (OA samples vs. RA samples); "ðg : pÞ" in the formula represents the association between gene expression and phenotype; and "ε" in the formula represents normal distribution. Then, the difference of correlation coefficients between the two groups ðjr OA − r RA jÞ > 1 and P value < 0.001 were selected as the cut-off values for DEG-metabolite interaction investigation.
2.6. Metabolic Pathway Enrichment Analysis. The clusterProfiler software (version: 3.2.11) [22] is an online tool that provides enrichment analyses including KEGG [23]. Based on the P value of DEG-metabolite interactions, the KEGG pathway enrichment analysis was used to investigate pathways enriched by the DEGs associated with metabolites. P value (the significance threshold of the hypergeometric test) < 0.05 was chosen as the cut-off criterion for the present enrichment analysis.

DEG and Metabolite Investigation.
After preprocessing, a total of 171 metabolites from metabonomic data were enrolled for further investigation. Meanwhile, a total of 20,192 genes were obtained from 54,675 probes in the cur-rent transcriptome data. Among these 20,192 genes, a total of 416 upregulated genes (such as hedgehog acyltransferase (HHAT)) and 669 downregulated genes (such as Glycogenin 1 (GYG1), Unc-51 Like Kinase 3 (ULK3), and breakpoint cluster region protein (BCR)) were revealed between OA samples and RA samples in transcriptome data. The heat map for all these DEGs is shown in Figure 1.

Discussion
Although metabolites are proven to be associated with the pathological process of OA and RA, the difficulty in interpreting metabonomic data impedes the understanding of the differences of the molecular mechanisms between these two diseases. The current metabonomic-transcriptome integration analysis revealed a total of 976 DEGs and 171 metabolites between OA samples and RA samples. The PCA and  The heat map: OA-the correlation between differentially expressed genes and metabolites in the osteoarthritis group; RA-the correlation between differentially expressed genes and metabolites in the rheumatoid arthritis group; the different colors represent the correlation difference. HHAT is an enzyme in the endoplasmic reticulum that palmitoylates hedgehog proteins [24]. HHATs participate in the expression of the sonic hedgehog signaling pathway [25]. The sonic hedgehog signaling pathway can regulate the neuronal-like differentiation of bone mesenchymal stem cells [26]. It also promotes carcinoma cells associated with bone destruction [27]. A previous study shows that sonic hedgehog signaling pathway-associated factors are upregulated in synovial tissues of RA [28]. Wang et al. indicated that genes such as HHAT in the sonic hedgehog signaling pathway are novel therapeutic targets for RA [29].  Interestingly, sonic hedgehog alleviates the inhibitory effects of high glucose on the osteoblastic differentiation of bone marrow stromal cells [30]. Glucose deficiency is closely related with the development of OA [31]. A clinical and epidemiological survey shows that the plasma glucose concentration is associated with the severity of OA [32]. Data from the OA investigation indicates that glucose homeostasis can influence the risk of knee OA. Moreover, GYG1 is an enzyme that is involved in the biosynthesis of glycogen [33]. It plays a role in glycogen metabolism regulation and in the maximal glycogen level attainment in skeletal muscle [34]. A previous study shows that the inhibition of GYG synthase kinase attenuates the glucocorticoid-induced bone loss [35]. A study of hepatic glycogen catabolism and glycogen levels in rats with chronic arthritis shows that GYG expression level is closely related with arthritis progression [36]. Meanwhile, a previous study reveals that asparagine-linked glycosylation of bone morphogenetic protein is required for secretion and osteoblast differentiation [37]. In a clinical study, Tong et al. indicated that asparagine levels were all below the lower limit of quantification in the bone marrow of leukemia patients [38]. In this study, the PCA and linear module analysis showed that GYG1-asparagine and HHAT-glucose were two of the most outstanding DEG-metabolites. Thus, we speculate that GYG1-asparagine and HHAT-glucose interaction might be implicated in the pathogenesis of OA and RA. A previous study shows that gap junctions have an important function in the control or coordination of bone cell activity [39]. The intercellular gap junctions play a vital role in skeletal physiology and bone cell mechanosensing [40]. The metabolite ethanol decreases gap junction permeability in primary cultures from defined brain regions [41]. Actually, ethanol had been successfully used in the treatment of distal tarsal joint OA [42]. Meanwhile, Jonsson et al. indicated that ethanol could prevent the development of destructive RA [43]. Although the increased intercellular communication through gap junctions may contribute to the progression of OA [44], so far, no study has shown that the gap junction pathway is involved in RA progression. Furthermore, a previous study indicated that phagosomelysosome fusion participates in the biological function of bone marrow macrophages [45]. Based on a fluorescence and electron microscope study, Zuckerfranklin showed that phagosomes participated in the regulation of synovial fluid leukocytes during RA development [46]. A previous study showed that some metabolites including tyrosine participated in the phagosome pathway [47]. An expression profile analysis of tyrosine genes in human OA proves the relationship between tyrosine expression and OA progression [48]. Meanwhile, Murakami et al. indicated that the tyrosine kinase promotes the development of RA through the activation of macrophages [49]. In this study, KEGG pathway analysis showed that the DEG-associated metabolites between OA samples and RA samples such as ethanol and tyrosine were mainly enriched in pathways like the gap junction and phagosome, respectively. Thus, we speculated that the gap junction pathway enriched by the ethanol and phagosome pathway enriched by tyrosine might participate differentially in OA and RA progression. (1) Adenosine_170   [53]. Zhang et al. indicated that overexpression of microRNA-125b could facilitate inflammation in RA by activating the NF-κB signaling pathway [54]. In this study, metabolite acetoacetate was significantly enriched in the NF-kappa B and IL-17 signaling pathway, and TRAF3-acetoacetate_35 interaction was identified. TRAF3 is an intracellular protein that belongs to the TNF receptorassociated factor protein family which is implicated in NF-kappa B activation [55,56]. Reportedly, the Tengmei Decoction could improve inflammatory injury of synovium in collagen-induced arthritis rats probably by regulating the TRAF3/NF-κB signaling pathway [57]. Liu et al. showed that miR-671-3p played a crucial role in the pathogenesis of OA by targeting TRAF3 and regulating chondrocyte apoptosis and inflammation [58]. Thereby, we speculated that the NF-kappa B and IL-17 pathway enriched by acetoacetate was implicated in regulating inflammation in OA and RA probably by targeting TRAF3. However, there were some limitations in the current study including small sample size and lack of verification analysis. Thus, a further verification study based on a large sample size is needed to confirm all speculations in this study.

Conclusions
In conclusion, genes such as HHAT, GYG1, and TRAF3, as well as metabolites including glucose, asparagine, and acetoacetate might be implicated in the pathogenesis of OA and RA. Metabolites like ethanol and tyrosine may participate differentially in OA and RA progression via the gap junction pathway and the phagosome pathway, respectively. TRAF3acetoacetate interaction may be involved in regulating inflammation in OA and RA by the NF-kappa B and IL-17 pathway.

Data Availability
The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.