MGST1 May Regulate the Ferroptosis of the Annulus Fibrosus of Intervertebral Disc: Bioinformatic Integrated Analysis and Validation

Background : The objective of this research was to identify differentially expressed genes (DEGs) related to ferroptosis in the annulus fibrosus (AF) during intervertebral disc degeneration (IDD). Methods : We analyzed gene data from degenerated and normal AF obtained from the GSE70362 and GSE147383 datasets. An analysis to determine the functional significance of the DEGs was conducted, followed by the creation of a network illustrating the interactions between proteins. We further analyzed the immune infiltration of the DEGs and determined the hub DEGs using LASSO regression analysis. Finally, we identified the hub ferroptosis-related DEGs (FRDEGs) and verified their expression levels using Real-time quantitative polymerase chain reaction (RT-qPCR), Western blot, Immunohistochemical Staining (IHC)


Introducation
Low back pain (LBP) represents a significant health concern and is a leading contributor to chronic disability worldwide [1,2].Intervertebral disc degeneration (IDD) is the main cause of LBP, and IDD is a condition that affects the spine and is characterized by the gradual deterioration of the intervertebral discs (IVD) that act as cushions between the vertebrae [2,3].This condition is a common musculoskeletal disorder that affects millions of people worldwide, particularly the elderly [4].Every IVD is composed of a central gel-like nucleus pulposus (NP), a surrounding annulus fibrosus (AF), and the cartilaginous endplates that secure the disc to the neighboring vertebral bodies [5].Intervertebral disc degeneration (IDD) typically arises from structural and biochemical alterations in the intervertebral disc (IVD), leading to a shift from an anabolic to a catabolic state.This transition is characterized by a reduction in key components of the extracellular matrix, such as Collagen II and Aggrecan, alongside an elevation in catabolic enzymes like MMP3 and ADAMTS5 [6,7].This transition impacts the production of the extracellular matrix (ECM), generation of enzymes, production of cytokines and chemokines, as well as the production of neurotrophic and angiogenic factors [8].In the development of IDD, the AF assumes a crucial function as it encircles the NP, thereby restricting its movement [9].The gradual deterioration and damage of the AF during IDD progression cause the NP to protrude, resulting in a diverse array of clinical symptoms, including LBP, stiffness, and reduced mobility [9,10].While the exact causes of IDD are not fully understood, factors such as aging, genetics, and lifestyle choices can contribute to its development [11,12].
The AF is a critical component of the IVD, providing structural support and maintaining the integrity of the disc.It primarily consists of obliquely crisscrossed fibrocartilaginous tissue that encircles the nucleus pulposus.Currently, the main approach to AF repair is surgical intervention, which involves the excision of the degenerated AF to alleviate clinical symptoms [1].However, the molecular changes occurring within the AF during the degeneration process are not yet fully understood.Recent studies suggest that TGF-β1 has the ability to alleviate autophagy and apoptosis in rat annulus fibrosus, which is prompted by oxidative stress, by means of the ERK signaling cas-cade [13], however, research into ferroptosis in AF remains insufficient.There is only a single report suggesting that tert-Butyl hydroperoxide (TBHP) may induce ferroptosis in AF cells in an autophagy-dependent manner [14].Specifically, the DEGs of ferroptosis-related genes between degenerative and normal AF has yet to be thoroughly investigated.Recent research has suggested that ferroptosis, a form of regulated cell death driven by lipid peroxidation, could be a contributing factor [15][16][17].Ferroptosis is associated with a multitude of health conditions, encompassing renal diseases, neoplastic disorders, and cardiovascular ailments [18][19][20].The accumulation of iron in the body can catalyze the production of reactive oxygen species (ROS), leading to cellular damage.The research conducted by Zhu et al. [21] underscores the significance of the interaction between USP11 and Sirt3 in the pathological process of IDD through the modulation of oxidative stress-induced ferroptosis.Discovering genes that are differentially expressed in relation to ferroptosis within the AF could shed light on the underlying mechanisms of IDD and unveil new targets for treatment.
In this study, we performed bioinformatics analysis of publicly available gene expression data to identify DEGs in degenerated AF compared to normal AF.Fig. 1 provides an exhaustive visual representation of the research process, presenting a step-by-step breakdown of the procedures followed in this cutting-edge investigation.
We then integrated DEGs from two datasets using robust rank aggregation and identified genes involved in ferroptosis.Furthermore, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed the function and signaling pathway of DEGs.The immune infiltration of DEGs was analyzed, and the hub genes were further identified by LASSO regression.At the same time, the hub DEGs and the ferroptosis data set were intersected to determine the FRDEGs.In our investigation, it has been conclusively determined that MGST1 (Microsomal Glutathione S-transferase 1) exhibits significant differential expression during the process of annulus fibrosus degeneration.Furthermore, we verified the gene expression through experiments.Our study provides a foundation for further investigating the mechanisms underlying IDD and the potential for drug screening.

