Biological Functions of Selenoprotein Glutathione Peroxidases (GPXs) and their Expression in Osteoarthritis

Introduction

Osteoarthritis (OA) is a disease with joint degenerative changes as the main clinical manifestation, causing pathological changes in articular cartilage. Worldwide, more than 300 million people were reported to be affected by OA.¹ OA mostly occurs in the elderly, especially after the age of 50. Among them, the incidence in women is higher than that in men.² The social cost of OA may account for 0.25% ~ 0.50% of a country’s GDP.³ OA has a high disability rate and seriously affects the quality of life of the elderly population, and it represents a major public health challenge.^4,5 OA has been reported to be associated with a number of risk factors such as age, joint injury, genetic predisposition, gender and obesity.⁶ Studies have also revealed that OA was associated with metabolic syndrome and vascular diseases.^7,8 The pathological process of OA was very complex and involved multiple signaling pathways, including Wnt/β-cat-enin, TGF-β and BMP, Indian Hedgehog, FGF, NF-κB, and Notch pathways.⁹ Meanwhile, oxidative stress and synovial inflammation have been shown to be important mechanisms in the development of OA.^10,11

Selenium is an essential micronutrient for mammals and humans, and it plays the roles of antioxidant, immune regulation and toxin antagonism in the body. The lack of selenium can lead to the occurrence of various diseases.^12,13 The primary biological function of selenium is mediated through selenoproteins. 25 selenoproteins have been found in the human genome so far. The glutathione peroxidase (GPX) family is one of the important components of the selenoprotein family. It is an important active oxygen free radical scavenger in the body, which can maintain the balance of oxygen metabolism in the tissue.^14,15 There are eight GPXs isoforms have been identified in the GPX family, namely GPX1-GPX8. Among them, GPX1-GPX4 and GPX6 were called GPXs containing selenocysteine because they used selenocysteine as their catalytic site. While GPX5, GPX7 and GPX8 used cysteine instead of selenocysteine and were called non-selenocysteine GPXs.¹⁶

Oxidative stress has been proved to be an important mechanism for the occurrence and development of OA.¹¹ Oxidative stress is a state of imbalance between reactive oxygen species (ROS) production and antioxidant defense systems. This imbalance can lead to damage to important molecules and cells, with potential consequences for the entire organism.¹⁷ ROS are involved in a variety of diseases including chronic inflammation,¹⁸ and even tumors.^19,20 In addition, studies have demonstrated that ROS are involved in inflammatory joint diseases.^21,22 Evidence has indicated that people with OA have higher levels of oxidative stress.²³ The synovial tissue of patients with arthritis has decreased levels of selenium and GPXs, decreased antioxidant capacity, and damaged chondrocyte membranes, which ultimately lead to joint degeneration.²⁴ Studies have confirmed that GPXs catalyze the reduction of H2O2 or organic hydroperoxides to water or corresponding alcohol by utilizing glutathione (GSH) as a reducing substrate. Therefore, it can further effectively reduce the level of oxidative stress in the body.^16,25

However, the association between oxidative stress in OA and the biological functions of GPXs remains unclear. The biological functions of GPXs have received more and more attention in recent years. In this study, microarray data of OA patients were used to explore GPXs expression level. Combined with the experimental study, the biological functions of GPXs and their expression in OA patients were comprehensively analyzed. We identified the core genes that GPXs mediate their biological functions, and discussed the correlation between these core genes and OA, in order to provide new ideas for the prevention and treatment of OA.

Materials and Methods

Participants and Raw Data Collection