Data Collection
Search the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) with the keyword "Annulus fibrosus" and find two datasets, GSE147383 and GSE70362, whose details are shown in Table 1.The dataset GSE147383 utilized the GPL570 platform, the Affymetrix Human Genome U133 Plus 2.0 Array, comprising a total of four nucleus pulposus samples and four annulus fibrosus samples, with an equal distribution of two normal and two degenerated samples for each type [22].GSE70362 uti-lizes the GPL17810 platform (Affymetrix Human Genome U133 Plus 2.0 Array) and consists of 24 AF samples and 24 NP samples [23].Based on the Thompson grade classification, normal groups were chosen from Thompson grades I, II, and I-II, while degenerative groups were chosen from grades IV and V.

Data Preprocessing and Identification of DEGs
We first acquired GSE70362 and GSE147383 datasets from the GEO database, which were then standardized using the Limma package in R (version 4.3.2,R Foundation for Statistical Computing, Vienna, Austria).A Quadrant bitmap was produced with the ggplot package.To identify DEGs, we again utilized the Limma package, setting the cutoff criteria at |log 2 (Fold-Change)| > 1 and p < 0.05.Heat maps and a volcano plot were generated using sangerbox [24].

Integrated Analysis of DEGs by RRA Method
Integrated analysis of DEGs is crucial in identifying common gene expression patterns across multiple datasets.We employed the Robust Rank Aggregation (RRA) method to integrate the DEGs identified from two independent datasets.RRA is a powerful algorithm that combines rankings from various studies to generate a comprehensive rank list, increasing the robustness of gene selection.The RRA method utilizes a probability model to combine sorted lists and identify genes that are differentially expressed across multiple datasets.This approach considers the sorting conditions of each microarray simultaneously and identifies genes that consistently rank at the top of each list.Using the R package, RobustRankAggreg, we integrated the DEGs identified from each dataset.After each group of DEGs was obtained, the RRA method was used to integrate and sequence the two groups of DEGs.

Functional Enrichment Analysis
Functional annotations of various genes in the integrated analysis were achieved through GO and KEGG enrichment assessments, elucidating their involvement in biological processes (BP), cellular components (CC), and molecular functions (MF).The KEGG approach serves to pinpoint overrepresented signaling pathways within the pool of genes that show differential expression.We perform the analysis using R packages such as clusterProfiler [25] and pathview [26], and visualize the results using http://www.bioinformatics.com.cn/.

Protein-Protein Interaction Network Analysis
Protein-protein interactions (PPIs) were analyzed using the STRING database (https://string-db.org),an online resource that offers insights into protein associations.We deemed a combined score exceeding 0.4 as statistically significant.For the visualization of the resulting PPI network, we utilized Cytoscape 3.10 (https://cytoscape.org/),which

GSE number Platform
Samples GSE147383 GPL570 2 degenerative and 2 normal annulus fibrosus GSE70362 GPL570 10 degenerative and 8 normal annulus fibrosus is adept at analyzing and illustrating networks.Additionally, to pinpoint key genes in the PPI network, the Cyto-Hubba plug-in within Cytoscape, employing its Maximum Common Subnetwork algorithm, was used.

Immune Infiltration Analysis
Upon thorough examination of the available datasets, it has been ascertained that GSE70362 encompasses a more extensive array of samples in comparison to GSE147383.To estimate the proportion of immune cells in the intervertebral disc (IVD) tissue, we performed immunoinfiltration analysis using the differentially expressed genes from the GSE70362 dataset.Specifically, we utilized the "IOBR" R package [27] to estimate the relative subsets of RNA transcripts (CIBERSORTx) associated with specific cell types.By leveraging the gene expression characteristics unique to immune cell types, we inferred the relative proportions of immune cells in the samples.To ensure the accuracy of the deconvolution algorithm, we enabled batch correction and only included samples with a CIBERSORTx p < 0.05 for subsequent analysis.This approach enabled us to uncover the differences in immune cell composition between individuals with IDD and healthy controls.

Identification of Hub DEGs by Machine Learning
We employed the LASSO regression machine learning method to identify crucial genes.The 'glmnt' R package was utilized, and a significance level of p < 0.05 was considered indicative of diagnostic value.Through feature selection, we effectively reduced the dimensionality using the LASSO regression algorithm, resulting in the identification of DEGs for constructing the diagnostic monitoring model of IDD.Subsequently, we performed receiver operating characteristic (ROC) curve analysis on the biomarkers using R software.

Construction of Primary Degenerated AF Cell Model in Rats
Three 6-week-old Sprague dawley rats (spfbiotech) were utilized to isolate AF cells from rat tail under aseptic conditions.The cells were digested using 0.2% collagenase I (Sangon Biotech) and incubated at 37 °C for 3 hours.The cultured cells were then passaged for 2-3 generations in Hy-Clone Dulbecco's Modified Eagle Medium (DMEM)/F12 (Cytiva) medium.Referring to the method of constructing the degeneration model of nucleus pulposus, we treated the annulus fibrosus cells with 200 µM H 2 O 2 for 3 hours to construct the degeneration model of annulus fibrosus [29].

Western Blotting
IVD tissue was obtained and homogenized in a loading buffer.The samples were boiled in a loading buffer for 10 minutes and subjected to 10% sodium dodecyl sulfate polyacrylamide gel electrophoresis to separate the proteins.The protein was transferred onto a nitrocellulose membrane, which was subsequently blocked with 5% non-fat milk in TBST for 1 hour at room temperature.The primary antibodies, β-actin (1:10,000 dilution, Proteintech, Wuhan, China) and MGST1 (1:1000 dilution, ABclonal, Wuhan, China), were applied to the membrane and left to incubate at 4 °C overnight.The membrane was then washed and further incubated with secondary antibodies-sheep antimouse and sheep anti-rabbit-for 1 hour at room temperature.Following this, protein bands were detected using an enhanced chemiluminescence system.

Immunofluorescence
Cells were inoculated into a 12-well plate Petri dish, washed by PBS and fixed with 4% paraformaldehyde solution.The cells were permeabilized with Triton-100, and then washed with PBS.The cells were blocked with 5% BSA solution for 30 min, and the primary antibody (Mgst1, 1:100 dilution, ABclonal, Wuhan, China) was incubated at 4 °C overnight.After PBS washing, the second antibody (Goat anti-Rat IgG, Alexa Fluor 488, ZSGB-BIO, Beijing, China) incubated at room temperature in the dark for 1 hour and examined by fluorescence microscope.

IDD Rat Model
The IDD rat model was established in Sprague-Dawley rats by AF acupuncture.In short, 1% pentobarbital sodium (100 mg/kg) was used for general anesthesia.We located the skin of the rat tail intervertebral disc with the aid of X-ray, then inserted a 21-gauge needle 3.0 mm parallel to the endplate, and then used X-ray again to determine the puncture position in the intervertebral disc.The control group chose other non-punctured positions without any treatment.Finally, the rats were examined by MRI 2 weeks after the operation.All procedures and agreements have been approved by Tianjin Nankai Hospital (No.2021-031).

Immunohistochemical Staining
The rat discs were fixed in 4% paraformaldehyde for 3 days, decalcified in 20% EDTA for 3 weeks, paraffinembedded, and then sectioned carefully to a thickness of 4 µm.After deparaffinization, antigen retrieval, and blocking with 5% goat serum, the slides were incubated with Mgst1 primary antibody (1:100, Abclonal).Post-primary antibody incubation, the slides are washed and the secondary antibody, horseradish peroxidase (HRP), is added.Following another series of washes to remove unbound antibodies, the sections were developed with DAB solution (GeneTech, Shanghai, China) and counterstained with hematoxylin.Histological images were acquired at 400× magnification using an Olympus BX63 microscope (Olympus, Tokyo, Japan) in randomly selected fields of each section.

Statistical Analysis
The data were analyzed using GraphPad Prism 8.0 (GraphPad Inc, San Diego, CA, USA) and presented as the mean ± standard deviation (SD) of four independent experiments.The normality of the data was assessed using the Shapiro-Wilk test.When the data were normally distributed, the difference between the two groups was evaluated using the Student's t-test.A p-value of less than 0.05 was considered statistically significant.

Analysis and Determination of DEGs from AF
Microarray data sets are downloaded and standardized as depicted in Supplementary Fig. 1.In the analysis of the annulus fibrosus dataset GSE70362, a total of 18,482 genes were analyzed, with the volcanic map illustrating 127 DEGs that were up-regulated and 67 that were down-regulated (Fig. 2A).The corresponding thermogram vividly displays the expression patterns of these 127 DEGs (Fig. 2C).In a parallel analysis of dataset GSE147383, a slightly larger cohort of 21,655 genes was identified, among which 368 were classified as DEGs, with an up-regulation in 216 genes and down-regulation in 152 genes, as detailed in Fig. 2B,D, respectively.

Integrated Analysis of DEGs by RRA Method
Using the R package, RobustRankAggreg, we integrated the DEGs identified from each dataset.The goal of the RRA method is to mitigate the impact of random noise and biases, resulting in a more robust and credible composite ranking.Our analysis revealed a total of 118 DEGs that  were present in both datasets.Of these genes, 68 were upregulated and 50 were downregulated (Supplementary Table 1).To visualize the expression patterns of these genes, we created a heatmap using the top 20 upregulated and top 20 downregulated genes, as determined by their RRA ranks (Fig. 3A).

PPI Network Analysis
To further identify the most important genes within the PPI network, we employed the Maximum Common Sub-network (MCS) algorithm of the CytoHubba plug-in in Cytoscape.Using this algorithm, we identified 57 candidate DEGs that were most strongly associated with the PPI network.These genes were ranked within the PPI network and their functional relationships were further analyzed to determine their potential biological significance (Fig. 3B).

GO and KEGG Analysis
The DEGs were subjected to functional enrichment analysis using the GO and KEGG databases.The enriched GO terms were classified into three categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF).Notably, GO analysis revealed significant enrichment of terms related to extracellular matrix organization, glycosaminoglycan binding, collagen-containing extracellular matrix, etc. (Fig. 4A-C).In addition to GO analysis, KEGG pathway analysis was conducted to identify enriched pathways associated with the identified DEGs.The most enriched KEGG pathways included the ECMreceptor interaction, Focal adhesion, Protein digestion and absorption (Fig. 4D).