This study explored the biological functions of five reported human selenium-containing selenoproteins GPX1, GPX2, GPX3, GPX4 and GPX6. The whole process of this study was showed in the flow chart, which was shown in Figure 1. The bioinformatic data in this study were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/gds/). These data were derived from the following datasets: GSE1919, GSE32317, GSE41038, GSE46750, GSE48556, GSE55235, GSE55457. All the above 7 datasets provided microarray data of human OA group and control group. Details of the data were listed in Supplementary Table 1. In addition, 16 OA patients and 16 healthy people from Xi’an Central Hospital were included in the experimental study of this study. According to the diagnostic criteria by American College of Rheumatology, we identified 16 patients with primary OA.²⁶ They were enrolled in the study before receiving treatment after admission to the hospital. All of the patients were selected by one surgeon according to exclusion criteria, and no blinding was applied. Exclusion criteria were patients with secondary post-traumatic OA, rheumatoid arthritis (RA), infectious and endocrine-related arthropathies, acute or chronic inflammatory diseases, malignancy, renal insufficiency, diabetes, other rheumatic diseases, congenital malformations or tumors. In addition, age-, sex-, and body mass index (BMI)-matched 16 control group members were from the same region. They were recruited by orthopedic surgeons in the orthopedics outpatient clinic. The exclusion criteria for the control group were as follows: any acute or chronic inflammatory disease, any persistent joint pain, visit to a physician for joint discomfort, diabetes, malignancy, and renal insufficiency. This study was conducted in accordance with the ethical principles contained in the current version of the Declaration of Helsinki. Informed consent was signed by all participants. This study was approved by the Medical Ethics Committee of Xi’an Central Hospital (LW-2022-016).

Figure 1 A flow chart illustrating the analysis process of this study.

GO and KEGG Pathway Enrichment Analysis of GPXs