Immune Infiltration Analysis
According to the results of CIBERSORTx, there is an obvious level of immune cell infiltration between degenerated AF and healthy samples from GSE70362 (Fig. 5A), and the increase of Resting memory CD4 + T cells is found in degenerated samples (Fig. 5B).

Identification of Hub DEGs by Machine Learning
We employed the LASSO regression machine learning method to identify crucial genes from 118 DEGs.We selected a value of 9 for λ, which represents the penalty parameter (Fig. 6A).In the coefficient diagram, variables with coefficients different from zero are selected for the target λ value (Fig. 6B).Through feature selection, we identified of 9 DEGs for constructing the diagnostic monitoring model of IDD, including IGFBP3, PMAIP1, MGST1, WNT16, CPAMD8, PRPH2, PCDHB15, SOD3.We generated ROC curves to assess the discriminatory ability of the model, which exhibited strong performance with an area under the curve (AUC) value of 0.762 (Fig. 6C).

Identification of Ferroptosis-Associated Genes
From Ferrdb database, we downloaded a total of 1090 genes that were previously identified as being associated with ferroptosis.We removed duplicates and obtain a more refined list of genes, resulting in the identification of 743 unique genes that were specifically implicated in ferroptosis.Next, we intersected this list of 743 genes with the list of 9 hub DEGs that we had obtained from our integrated analysis.This process allowed us to narrow down the list of potential ferroptosis-associated genes in AF to a specific DEGs: MGST1 (Fig. 6D).Through immunoinfiltration analysis, the distribution of MGST1 in different immune cell expression groups is mainly related to macrophage M0 (Fig. 6E).

In Vivo and in Vitro Verification of MGST1 Expression
We isolated annulus fibrosus cells from healthy SD rats and employed H 2 O 2 to induce degeneration, thus establishing a rat model of degenerated annulus fibrosus cells.This model was subsequently validated through experimental analysis.Comparing with the control group, RT-qPCR demonstrated an increase in the expression of MGST1 in degenerated annulus fibrosus cells (Fig. 7A).Consistent, Western blot and Immunofluorescence (IF) confirmed that the expression of MGST1 in degenerated AF cells increased (Fig. 7B,C).Our IHC of SD rats after acupuncture revealed a significant increase in MGST1 in the degeneration group of fine acupuncture (Fig. 7D-F).

Discussion
IDD is a complex pathological change that is characterized by the aging process and damage of the IVD.It is caused by a series of complex molecular mechanisms that ultimately lead to serious clinical symptoms [30].LBP is one of the most common clinical symptoms of IDD [31].It is not only a common reason for patients to seek medical attention, but also one of the leading causes of disability worldwide [32].The IVD is a complex structure that is composed of the NP, AF and cartilaginous endplate.The AF is a fibrous ring that surrounds the NP and provides structural support to the IVD [6].The thickness and intensity of various regions within the AF may vary, the posterior part of the IVD is thinner and lacks interweaving collagen fibers, resulting in a lower intensity compared to other regions [33].This may explain the higher incidence of disc herniation at the posterior part of the IVD, and this can cause an inflammatory response and a decrease in the height of the IVD [7].At present, there are few articles devoted to the changes in the transcriptome level of the AF during the process of IDD.However, the annulus fibrosus plays an extremely important role in the progression of IDD.Zhou et al. [34] identified the pivotal genes of IDD related to basement membrane (BMs) through bioinformatics analysis, and found out the potential molecular targets and drugs of BMs-related AF degeneration based on bioinformatics analysis and molecular methods.Zhu et al. [35] discovered that the inhibition of PP2A can decrease the apoptosis and degeneration of AF cells via the p38/MAPK pathway, making it a potential therapeutic target for IDD.In the present study, we conducted a comprehensive bioinformatics analysis to identify the DEGs of AF after IDD.
At present, the pathological mechanism leading to IDD remains to be studied, and ferroptosis is one of the pathological mechanisms leading to IDD [36,37].Ferroptosis is a newly discovered form of cell death that is characterized by the accumulation of lipid peroxidation and irondependent reactive oxygen species (ROS) [38].Many studies have shown that ferroptosis is involved in the antioxidant mechanism and regulation, including the regulation of glutathione metabolism in cells.The system Xc-glutamate/cystine antiporter is a heterodimeric protein complex composed of two subunits: Solute carrier family 7 member 11 (SLC7A11) and Solute carrier family 3 member 2 (SLC3A2).SLC7A11, which spans the membrane multiple times, plays a crucial role in controlling the exchange mechanism of the Xc-antiporter, while SLC3A2 functions to stabilize SLC7A11 and ensure its proper anchoring to the cell membrane.This complex facilitates the bidirectional transport of cysteine into the cell and glutamate out, thereby contributing to the cellular cysteine supply for biosynthetic processes [39][40][41].In addition to promoting the transmembrane exchange of external cysteine and internal glutamic acid, transporters also stimulate the uptake of cysteine into cells for production.Previous studies have shown that ferroptosis plays a crucial role in the pathogenesis of various diseases, including cancer, neurodegenerative diseases, and cardiovascular diseases [17,42].However, the role of ferroptosis-related genes in the development of AF degeneration remains poorly understood.
Through bioinformatics analysis, we found that the expression of 68 DEGs, including MGST1, was significantly upregulated in degenerative AF compared to healthy AF.To further validate our findings, we performed experiments to confirm that the expression of MGST1 in degenerative AF.Previous studies have shown that MGST1 plays an important role in maintaining redox homeostasis and protecting cells from oxidative stress-induced damage [43].Therefore, the upregulation of MGST1 in degenerated AFs is likely to reduce to the accumulation of ROS and oxidative damage, which are known to play a key role in the pathogenesis of IDD.Interestingly, with the development of IDD, MGST1, which plays an antioxidant role, should be gradually depleted.However, contrary to this idea, MGST1 actually shows a significant increase and plays an important role as an antioxidant.This compensatory mechanism may be triggered by the initial oxidative stress response, which activates various antioxidant pathways, including the upregulation of MGST1 expression.
Based on our results, we propose that increasing the expression of MGST1 may be a potential treatment for IDD.MGST1 is a member of the glutathione S-transferase (GST) family, which plays a critical role in the detoxification of lipid peroxidation products and maintenance of redox homeostasis [44].Peng and colleagues [44] demonstrated that heightened levels of MGST1 stimulate the proliferation of gastric carcinoma cells through the activation of the Akt/GSK-3β pathway.Wang and associates discovered that NMN enhances the cytochrome P450-dependent antioxidant pathway and glutathione metabolism.This is achieved by facilitating the translocation of Nrf2 to the nucleus and increasing the levels of Gstm1, Gstm2, and Mgst1, which in turn boosts GSH production and clears ROS, thus mitigating lung damage caused by silica [43].
However, to thoroughly substantiate the therapeutic potential of MGST1 modulation in IDD, there are several limitations and avenues for future research that warrant discussion.Firstly, the current study's findings would benefit from further validation through both knockdown and overexpression experiments.Such investigations would not only confirm the specificity of MGST1's role in IDD but also delineate the extent to which its modulation can influence disease progression or amelioration.Secondly, the mechanisms underlying ferroptosis induction in relation to MGST1 expression remain to be elucidated.Understanding the interplay between MGST1 and iron metabolism, lipid peroxidation, and glutathione levels will be critical in painting a complete picture of how MGST1 manipulation could mitigate IDD.