The STRING 11.0 (https://string-db.org/) is a comprehensive and objective database that integrated all publicly available protein interaction information covering 5090 organisms and 24.6 million proteins. Its important function is to allow users to visualize gene data as interactive networks and perform gene enrichment analysis.²⁷ In this study, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed using the online software STRING to further understand the biological functions of GPXs and their associated signaling pathways. GO functional enrichment analysis, which includes molecular function (MF), biological process (BP) and cell component (CC), is one of the most important bioinformatic methods for identifying the biological functions of target genes. KEGG enrichment analysis used gene expression information to identify key pathways that were associated with target genes.^28,29 False discovery rate (FDR)-corrected P-values less than 0.05 in GO and KEGG terms were considered significantly enriched. This study presented the top 15 GO terms and the top 20 KEGG pathways in the form of a bubble plot.

Protein-Protein Interaction (PPI) Network Was Constructed

The confidence and other key parameters of the PPI networks of GPX1, GPX2, GPX3, GPX4 and GPX6 were adjusted using the online software STRING 11.0, and then the PPI network of GPXs were visualized using Cytoscape 3.8.0 software.^30,31 The importance of each protein in the network can be determined according to the interactions among various PPI nodes. The more interacting proteins a node has, the more important the node is in the network. Therefore, the core proteins can be identified by the PPI network of GPXs.

Establishment of microRNAs-GPXs Interaction Network

MicroRNAs (miRNAs) plays an important role in the regulation of bodily functions. In order to determine the upstream regulatory role of microRNAs, a miRNA-GPXs regulatory network for GPX1, GPX2, GPX3, GPX4 and GPX6 was constructed in the present study. The regulatory functions of miRNAs in GPXs were analyzed by Network analyst (https://www.networkanalyst.ca/)^32,33 and Cytoscape 3.8.0 after adjusting the filtration parameters. An interaction network was created to understand the upstream miRNA genes of GPXs.

Establishment of Transcription Factor-GPXs Interaction Network

Transcription factor (TF) is a protein molecule with a unique structure that functions to regulate gene expression. We constructed a TF-GPXs regulatory network to examine the role of TFs in the regulation of GPXs and to identify TFs that specifically regulate GPXs by Network analyst and Cytoscape 3.8.0 in this study.

Analysis of GPXs Expression in OA Patients

We conducted a literature search in the GEO public database PubMed (National Center for Biotechnology Information, NCBI) using the term “OA” obtained the “GSE1919” dataset submitted by Ungethuem U and Haeupl T et al (https://www.ncbi.nlm.nih.gov/gds/?term=GSE1919) as well as other six datasets (GSE32317, GSE41038, GSE46750, GSE48556, GSE55235 and GSE55457).^34–37 Inclusion criteria of the eligible datasets included: 1) Case-control study; 2) Diagnosis based on the OARSI classification criteria; 3) Specimen tissues including synovial tissue, peripheral blood, chondrocytes etc. from homo sapiens; Expression profiling data by array source data were extracted and transformed into logarithmic values as needed to ensure the reliability of subsequent data analysis. The mean and standard deviation of GPXs expression in the OA and the control group were calculated for each dataset.

We used Review Manager 5.3 for meta-analysis. A test for heterogeneity (Q test) was conducted on all 7 datasets before the meta-analysis. If I²<50% or P>0.05 was met, then the datasets were considered homogeneous and were analyzed using the fixed effect model. On the contrary, if the datasets were not homogeneous, the cause for heterogeneity must be identified as soon as possible and the datasets were then analyzed using the random effect model. The Combined Effect Size: standard mean difference and its 95% confidence interval (95% CI.) were subsequently calculated. P < 0.05 was considered to be statistically significant.

qRT-PCR for Validation the Results of Meta-Analysis

To validate the above results obtained by bioinformation and meta-analysis, qRT-PCR experiment was used to confirm the expression levels of main GPXs members (GPX1 and GPX4), which were randomly selected from GPX1-GPX8. The qRT-PCR experiment was conducted based on total RNA from peripheral blood mononuclear cells (PBMCs) samples of sixteen OA patients and sixteen controls. Then the total RNA was converted into cDNA with Revert Aid RT Reverse Transcription Kit (Thermo Fisher Scientific) and qRT-PCR experiment was conducted by the platform of the ABI7500 Real-Time PCR system (Applied Biosystems, Foster City, CA) according to the manufacturer’s instructions carefully. The qRT-PCR cycling conditions were as follows: 95°C for 30s, followed by 38 cycles of 95°C for 5s and 57°C for 30s.

The primer sequences of GPX1 and GPX4 were as follows:

GPX1 forward, 5’-AAGGTACTACTTATCGAGAATGTG-3’ and reverse, 5’-GTCAGCCTCGATGTCAATGGTCTG-3’;

GPX4 forward, 5’-ACAAGAACGGCTGCGTGGTGAA-3’ and reverse, 5’- GCCACACACTTGTGGAGCTAGA-3’;

β‐actin forward, 5′‐ATTGCCGACAGGATGCAGA‐3′ and reverse, 5′‐GAGTACTTGCGCTCAGGAGGA‐3′.

All the above primers were synthesized by Beijing Huada Genetic Engineering Company (Beijing, China). The relative expression levels of the two genes were calculated based on the comparative cycle threshold (Ct) equation: 2^−ΔΔCt, in this equation ΔCt = Ct (target gene) − Ct (β‐ACTIN) and ΔΔCt = Mean ΔCt_OA − Mean ΔCt_control. Two sampling t test was selected to identify the significance of differences for expression of the target genes between OA and controls.

Results

GO and KEGG Pathway Enrichment Analysis of GPXs

Enrichment analysis of a series of pathways involved in GPXs is helpful to further evaluate the molecular mechanism of GPXs in diseases. The enrichment analysis results were shown in Figure 2. The GO enrichment analysis demonstrated that the GPXs were primarily involved in biological processes including response to oxidative stress, cellular response to chemical stimulus, response to chemical and metabolic process. The molecular functions of GPXs included antioxidant activity, peroxidase activity and glutathione peroxidase activity. The cellular components of GPXs included nuclear chromatin, nucleoplasm and intracellular organelle lumen (Figure 2A). The results of KEGG pathway analysis showed that most of the terms were enriched in glutathione metabolism, hippo signaling pathway, arachidonic acid metabolism, Huntington disease and amyotrophic lateral sclerosis (Figure 2B).

Figure 2 GO and KEGG pathway enrichment analysis of GPXs. (A) GO terms for GPXs enrichment analysis. The X-axis represents gene proportions, and the Y-axis represents enrichment terms. The three shapes of circle, triangle and square represent BP, CC and MF, respectively. The size of the shapes shows the level of enrichment. (B) KEGG pathway for GPXs enrichment analysis. The X-axis represents gene proportions, and the Y-axis represents enrichment terms. Circles represent the KEGG pathway. The size of the shape shows the level of enrichment.

PPI Network of GPXs

The PPI network was shown in Figure 3. The network constructed according to the biological functions of GPXs showed that Superoxide dismutase 1(SOD1), Superoxide dismutase 2(SOD2) and catalase (CAT) interact with five proteins in the PPI network. Deleting the above 3 nodes will destabilize the entire PPI network, so they were considered the core nodes of the PPI network.

Figure 3 PPI network analysis of GPXs. Pink rhombuses represent GPXs and purple ovals represent interacting genes.

The Establishment of miRNAs-GPXs Interaction Network

The microRNA interaction network of GPXs was shown in Figure 4. 18 miRNAs (including hsa-mir-424, hsa-mir-509-3-5p, hsa-mir-15b, hsa-mir-103, hsa-mir-491-5p, hsa-mir-516a-5p, hsa-mir-107, hsa-mir-609 and hsa-mir-450-3p) were involved in the regulation of GPX1; 13 miRNAs were involved in the regulation of GPX2 (hsa-mir-325, hsa-mir-149, hsa-mir-34c-3p, hsa-mir-760, hsa-mir-629 and hsa-mir-34b), and 9 miRNAs (including hsa-mir-608, hsa-mir-524-3p, hsa-mir-498, hsa-mir-873 and hsa-mir-525-3p) were involved in the regulation of GPX3, and 6 miRNAs (including hsa-mir-371-3p, hsa-mir-346, hsa-mir-657, has-mir-324-3p, has-mir-634 and has-mir-214) were involved in the regulation of GPX4. We found that miRNA-mediated regulation was more potent for GPX1 and GPX2.

Figure 4 GPXs-miRNA interaction network. Red rhombuses represent GPXs. Navy blue ovals represent miRNA. Light blue ovals represent predicted proteins that co-regulate with miRNA.

Role of TFs in the Regulation of GPXs

The role of TFs in the regulation of GPXs was shown in Figure 5. 137 TFs were found to be involved in the regulation of the five GPXs. In particular, most TFs focus on the regulation of GPX4.

Figure 5 Regulation of GPXs by TFs. Pink rhombuses represent GPXs and purple ovals represent interacting TFs.

Analysis of Expressions of GPXs in OA Patients

The results of the meta-analysis were shown in Figure 6. Given the selected datasets in which GPX1 expression included were heterogeneous (I²=76%, P<0.05), the random effects model was used for the meta-analysis. The results of Meta-analysis showed that GPX1 expression was significantly increased in the OA group than that in the healthy controls (95% CI=0.01–1.43, P≤0.05) (Figure 6A).

Figure 6 The results of microarray-based meta-analysis. (A) Forest plot of GPX1 expression between OA group and control group. (B) Forest plot of GPX2 expression between OA group and control group. (C) Forest plot of GPX3 expression between OA group and control group. (D) Forest plot of GPX4 expression between OA group and control group. (E) Forest plot of GPX6 expression between OA group and control group.

Similarly, the selected datasets including GPX2 expression were heterogeneous (I²=79%, P<0.05), and the random effects model was used for Meta-analysis, which showed a decreasing trend of GPX2 expression in the OA group compared with the control group (95% CI=−1.48–0.02, P>0.05) (Figure 6B).

The selected GPX3 datasets were heterogeneous (I²=56%, P<0.05), and hence the random effects model was used for analysis. Meta-analysis showed that GPX3 expression was significantly lower in the OA group than that in the control group (95% CI=−0.98– −0.03, P≤0.05) (Figure 6C).

The selected GPX4 and GPX6 datasets were not heterogeneous (I²=43%, P>0.05; I²=54%, P>0.05), and fixed effects models were used for analyses. Meta-analysis showed that the results were not statistically significant. GPX4 expression showed a trend of increase in the OA group compared with the control group (95% CI=−0.17–0.39, P>0.05) (Figure 6D). However, GPX6 expression showed a trend of decrease in the OA group compared with the control group (95% CI=−0.37–0.99, P>0.05) (Figure 6).

Validation by Expression Analysis of GPXs

GPX1 and GPX4 were selected for qRT-PCR analysis to validate the above analysis, and the analysis results were shown in Figure 7. The results showed that the expression level of GPX1 in OA patients was significantly higher than that of the normal controls (Figure 7A, P<0.05). And the expression level of GPX4 of OA patients was not significantly different with that in the control group (Figure 7B, P>0.05). The results were consistent with those obtained by meta-analysis.

Figure 7 Evaluating GPXs expression in OA patients using qRT-PCR. (A) Relative expressions of GPX1 in OA patients and controls. (B) Relative expressions of GPX4 in OA patients and controls. *p < 0.05. Abbreviation: NS, Non-significant.

Discussion

OA is a progressive degenerative joint disease. The pathological landscape of OA remains unclear until now, but it is certain that various pathologies associated with OA could make it worsen in combination, including cartilage degradation, subchondral bone sclerosis, angiogenesis and nerve innervation. Current therapies for symptomatic OA include oral drug administration, intraarticular drug injection, intravenous drug administration and surgical operations. Thus, exploring the pathogenesis of OA is important for both therapeutic approaches and bioengineering research.^38,39 OA is characterized by progressive articular cartilage degeneration, osteophyte formation and subchondral bone sclerosis. Its core features are degeneration and damage of articular cartilage. Cartilage degeneration in OA is due to the disruption of the dynamic balance caused by the activation of chondrocytes by multiple factors, in which the production of matrix-degrading enzymes exceeds the ability of chondrocytes to secrete matrix components.⁴⁰ Studies have shown that changes in the microenvironment of subchondral bone may directly or indirectly affect cartilage metabolism in patients with OA.⁴¹ Oxidative stress was a well-established pathogenesis of OA. Elevated levels of ROS and oxidative stress in chondrocytes play a role in the formation of OA.¹¹ The incidence of OA patients is gradually increasing and poses a huge burden to human beings. Therefore, we explored the relationship between OA and GPXs in order to provide new ideas for etiological studies and the development of therapeutic measures.

In this present study, enrichment analysis revealed that GPXs were mainly involved in biological processes such as oxidative stress response and redox response. The molecular functions of GPXs include peroxidase activity and antioxidant activity. These findings demonstrated that GPXs play important antioxidant function in the body. ROS in the body is usually kept at a low level in cells and is crucial for maintaining cellular homeostasis and function. Oxidative stress occurs in the body when ROS production and antioxidant processes are out of balance.⁴² There is increasing evidence that oxidative stress is tightly coupled with OA pathology.²³ Chondrocyte death is thought to be a major factor in the pathogenesis of OA, and oxidative stress appears to play an important role through chondrocyte apoptosis.^40,43 As one of the three most important families of antioxidant enzymes, GPXs usually play a key role in scavenging hydrogen peroxide (H2O2), organic peroxides and lipid peroxides.¹⁶ Studies also indicated that GPX1 can be activated in the process of cartilage formation and ROS in the cartilage microenvironment can be eliminated by GPX1.⁴⁴ Meanwhile, our study found that the expression of GPX1 was elevated in OA patients. GPX1 is a key enzyme involved in the redox metabolic cycle of glutathione. At the same time, it is also one of the critical free radical scavenging enzymes in the human body. It exists in the cytoplasm and mitochondria of eukaryotic cells.⁴⁵ We speculate that the level of oxidative stress is elevated during the occurrence and development of OA, and the compensatory increase of GPX1 further plays an antioxidant role.

Studies have shown that inflammation of the synovium is an emerging pathological mechanism in OA. Synovium-synovial fluid-cartilage has been identified as the central axis of osteoarthritis pathology. A vicious circle-ROS, pro-inflammatory cytokines and cartilage degradation products that stimulate TLRs (immune responses) and activation of complementary system in OA.¹¹ Synovial inflammation and increased angiogenesis are associated with oxidative stress.^11,46 ROS and oxidative stress are elevated in patients with OA. Systemic oxidative stress induces synoviocytes and chondrocytes to produce IL-6, IL-8 and prostaglandins. In the inflammatory environment, the imbalance of synovial microvascular network and the increased energy demand of activated infiltrating immune cells and inflammatory cells lead to the decreased efficiency of synovial oxygenation, which further contributes to the hypoxic microenvironment. Hypoxia of synovial cells in a hypoxic microenvironment leads to mitochondrial damage, which in turn leads to elevated ROS levels and over-activated oxidative stress, further aggravating synovitis. The production of cytokines and growth factors by the inflamed synovium may affect the production of degrading enzymes and cartilage loss.^46,47 Studies have shown that the viscosity of synovial fluid is reduced and lipid peroxidation is increased in OA patients, suggesting a greater oxidative modification of synovial fluid components by ROS associated with OA.⁴⁸ Therefore, we speculate that there is a tendency to change the antioxidant status of synovial fluid in OA patients with increased GPXs activity. Another study showed that icariin can inhibit the sagging of synovial cells through the Xc/GPX4 axis, and further protect synovial cells. Meanwhile, it is also confirmed that the Xc/GPX4 axis is the mediator of the protective effect of icariin on synovitis.⁴⁹ These studies suggested that GPXs exhibit important antioxidant functions in humans and may be involved in the occurrence of certain bone diseases, which were consistent with the findings of this study.

Our study found that GPXs were involved in arachidonic acid metabolism. In patients with RA and OA, synovial tissue proliferates and hypertrophies due to the accumulation of neovascularization and inflammatory cells, including lymphoplasmacytic cells and macrophages. These cells were thought to be a key factor in articular cartilage and bone. Cytokines and metabolites of arachidonic acid play a potential role in this process. Arachidonic acid metabolites were able to inhibit cell proliferation in the inflamed synovium, thereby limiting arthritis inflammation and pannus formation.^50,51 The enrichment of GPXs in the arachidonic acid metabolic pathway suggested that GPXs may reduce the inflammatory response to OA treatment through the arachidonic acid pathway.

In addition, the PPI network of GPXs revealed that SOD1, SOD2, CAT and Glutathione-disulfide reductase (GSR) were important node of this network. GSR was also a key enzyme involved in the redox metabolic cycle of GSH and was a glutathione reductase. It can effectively remove excess ROS in cells.^45,52 Superoxide dismutase (SOD) and CAT were both intracellular and extracellular with several antioxidant defense lines.⁵³ Studies have shown that the levels of SOD, CAT, GPXs and other antioxidant enzymes in OA patients were reduced, and the role of oxidative stress in the pathogenesis of OA has been confirmed.⁴⁸ Collectively, the GPXs antioxidant system can regulate ROS levels by regulating related signaling pathways to maintain redox balance and the integrity of cellular components. Therefore, we speculate that antioxidant therapy will provide a promising approach for the prevention and treatment of OA caused by overexposure to ROS.

Meta-analysis of GPXs expression in OA showed that GPX1 and GPX3 expression were significantly different between the OA group and control group. In particular, the expression of GPX1 was significantly increased, which was consistent with our experimental verification results. Therefore, the high expression of GPX1 has a significant impact on the occurrence and development of OA patients. On the other hand, GPX2, GPX4 and GPX6 expression was not significant different between the OA group and control group. Our experiments also yielded similar results to the meta-analysis.

GPXs was a member of the antioxidant system. The decreased levels of GPX2, GPX3 and GPX6 in OA patients suggested that the body’s antioxidant system capacity was reduced in the OA environment, further providing direct evidence that oxidative stress may play a key role in the occurrence of OA.⁵⁴ In clinical practice, doctors should pay more attention to the patients with decreased expression of GPX2, GPX3 and GPX6. The antioxidant properties of these patients may have a downward trend, which may affect the treatment of diseases to a certain extent, and even affect the prognosis of patients. The clinical applications of the GPXs and whether they can serve early biomarkers for OA development will be continue to be studied and validated. However, the results of this study have further expanded our understanding of the relationship between the GPXs and OA, as well as the biological functions, signaling pathways related to GPXs and their expressions in OA patients. It provided new evidence for better practice of preventing and treating OA in the future by GPXs.

There are some limitations in our present study. First, some of the results of this study were derived from bioinformatics analysis, and these results need to be further screened and verified by independent experiments. Second, the expression of GPXs in OA patients was obtained by meta-analysis. If on the basis of meta-analysis, the expressions of GPXs in OA patients were detected by qRT-PCR with large population, the results of the present study will be more comprehensive.

In future studies, our results will further be screened and validated by independent experiments. The implication of key members of GPXs in OA occurrence and development will be further explored. What’s more, the expressions of GPXs in OA patients were detected by qRT-PCR with large population. Based on the above studies, evidence-based potential biomarkers that keep pace with clinical need will be developed.

Conclusion

In conclusion, GPXs play an important antioxidant role in oxidative stress response. The expression of GPX1 was elevated in PBMCs of OA patients. The changes in GPXs in patients with OA provided direct evidence that oxidative stress may play a key role in the development of OA. The changes of GPXs in OA patients may regulate the level of oxidative stress, which may influence synovial lesions and chondrocyte apoptosis.