Conclusion
To identify the FRDEGs in AF, we performed bioinformatics analysis and observed the significantly high expression of MGST1 in degenerative AF.Subsequent experiments revealed the potential significance of MGST1 in the progression of IDD.It appears that MGST1 may exert its influence by regulating various processes involved in this degenerative condition.This offers a fresh perspective for the development of novel therapeutic strategies and medications.

Fig. 3 .
Fig. 3.The integrated analysis results in the identification of common DEGs and the construction of a PPI interaction network.(A) Heatmap plot of the top 20 upregulated and top 20 downregulated genes.(B) PPI network of DEGs.

Fig. 4 .
Fig. 4. The functional analysis of common DEGs was performed using GO and KEGG methods.(A-C) DEGs obtained by RRA method was enriched by GO and bubble diagram was drawn.(D) The enrichment circle diagram shows KEGG-related items.

Fig. 5 .
Fig. 5.An analysis of immune infiltration was conducted on the dataset GSE70362, yielding two visual representations.(A) Heatmap of the proportions of 22 immune cell types.(B) The stacked bar chart depicts the distribution of different immune cells in GSE70362.Resting memory CD4 + T cells, shown in yellowish-brown, make up the majority of the immune cell composition in each sample.

Fig. 6 .
Fig. 6.An analysis of immune infiltration was conducted on the dataset GSE70362.(A) The penalty parameter graph was plotted, with the logarithm of λ on the x-axis and the degree of freedom on the y-axis.(B) coefficient diagram.(C) The results of the ROC verification of LASSO analysis.(D) Venn diagram: the intersection of genes from Ferrdb database and hub genes obtained by LASSO regression analysis.(E) Immune cell type-specific expression levels of MGST1.

Fig. 7 .
Fig. 7. Establishment of a degeneration model for rat tail AF. (A) RT-qPCR revealed that the RNA expression of AF cells treated with differed from that of the control group MGST1.(B) Western blot demonstrated the disparity in MGST1 expression in AF cells treated with H2O2.(C) IF exhibited the variation in MGST1 expression in AF cells treated with H2O2, with MGST1 appearing in green and the nucleus in blue.(D) MRI displayed the degeneration model of SD rats, while X-ray indicated the puncture site.(E,F) IHC indicated a significant increase in MGST1 expression in the AF puncture group (AFP).*p < 0.05, **p < 0.01, ***p < 0.001, n = 3.