GREM1,LRPPRC and SLC39A4 as Potential Biomarkers of Intervertebral Disc Degeneration:A Bioinformatics Analysis based on Multiple Microarray and Single-cell Sequencing Data

doi:10.21203/rs.3.rs-2596952/v1

Download PDF

Research Article

GREM1,LRPPRC and SLC39A4 as Potential Biomarkers of Intervertebral Disc Degeneration:A Bioinformatics Analysis based on Multiple Microarray and Single-cell Sequencing Data

https://doi.org/10.21203/rs.3.rs-2596952/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 12 Sep, 2023

Read the published version in BMC Musculoskeletal Disorders →

You are reading this latest preprint version

Background:The issue of low back pain (LBP) has received considerable critical attention and has been a worldwide health problem. Intervertebral disc degeneration (IVDD) is always the subject of many classic studies in this field. The mechanistic basis of IVDD is poorly understood and has produced equivocal results.

Methods: Gene expression profiles (GSE34095,GSE147383) of IVDD patients together with control groups were analyzed in order to identify differentially expressed genes (DEGs) in GEO database.GSE23130 and GSE70362 were applied to validate the obtained key genes from DEGs by means of a best subset selection regression. Four machine-learning models were established to assess their predictive ability. Single-sample gene set enrichment analysis (ssGSEA) was used to profile correlation between overall immune infiltration levels with pfirmann grades and key genes. We also analyzed the upstream targeting miRNAs of key genes (GSE63492).We used single-cell transcriptome sequencing data (GSE160756) to define several cell clusters of nucleus pulposus (NP),annulus fibrosus (AF) and cartilaginous endplate (CEP) of degenerated disc and obtained the distribution of key genes in different cell clusters.

Results: By developing appropriate p-values and logFC values, we obtained a total of 6 DEGs. We validated 3 key genes (LRPPRC, GREM1 and SLC39A4) by an externally validated predictive modeling method. The ssGSEA results indicated that key genes were correlated with the infiltration abundance of multiple immune cells, such as dendritic cells and macrophages. Accordingly these 4 key miRNAs (miR-103a-3p,miR-484,miR-665,miR-107)were identified as upstream regulators targeting key genes using miRNet database and external GEO datasets. Finally, we plotted the spatial distribution of key genes in AF, CEP and NP.

Conclusions: Our study offered a new perspective to identify the creadible and effective gene therapy targets in IVDD.

Intervertebral disc degeneration

Transcriptomes

Immuno-infiltration

MicroRNAs

GREM1

SCL39A4

LRPPRC

Single-cell sequencing

Low back pain (LBP) is the leading cause of disability in the world today. According to GBD [1] (Global Burden of Disease), the average population growth rate of low back pain is approximately 7.8% annually worldwide, and the trend is accelerating with the aging of the society and the promotion of poor lifestyle habits. Intervertebral disc degeneration (IVDD) is the main cause of lower back pain at the the anatomic level, and its basic lesions include fibrous tissue degeneration, disc height reduction, fibrous annulus rupture, cartilage endplate loss, annulus fibrosus mucosus, small synovial bulge formation, and disc calcification[2].IVDD is a common chronic disease in orthopedics, and the procedures and approaches for disc degeneration are quite mature including interventional and minimally invasive endoscopic procedures. The academic consensus on this lesion lies in the difficulty of knowing the natural history of the disease and the lack of specific diagnosis and treatment options for different populations.

The basic structure of the intervertebral disc consists of the central nucleus pulposus (NP), the peripheral fibrocartilaginous annulus (AF) along with the cartilaginous endplate (CEP). A growing number of studies indicate that, in addition to mechanical loading, inflammatory substances, and nutritional disorders, genetic factors are considered to be another important underlying risk factor[3,4]. The degenerative course of the intervertebral disc has been considered to be a process with polygenic involvement, so targeted gene therapy may be proven to be an extremely promising treatment tool for IVDD in the future[4]. Recent studies suggest that immune infiltration may have a key role in the course of IVDD. Due to the unique structure of intervertebral disc tissue, the NP and AF are isolated from the host's immune system and are therefore considered to be immune-exempt organs[5,6]. Degenerated disc tissue is thought to have a specific pattern of immune cell infiltration. For example, Lan T [7] et al. analyzed the association between inflammatory response-related features and IVDD immune infiltration, revealing that IL-1β, LYN and NAMPT have potential as biomarkers and therapeutic targets for IVDD. Silva AJ [8] et al. assessed how macrophages affected IVDD cell gene expression and provided a new therapeutic target. Gene therapy targeting immune aspects offers a positive perspective in alleviating the development of IVDD, and deregulation of the disc inflammatory response is a promising strategy for the treatment of IVDD.

MicroRNAs (miRNAs) are non-coding RNAs with endogenous regulatory functions found in eukaryotes that are about 19–25 nt in length[9]. According to competing endogenouse RNA (ceRNA) theory, miRNAs can be conjugated to their targeted mRNAs to inhibit their translation or result in mRNA degradation, thus accomplishing the function of post-transcriptional regulation of gene expression[10]. Due to the prevalent effect on mRNA expression of miRNA and its potential utility as a disease biomarker or therapeutic target, miRNAs have become an influential area of basic biological and translational research, whose main application prospective focuses on disease therapy, molecular markers, and artificial miRNAs [10,11,12,13]. With the improvement of bioinformatics and related tool-based databases, miRNA sequencing technology has become a powerful approach for identifying and quantitatively resolving miRNAs, which breaks through the limitations of other measures. To date, some studies have revealed the IVDD-associated ceRNA regulatory network. Hu P [14] et al. used miRNA microarray data from GEO database to select 9 differentially expressed microRNAs in IVDD and obtained a regulatory network containing 3 up-regulated microRNA target gene pairs and 4 down-regulated microRNA target gene pairs based on TargetScan database. Cao S [15]et al. used AUC validation and lasso regression together with SVM screening features to construct key ceRNA regulatory networks related to oxidative stress and identified two pairs of ceRNA regulatory axes (PKD1-miR-20b-5p-AP000797 and CCNB1-miR-212-3p-AC079834). However, a common problem in related studies is that only miRNA-mRNA pairs are predicted, and the corresponding external validation is absent. Our study not only predicted the target miRNAs of key mRNA sets but also validated them accordingly by miRNA microarray data in GEO, leading to more credible conclusions.

In summary, IVDD intervention and treatment based on the gene-molecular level is currently a trend in this area of research, which focuses on how to screen for prevalent and important genetic markers in IVDD, while the immune microenvironment in IVDD is another immune target for IVDD treatment from the perspective of immune infiltration. In our study, this will be used as the starting point. First, we used the linear model to initially screen out the differential gene sets and then used external data to screen out the key gene sets with diagnostic value using the best subset regression method. The ssGSEA method was used to obtain related immune infiltration scores, and we established four machine learning models to demonstrate predictive ability for IVDD pfirmann grades. Finally, we analyzed the space distribution of key genes in IVDD disc tissue using single-cell transcriptome sequencing data, which will provide reliable evidence for potential future therapeutic targets for IVDD.

Acquisition and Preprocessing of Data

The mRNA and miRNA expression data of IVDD patients were obtained from Gene Expression Omnibus (GEO) database [16] (http://www.ncbi.nlm.nih.gov/geo/), which serves as a functional public genomics database thta includes high-throughput sequencing data, single-cell sequencing data, and microarray data.GSE147383 (4 young controls with 4 elderly patients) and GSE34095(3 young controls with 3 elderly patients) were for screening and enrichment analysis of DEGs. GSE23130 (23 disc tissue samples with different pfirmann grades) and GSE70362(48 disc tissue samples with different pfirmann grades) were combined into one expression matrix as an external data for regression analysis, machine learning models and ssGSEA immune infiltration analysis. GSE63492 (5 controls versus 5 patients) was subjected to microRNA expression profiling to filter out miRNAs with meaningful differential expression. GSE160756 is a single-cell sequencing data containing 3 NP,2 CEP and 2 AF samples from an IVDD patient using the 10X Genomics Kit, in which a total of 91,295 cells were sequenced. Preprocessing and quality control included the following steps: i) The “GEOQuery” [17] and “Bioconductor” [18] packages in R (4.2.1) were used to directly obtain the expression matrix, gene annotation files and corresponding clinical profiles. Probe sets were converted to gene symbols according to the platform's annotation files. Removing probe sets without corresponding gene symbols and retaining the maximum expression values for different probe sets targeting the same gene. ii) Plotting boxplots for each expression matrix to recognize batch effects. Principal component analysis (PCA) along with hierarchal clustering was carried out to identify intragroup differences and remove outlier samples based on the “ggplot2” package. iii) For the batch effect between samples, the “NormalizeBetweenArrays” function was used to remove them; For combined GSE23130 and GSE70362, the batch effect was eliminated using the “combat” function.Both functions were from the “limma” package [19]. iv) For GSE160756, we used the “Seurat” package [20] to perform the analysis and we merged the counts together of same tissue using the “harmony” package [21]. We filtered according to the distribution of gene expression (at least 300 genes detected per cell) and mitochondrial gene expression (max 10%). The data were log-normalised using the “NormalizeData” function with a scale factor of 10000 and were then normalized across all cells using the “ScaleData” function.

Identity DEGs and Functional Enrichment Analysis

The limma package [19] was used to construct a generalized linear model to screen the differential genes. For GSE147383 and GSE34095, the cut-off thresholds were |log2FC| >1 and p-value < 0.05. Heatmaps and volcano maps of DEGs were plotted using the “ggplot2”, the “pheatmap” packages; The intersection of DEGs from the two array datasets was taken to obtain the final DEG set. GO terms composed of molecular functions (MFs) and biological processes (BPs) were applied based on the “clusterProfiler” package [22] of R. Statistical significance was set at P < 0.05.

Definition of Key Genes Set Using a Best Subset Selection Regression

The external validation data came from GSE23130 and GSE70362. Since both chips are from Affymetrix, the data were combined for analysis. The expression matrix of 72 IVDD patients containing Pfirmann classification was obtained. Pfirmann grades I-II were defined as low degeneration grades, while grades III-V were considered as high degeneration grades as a primary outcome variable. Our subsequent analysis would be based on these data. Since there may be multicollinearity among the initially screened DEGs and may lead to overfitting of the model if all of them are included in the modeling, we used the best subset selection for variable screening. The basic idea consists of fitting a model to all possible combinations of predictor variables, and then screening the best model for the given variable based on some criteria (e.g. R2, corrected R2, MSE, Cp, AIC, SBIC, etc.), whose advantages included traversing all possible feature combinations, so the filtered features must be optimal. [23,24] By doing so, we selected genes with predictive value in DEGs, which was defined as the Key Genes. The R packages “glmnet” [25] and “olsrr” were used in this step.

Machine Learning Modeling

Using the key gene set as the independent variable, 4 supervised machine learning methods, namely, logistic regression, support vector model (SVM, with a nonlinear kernels radial basis function), randomForest (RF), and extreme gradient boosting (XGBoost), were applied to build predictive models. We randomly divided the combined data into a training set and a validation set at a ratio of 7:3. The training set was used to train the model, which subsequently was used for predictions using the validation set. In this model, gene expression values are used as continuous predictor variables. These supervised learning methods included nonlinear classification and linear regression. The model was established using the following R packages: “e1071”, “randomForest” and “XGBoost”. We chose a 3-fold cross-validation to train the XGBoost classifier. We used the validation set data to assess the reliability of the model based on the area under the curve (AUC) values of the receiver operating characteristic (ROC) curve with the “pROC” package.

Immune Infiltration Analysis based on ssGSEA

We used single sample enrichment analysis (ssGSEA) based on the “GSVA” package [26] to calculate the enrichment fraction of immune cells. Each enrichment fraction represents the extent to which genes in a particular gene set are up- or down-regulated in our external validation data. Marker information for immune cells from the article of Pornpimol Charoentong et al. [27], which revealed the related genetic phenotypes of immune cell infiltration. The “pheatmap” package was used to obtain the abundance matrix of immune cells, the correlation matrix of immune cells with degenerate grades and key genes and to plot heatmaps for visual analysis.

Targeted miRNA Selection and Validation

First, we identified the target miRNAs of key genes by miRNet (http://www.mirnet.ca/) database,which was developed to build miRNA–mRNA target networks. GSE63492 contains miRNA expression information of disc tissues from patients with IVDD and controls. We used the same method to screen for differentially expressed miRNAs in IVDD patients using the “limma” package. Cut-off threshold for this purpose was p-value < 0.05. Finally, we took the intersection of the above two miRNA sets and extracted miRNAs whose expression trends were opposite to those of their targeted mRNAs. It could be assumed that the miRNA-mRNA pairs that we predicted in this way were more reliable.

Cellular Fractionation and Analysis in Different Tissues of IVDD

After preliminary processing, we obtained single-cell sequencing data for NP, AF, and CEP of the degenerated discs in GSE160756 [28]. After using the uniform manifold approximation and projection (UMAP) method to confirm that the batch effects within the same tissue were initially reduced, a two-step clustering method was used to cluster the cells. Linear PCA clustering followed by UMAP clustering, and the final results were visualized using UMAP plots. Cell type annotations are defined by the “SingleR” package [29]. Visualization of the space expression distribution of key genes using violin plots and UMAP plots. We used the “monocle2” package to analyze cell lineage relationships.

The flow chart for our study is listed as follows (Fig. 1).

1. Screening of DEGs

Based on the results of PCA and hierarchical clustering (Fig 2), we removed the following samples: GSM841717 in GSE34095 and GSM4429813 in GSE147383. The boxplot results show an essentially uniform distribution within the normalized expression matrix (Fig 2). For GSE147383, we screened 299 differentially expressed genes relative to normal controls, of which 173 were up-regulated and 126 were down-regulated (Fig 3). As for GSE34095, we screened for 119 differentially expressed genes, of which 103 were up-regulated and 15 were down-regulated (Fig 3). Then we took the intersection of the both DEGs and 6 common genes (LRPPRC, GREM1, XPO1, HNRNPA2B1, UGP2, SLC39A4) were used in the following of the analysis. LRPPRC, GREM1, XPO1, HNRNPA2B1 and UGP2 were up-regulated while SLC39A4 was down-regulated.

2. Functional Enrichment of Both DEGs

GO enrichment analysis of DEGs in GSE34095 and GSE147383 was carried out. The GO terms whose p-value < 0.05 were selected and arranged in descending order according to gene counts. We listed the terms contained in the DEGs separately and visualized them using bubble plots (Fig 4). We found 3 terms in conmmon:GO:0062023, GO:0008278 and GO:0010008 which were on behalf of collagen-containing extracellular matrix, cohesin complex and endosome membrane.

3. Results of Key Genes and Corresponding External Validation

The dataset conbined by GSE23130 and GSE70362 was regarded as an independent external expression data for further analysis. The boxplot showed that the batch effects were moved mostly (Fig 5a). The 6 obtained DEGs were incorporated into the best subset regression model, and the correlation indexes of the accuracy and fit of the model with different numbers of variables were obtained (Fig 5b, Fig 5c and Fig 5d). It could be concluded that the model fit best when 3 variables were included, and had a pretty good accuracy. Finally, we selected LRPPRC, GREM1 and SLC39A4 as key genes.

The predictive capability of key genes is evaluated by 4 types of supervised machine learning methods. The ROC plots showed different results in the validation set when using different classifiers (Fig 6). The AUC value of the logistic regression classifier reached at 0.7857143 while the same for SVM, RandomForest and XGboost were 0.6938776, 0.6071429 and 0.622449 respectively, which showed that key genes had some predictive power for the defined group divided by pfirmann grades in either linear or non-linear models. We also plotted the nomogram of logistics and ranked the incnodepurity and gain values in RandomForest and XGboost, which could be considered proportional to the significance of the variable (Fig 6). It could be concluded that GREM1 had the greatest weight among all models.

4. Screening for Targeting miRNAs

An interplaying miRNA-mRNA network was built using miRNet database regarding to LRPPRC, GREM1 and SLC39A4 including 139 miRNAs and 153 nodes (Fig 7a). We followed the same method as screening for DEGs to handle with GSE63492. We deleted GSM1551029 in GSE63492 because they were identified as outlier samples in the hierarchical clustering (Fig 7b). We identified 124 miRNAs with significantly different expression (p-value <0.05). We obtained 7 common items by comparing the result predicted by miRNet with the miRNAs obtained from the dataset (Fig 7c). We selected mRNA-miRNA pairs with expression trends opposite to the target mRNA or with multiple opposite targets in the degenerate disc tissues and found statistically significant differences in the down-regulated expression of miR-665 (p-value =0.0023) and miR-107 (p-value =0.0202) in degenerating discs comparing to normal tissues (Fig 7d, Fig 7f), which were considered to be targets for GREM1 and LRPPRC. Interestingly, we found an additional fact: For LRPPRC and SLC39A4 whose expression trends were opposite, miR-484 (p-value =0.0493) and miR-103a-3p (p-value =0.0408) were thought to be regulatory to both genes (Fig 7e, Fig 7g). These 4 miRNAs may have critical modulatory roles for key genes.

5.ssGSEA Immune Infiltration Results

We obtained an infiltration score scale for immune cells based from the paper (Supplementary material 1 in supplementary materials). Heatmap revealed the level of infiltration of immune cells disc samples with different Pfirmann grades (Fig 8a). In order to provide a better contrast between the infiltration of different immune cells in the samples, we grouped the samples into low-level (Pfirmann grade I-II) and high-level (Pfirmann grade III-V) and visualized them as box plots (Fig 8b). The two plots presented here confirmed that the distribution of effector memeory CD4 T cell (p-value＜0.05), immature B cell (p-value＜0.05), type 1 helper T cell (TH1, p-value＜0.05), type 17 helper T cell (TH17, p-value＜0.05), myeloid-derived suppressor cells (MDSC, P-value <0.01), immature dendritic cell (hiDCs, p-value <0.01), macrophage (p-value＜0.01), plasmacytoid dendritic cells (pDC, p-value＜0.01) was statistically significantly different between the two groups. Another heatmap also showed the links between our key genes and immune cells (Fig 8c). The role of key genes in immune infiltration was confirmed. Correlation analysis showed that all 3 key genes were negatively correlated with the proportion of infiltrating immune cells. It is clear that the trends of GREM1 with pDC, hiDCs and SLC39A4 with immature B cells are consistent with the development of IVDD.

6. Spatial Distribution of Key Genes

We filtered a portion of the low-quality cells with the expected method and successfully removed the batch effect from the combined samples in GSE160756 (See supplementary material 2 and 3 for details about the quality control step). We automatically annotated the cell clusterings in the three UMAP maps according to the “SingleR” package (Fig 9a, Fig 9b, Fig 9c). Accordingly, we visualized the spatial expression of LRPPRC, GREM1 and SLC39A4 in these three tissues (Fig 9). It is clear that these three genes have different distributions within different cell clusters.

Using UMAP and violin plots, we found a significant spatial difference in GREM1 in chondrocytes of NP. Accordingly, we classified chondrocytes into 3 subgroups based on the expression levels of GREM1: Chondrocytes 1 to 4 were defined as “High GREM1” group, chondrocytes 5 and 6 were defined as “Medium GREM1” and chondrocytes 7 and 8 were defined as “Low GREM1”. The R package “monocle2” (v.2.18.0) was used to conduct the pseudo-time series analysis between those subgroups (Fig 10). The results suggested that the level of GREM1 and the differentiation trajectory of NP chondrocytes were generally consistent. GREM1 may exacerbate the degeneration of NP cells in IVDD. We obtained the top 50 genes that had the most significant variation with pseudotime based on this and plotted the corresponding heat map (Fig 10d). Based on the results of the clustering, the genes can be broadly classified into up-regulated and down-regulated types. We have done enrichment analysis for both types of those genes (Fig 10e, Fig 10f). It can be seen that down-regulated genes are mainly associated with Amoebiasis pathway, TGF-β signaling pathway fibronectin binding etc. while up-regulated genes are associated with ribosome metabolism, 5S rRNA binding and ubiquitin−protein transferase regulator activity and so on.

IVDD is a common degenerative disease with a complex etiology that poses many challenges to medicine and society. Exploring the potential mechanisms of IVDD and defining it with more specific and sensitive features for clinical diagnosis and early intervention is imperative [30]. As research progresses, we will not only be able to target disease markers based on sequencing data, but we can fully investigate their hidden mechanisms and operational characteristics. In our study, we used two independent GEO microarray datasets (GSE34095 and GSE147383) for the screening of DEGs and hundreds of DEGs were screened and identified. A large number of relevant entries were enriched using GO enrichment analysis. Among the overlapping GO terms were collagen-containing extracellular matrix, cohesin complex and endosome membrane. The extracellular matrix (ECM) [31,32] is defined as a complicated macromolecular network composed of collagen, proteoglycan, cohesin and several other glycoproteins that not only plays a supportive, protective and nutritional role for tissue cells, but is also closely related to basic life activities such as cell proliferation, differentiation, metabolism, recognition, adhesion and migration. Degradation of the intervertebral disc ECM is one of the main pathological changes of IVDD in several studies. [33,34] Yao M [35] et al. used Marein as a protective agent against intervertebral disc ECM degeneration and effectively resisted apoptosis of NP cells. Neidlinger-Wilke C [36] et al. revealed that ECM remodeling altered the mechanical microenvironment of the IVDD thereby compromising disc function. Endosomes are small vesicles produced by cells through endocytosis and this sorting function is essential for critical cellular functions [37]. Our enrichment results confirm that the biological processes of most DEGs are consistent with previous evidence and provide a theoretical basis for subsequent analysis.

When exploring the genetic biomarkers of diseases, previous studies often only used the microarray data from the disease group and the control group for simple difference comparison analysis. Future bioinformatics research needs to dig deeper into statistical results based on existing clinical and biologic theory. We used two other datasets (GSE23130 and GSE70362), that included pfirmann grades of IVDD to validate the DEGs. Pfirmann grading system is a common and well-established tool to evaluate degenerative abnormalities on T2-weighted sagittal MRI in daily clinical practice. It differentiates IVDD into 5 levels based on the signal of the central NP and peripheral NF and the height of the intervertebral space. We included pfirmann grades as outcome variables for assessing DEGs. On the one hand, the introduction of external data can avoid false positive or negative results. On the other hand, using pfirmann grades also makes the conclusions more relevant to clinical practice. And then, we employed the idea of variable screening of prediction models to further screen out genes with predictive value. With fewer variables included, we used a best subset selection regression and identified our final 3 key genes based on Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). The ROC plots of prediction models incorporating these 3 genes show that the maximum AUC can reach 0.78, suggesting that the key genes in our screen are predictive of the progression of the disease course of IVDD.

We analyzed the key genes based on existing databases such as GeneCards (http://www.genecards.org/). leucine rich pentatricopeptide repeat containing (LRPPRC) encodes a kind of leucine-rich protein that has multiple pentatricopeptide repeats (PPR) [38]. The role of this protein is unknown and may play a role in cytoskeletal organization, vesicular transport, or in transcriptional regulation of both nuclear and mitochondrial genes. Various studies demonstrate that high expression of LRPPRC is associated with poor prognosis in a variety of cancers, such as bladder urothelial carcinoma [39], lung cancer [40] and pancreatic cancer [41]. Maimaiti [42] et al. defined LRPPRC as a N6-methyladenosine (m6A) modification for intracranial aneurysms. Ghavami S [43] et al. reported that LRPPRC was associated with multiple neurodegenerative diseases via autophagy and apoptosis. In our study, LRPPRC is considered to be upregulated in the degenerated disc tissue and is expressed in endothelial cells and fibroblasts in addition to different types of chondrocytes. Accordingly, we speculate that LRPPRC may also be involved in the pathological changes of IVDD in an autophagic or apoptotic pathway, with the involvement of multiple cells.

Gremlin1 (GREM1) is known as a member of the transforming growth factor-β (TGF-β) signalling family encoding antagonists to bone morphogenetic protein (BMP) which may play a part in regulating organogenesis, body patterns and tissue differentiation and its related pathways include mainly angiogenesis (CST) and BMP signalling [44,45]. HKišonaitė M [46] et al. concluded that GREM1 can inhibit BMP2-mediated osteoblast differentiation from in vitro studies. Kobayashi H [47] et al. identified GREM1 and Islr as CAF-specific genes involved in BMP signaling and drived Colorectal Carcinogenesis. Shunlun Chen [48] et al. found that GREM1 promoted myeloid apoptosis and IVDD through inhibition of TGF-β-mediated Smad2/3 phosphorylation. This is coherent with our findings. The role BMP plays in IVDD is unclear though and it is characterized by a wide range and versatility of BMP and its signal pathways and inhibitors. Hollenberg AM [49] et al. demonstrated a correlation between the expression level of BMP-2 and Pfirmann grades, while Haschtmann D [50] et al. concluded that BMP-2 resulted in an up-regulation of Col I and type II, and of aggrecan gene expression. Most current studies believe that BMP can promote the ossification of AF. However, it is not clear whether it can be a unique treatment for IVDD and the effect of BMP in overall disc tissue remains to be investigated. We learned from the prediction models that GREM1 was the strongest predictor among the 3 genes so it will also be the centre of our research and discussion afterwards. Our study firstly proved the potential of GREM1 as a gene marker for IVDD, which is also in line with several previous studies [51,52]. GREM1 can affect cellular function by multiple mechanisms in addition to acting as a type of BMP inhibitor [53,54,55]. Whether GREM1 could also ameliorate the development of IVDD is another worthwhile studying field.

SLC34A4 encodes a member of the zinc/iron-regulated transporter-like protein (ZIP) family which localizes to cell membranes and is mostly required for zinc uptake in the intestine. Among its related pathways are Metal ion SLC transporters and Transport of inorganic cations/anions and amino acids/oligopeptides (http://www.genecards.org/). The main diseases associated with SLC39A4 are Acrodermatitis Enteropathica [56], Zinc-Deficiency Type Acrodermatitis [57], other than that there are relatively few studies. Zinc is an essential element for human and is involved in the maintenance of many metabolic functions. Recent studies have linked metabolic changes to cell fate. Liang J [58] et al. discovered that type II alveolar epithelial cell ZIP8 deficiency in young mice results in reduced precursor cell function and impaired self-renewal. Chen PH [59] et al. identified an unexpected role for ZIP7 in iron death by maintaining endoplasmic reticulum homeostasis, and these findings may have therapeutic implications for diseases involving iron death and zinc dysregulation. Few studies have reported the relationship between genes related to zinc metabolism and IVDD. Staszkiewicz R [60] et al. investigated the relationship between the concentration of local metal ions in the patient's disc tissue and the degree of progression of IVDD and found that the strongest relationships were noted between the concentrations of zinc. Our study suggests that down-regulation of SLC39A4 may be another important feature of IVDD and we hypothesize that abnormalities in SLC39A4 lead to an imbalance in zinc metabolism inside and outside the disc tissue cells which may promote IVDD through several pathways such as programmed cell death. The exact mechanism and principle of this is subject to further experimental verification.

Current studies on miRNAs in IVDD have confirmed that a variety of miRNAs play critical roles in the process of IVDD through apoptosis, aberrant proliferation, inflammatory response and ECM degradation [61,62,63]. Another major research interest in miRNA is its important role in exosome therapy [64]. Exosome therapy is achieved by direct in vitro injection of extracellular vesicles (EVs) containing miRNAs or by building vectors to intervene in cellular metabolism using paracrine signaling regulation [65,66]. However, caution must be exercised when studying the application of miRNAs, as miRNA-mRNA and miRNA-miRNA regulatory pathways may not be the same in different tissues [67], and misuse may lead to imbalance in ceRNA regulatory network. In this study, we successfully screened meaningful mRNA-miRNA pairs (GREM1-mir-665, LRPPRC-mir-107, LRPPRC/SLC39A4-mir-484, LRPPRC/SLC39A4-miR-103a-3p) for these 3 key genes by combining the miRNA target prediction database miRnet and external miRNA microarray dataset GSE63492. The application and validation of these miRNAs need further experimental proof.

The pattern of immune infiltration is another crucial issue. Our results on immune infiltration by ssGSEA in patients under different pfirmann grades show a decreased infiltration of effector memeory CD4 T cell, pDC,hiDCs and a increased infiltration of TH1, TH17, MDSC, macrophage and pDCs in disc tissues of high degeneration grade relative to low degeneration grade. Interesting is the simultaneous association of GREM1 with pDCs, hiDCs and SLC39A4 with immature B cells. On the one hand, infiltration of immune cells in degenerated discs can further amplify the inflammatory cascade response [68,69], and on the other hand, immune cells may lead to vascular invasion and release of neurogenic factors [70]. Our study confirmed the infiltration of some immune cells in IVDD and the possible mechanisms [71,72]. However, the role of dendritic cells (DCs) in IVDD has been poorly studied. We hypothesize that the migration and antigen-presentation functions of DCs are critical for disc tissue inflammation initiation and tolerogenic immune responses [73]. However, this must also be tested in future experimental studies.

The current study still has many limitations. For example, we depended on our own clinical experience to devide pfirmann grades into low degeneration and high degeneration groups, which was a reluctant move limited by the sample size and made the outcome variables less detailed. Secondly, the analysis of single-cell transcriptome data lacked the corresponding control group. We only used p-values to screen for significant miRNAs, which may result in a degree of false positives. We relied on available data that also lacked other clinical phenotypes of the patients, such as comorbid underlying diseases, which may have led to errors in the results of the difference analysis. In addition, in the single-cell sequencing clustering step, we divided the cells mainly into chondrocytes and non-chondrocytes. Our future studies will incorporate other omic approaches that incorporate genomics, proteomics and metabolomics, complemented by appropriate experiments, to more fully address the relevant role of IVDD biomarkers.

In conclusion, we identified 6 DEGs of IVDD by differential analysis of microarray data, and finally identified 3 key genes (GREM1, SLC39A4, LRPPRC) by external data validation. We validated the prediction models constructed using four machine learning methods: SVM, RF, XGBoost, and Logistic Regression. Finally, we analyzed the immune infiltration of key genes using the method of ssGSEA and predicted the immune infiltration pattern of IVDD in combination with pfirmann grades;We predicted the upstream miRNA of key genes using miRNet and external data and analyzed the distribution of key genes in NP, AF, CEP using single cell sequencing data. We also compared the change in GREM1 expression for the distinct transcriptional states using pseudo-time analysis. Our study offered a new perspective to identify the creadible and effective gene therapy targets in IVDD.

IVDD: Intervertebral disc degeneration

NP: Nucleus pulposus

AF: Annulus fibrosus

CEP: Cartilaginous endplate

ssGSEA: Single-sample gene set enrichment analysis

ceRNA: Competing endogenous RNAs

DEGs: Differentially expressed genes

PCA: Principal component analysis

GO:Gene Ontology

SVM: Support vector model

RF: RandomForest

XGBoost: Extreme gradient boosting

AUC: Area under curve

ROC: Receiver operator characteristic curve

UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction

LRPPRC: leucine rich pentatricopeptide repeat containing

GREM1: Gremlin1

SLC39A4: Solute Carrier Family 39 Member 4

BMP: Bone morphogenetic protein

ECM: Extracellular matrix

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The datasets analysed during the current study are available in Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/): GSE34095, GSE147383, GSE70632, GSE23130, GSE63492, GSE160756.

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by ShanXi Applied Basic Research Programme (nos.201901D111409) and ShanXi Health Care Commission Subjects (nos.2020011). All of the funding were gratefully acknowledged.

Authors' contributions

ZhaoLiang Zhangand XingHua Ji were responsible for the collection and preparation of data. LinDong Wei wrote the original draft. Jinfeng Zhang did some visualisation of the icons. Supervision and preparation of all work was the responsibility of JianZhong Huo. ZhaoLiang Zhang and XingHua Ji contributed equally to this work.

Acknowledgements

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

Institute for Health Metrics and Evaluation (IHME).GBD Compare Data Visualization. Seattle, WA:IHME,University of ashington, 2020.Available from http://vizhub.healthdata.org/gbd-compare. (Accessed [INSERT DATE])
Ge J, Cheng XQ, Yan Q, et al. Calcitonin inhibits intervertebral disc degeneration by regulating protein kinase C[J]. J Cell Mol Med, 2020, 24(15): 8650-8661.
Livshits G, Popham M, Malkin I, Lumbar disc degeneration and genetic factors are the main risk factors for low back pain in women: the UKtwin spine study. Ann Rheum Dis. 2011;70:1740-1745.
Sampara P, Banala RR, Vemuri SK, et al. Understanding the molecular biology of intervertebral disc degeneration and potential gene therapy strategies for regeneration: a review. Gene Ther. 2018 Apr;25(2):67-82. doi: 10.1038/s41434-018-0004-0. Epub 2018 Mar 22. PMID: 29567950.
Di Martino A, Merlini L, Faldini C. Autoimmunity in Intervertebral Disc Herniation: From Bench to Bedside. Expert Opin Ther Targets (2013) 17 (12):1461–70. doi: 10.1517/14728222.2013.834330
Sun Z, Liu B, Luo ZJ. The Immune Privilege of the Intervertebral Disc: Implications for Intervertebral Disc Degeneration Treatment. Int J Med Sci. 2020 Feb 24;17(5):685-692. doi: 10.7150/ijms.42238. PMID: 32210719; PMCID: PMC7085207.
Lan T, Hu Z, Guo W, et al. Development of a Novel Inflammatory-Associated Gene Signature and Immune Infiltration Patterns in Intervertebral Disc Degeneration. Oxid Med Cell Longev. 2022 Sep 22;2022:2481071. doi: 10.1155/2022/2481071. PMID: 36193061; PMCID: PMC9526649.
Silva AJ, Ferreira JR, Cunha C, et al. Macrophages Down-Regulate Gene Expression of Intervertebral Disc Degenerative Markers Under a Pro-inflammatory Microenvironment. Front Immunol (2019) 10:1508. doi: 10.3389/fimmu.2019.01508
Zhan S, Wang K, Xiang Q, et al. lncRNA HOTAIR upregulates autophagy to promote apoptosis and senescence of nucleus pulposus cells. J Cell Physiol. 2020 Mar;235(3):2195-2208. doi: 10.1002/jcp.29129. Epub 2019 Sep 3. PMID: 31478571.
Zhang P, Wu W, Chen Q, Chen M. Non-Coding RNAs and their Integrated Networks. J Integr Bioinform. 2019 Jul 13;16(3):20190027. doi: 10.1515/jib-2019-0027. PMID: 31301674; PMCID: PMC6798851.
Citron, F., Armenia, J., Franchin, G., et al. An integrated approach identifies mediators of local recurrence in head and neck squamous carcinoma. Clin. Cancer Res. 23, 3769-3780. 10.1158/1078-0432.CCR-16-2814
Calloni R, Bonatto D. Scaffolds for Artificial miRNA Expression in Animal Cells. Hum Gene Ther Methods. 2015 Oct;26(5):162-74. doi: 10.1089/hgtb.2015.043. Epub 2015 Sep 25. PMID: 26406928.
He B, Zhao Z, Cai Q, et al. miRNA-based biomarkers, therapies, and resistance in Cancer. Int J Biol Sci. 2020 Jul 19;16(14):2628-2647. doi: 10.7150/ijbs.47203. PMID: 32792861; PMCID: PMC7415433.
Hu P, Feng B, Wang G, et al. Microarray based analysis of gene regulation by microRNA in intervertebral disc degeneration. Mol Med Rep. 2015 Oct;12(4):4925-30. doi: 10.3892/mmr.2015.4022. Epub 2015 Jul 2. PMID: 26134418; PMCID: PMC4581765.
Cao S, Liu H, Fan J, et al. An Oxidative Stress-Related Gene Pair (CCNB1/PKD1), Competitive Endogenous RNAs, and Immune-Infiltration Patterns Potentially Regulate Intervertebral Disc Degeneration Development. Front Immunol. 2021 Nov 9;12:765382. doi: 10.3389/fimmu.2021.765382. PMID: 34858418; PMCID: PMC8630707.
Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 2002 Jan 1;30(1):207-10. doi: 10.1093/nar/30.1.207. PMID: 11752295; PMCID: PMC99122.
Davis S, Meltzer PS. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics. 2007 Jul 15;23(14):1846-7. doi: 10.1093/bioinformatics/btm254. Epub 2007 May 12. PMID: 17496320.
Huber W, Carey VJ, Gentleman R, et al.Orchestrating high-throughput genomic analysis with Bioconductor. Nat Methods. 2015 Feb;12(2):115-21. doi: 10.1038/nmeth.3252. PMID: 25633503; PMCID: PMC4509590.
Ritchie ME, Phipson B, Wu D, et al.limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015 Apr 20;43(7):e47. doi: 10.1093/nar/gkv007. Epub 2015 Jan 20. PMID: 25605792; PMCID: PMC4402510.
Butler A, Hoffman P, Smibert P, et al. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. 2018 Jun;36(5):411-420. doi: 10.1038/nbt.4096. Epub 2018 Apr 2. PMID: 29608179; PMCID: PMC6700744.
Korsunsky I, Millard N, Fan J, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019 Dec;16(12):1289-1296. doi: 10.1038/s41592-019-0619-0. Epub 2019 Nov 18. PMID: 31740819; PMCID: PMC6884693.
Yu G, Wang LG, Han Y, et al.clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012 May;16(5):284-7. doi: 10.1089/omi.2011.0118. Epub 2012 Mar 28. PMID: 22455463; PMCID: PMC3339379.
Goldstein BA, Navar AM, Carter RE. Moving beyond regression techniques in cardiovascular risk prediction: applying machine learning to address analytic challenges. Eur Heart J. 2017 Jun 14;38(23):1805-1814. doi: 10.1093/eurheartj/ehw302. PMID: 27436868; PMCID: PMC5837244.
Rustichini A, Conen KE, Cai X, et al.Optimal coding and neuronal adaptation in economic decisions. Nat Commun. 2017 Oct 31;8(1):1208. doi: 10.1038/s41467-017-01373-y. PMID: 29084949; PMCID: PMC5662730.
Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw. 2010;33(1):1-22. PMID: 20808728; PMCID: PMC2929880.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013 Jan 16;14:7. doi: 10.1186/1471-2105-14-7. PMID: 23323831; PMCID: PMC3618321.
Charoentong P, Finotello F, Angelova M, et al.Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017 Jan 3;18(1):248-262. doi: 10.1016/j.celrep.2016.12.019. PMID: 28052254.
Gan Y, He J, Zhu J, et al. Spatially defined single-cell transcriptional profiling characterizes diverse chondrocyte subtypes and nucleus pulposus progenitors in human intervertebral discs. Bone Res 2021 Aug 16;9(1):37.
Aran D, Looney AP, Liu L, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 2019 Feb;20(2):163-172. doi: 10.1038/s41590-018-0276-y. Epub 2019 Jan 14. PMID: 30643263; PMCID: PMC6340744.
Ohnishi T, Iwasaki N, Sudo H. Causes of and Molecular Targets for the Treatment of Intervertebral Disc Degeneration: A Review. Cells. 2022 Jan 24;11(3):394. doi: 10.3390/cells11030394. PMID: 35159202; PMCID: PMC8834258.
Karamanos NK, Theocharis AD, Piperigkou Z, et al. A guide to the composition and functions of the extracellular matrix. FEBS J. 2021 Dec;288(24):6850-6912. doi: 10.1111/febs.15776. Epub 2021 Mar 23. PMID: 33605520.
Theocharis AD, Skandalis SS, Gialeli C, et al. Extracellular matrix structure. Adv Drug Deliv Rev. 2016 Feb 1;97:4-27. doi: 10.1016/j.addr.2015.11.001. Epub 2015 Nov 10. PMID: 26562801.
RZhang XB, Hu YC, Cheng P, et al. Targeted therapy for intervertebral disc degeneration: inhibiting apoptosis is a promising treatment strategy. Int J Med Sci. 2021 May 27;18(13):2799-2813. doi: 10.7150/ijms.59171. PMID: 34220308; PMCID: PMC8241771.
Kim JH, Ham CH, Kwon WK. Current Knowledge and Future Therapeutic Prospects in Symptomatic Intervertebral Disc Degeneration. Yonsei Med J. 2022 Mar;63(3):199-210. doi: 10.3349/ymj.2022.63.3.199. PMID: 35184422; PMCID: PMC8860939.
Yao M, Zhang J, Li Z, et al. Marein protects human nucleus pulposus cells against high glucose-induced injury and extracellular matrix degradation at least partly by inhibition of ROS/NF-κB pathway. Int Immunopharmacol. 2020 Mar;80:106126. doi: 10.1016/j.intimp.2019.106126. Epub 2020 Jan 10. PMID: 31931363.
Neidlinger-Wilke C, Galbusera F, et al. Mechanical loading of the intervertebral disc: from the macroscopic to the cellular level. Eur Spine J. 2014 Jun;23 Suppl 3:S333-43. doi: 10.1007/s00586-013-2855-9. Epub 2013 Jun 21. PMID: 23793454.
Scott CC, Vacca F, Gruenberg J. Endosome maturation, transport and functions. Semin Cell Dev Biol. 2014 Jul;31:2-10. doi: 10.1016/j.semcdb.2014.03.034. Epub 2014 Apr 4. PMID: 24709024.
Tian T, Ikeda J, Wang Y, et al.Role of leucine-rich pentatricopeptide repeat motif-containing protein (LRPPRC) for anti-apoptosis and tumourigenesis in cancers. Eur J Cancer. 2012 Oct;48(15):2462-73. doi: 10.1016/j.ejca.2012.01.018. Epub 2012 Feb 10. PMID: 22326293.
Wei WS, Wang N, Deng MH, et al. LRPPRC regulates redox homeostasis via the circANKHD1/FOXM1 axis to enhance bladder urothelial carcinoma tumorigenesis. Redox Biol. 2021 Nov 27;48:102201. doi: 10.1016/j.redox.2021.102201. Epub ahead of print. PMID: 34864630; PMCID: PMC8645923.
Hu Y, Cui J, Jin L, et al. LRPPRC contributes to the cisplatin resistance of lung cancer cells by regulating MDR1 expression. Oncol Rep. 2021 Apr;45(4):4. doi: 10.3892/or.2021.7955. Epub 2021 Mar 2. PMID: 33649818.
Wang L, Luo J, Li Y, et al.Mitochondrial-Associated Protein LRPPRC is Related With Poor Prognosis Potentially and Exerts as an Oncogene Via Maintaining Mitochondrial Function in Pancreatic Cancer. Front Genet. 2022 Feb 14;12:817672. doi: 10.3389/fgene.2021.817672. PMID: 35237297; PMCID: PMC8885106.
Maimaiti A, Turhon M, Cheng X, et al. m6A regulator-mediated RNA methylation modification patterns and immune microenvironment infiltration characterization in patients with intracranial aneurysms. Front Neurol. 2022 Aug 5;13:889141. doi: 10.3389/fneur.2022.889141. PMID: 35989938; PMCID: PMC9389407.
Ghavami S, Shojaei S, Yeganeh B, et al. Autophagy and apoptosis dysfunction in neurodegenerative disorders. Prog Neurobiol. 2014 Jan;112:24-49. doi: 10.1016/j.pneurobio.2013.10.004. Epub 2013 Nov 6. PMID: 24211851.
Brazil DP, Church RH, Surae S, et al. BMP signalling: agony and antagony in the family. Trends Cell Biol. 2015 May;25(5):249-64. doi: 10.1016/j.tcb.2014.12.004. Epub 2015 Jan 12. PMID: 25592806.
O'Reilly S. Gremlin: a complex molecule regulating wound healing and fibrosis. Cell Mol Life Sci. 2021 Dec;78(24):7917-7923. doi: 10.1007/s00018-021-03964-x. Epub 2021 Nov 3. PMID: 34731251.
HKišonaitė M, Wang X, Hyvönen M. Structure of Gremlin-1 and analysis of its interaction with BMP-2. Biochem J. 2016 Jun 1;473(11):1593-604. doi: 10.1042/BCJ20160254. Epub 2016 Apr 1. PMID: 27036124; PMCID: PMC4888461.
Kobayashi H, Gieniec KA, Wright JA, et al. The Balance of Stromal BMP Signaling Mediated by GREM1 and ISLR Drives Colorectal Carcinogenesis. Gastroenterology. 2021 Mar;160(4):1224-1239.e30. doi: 10.1053/j.gastro.2020.11.011. Epub 2020 Nov 14. Erratum in: Gastroenterology. 2021 Nov;161(5):1728. PMID: 33197448.
Chen S, Lei L, Li Z,et al. Grem1 accelerates nucleus pulposus cell apoptosis and intervertebral disc degeneration by inhibiting TGF-β-mediated Smad2/3 phosphorylation. Exp Mol Med. 2022 Apr;54(4):518-530. doi: 10.1038/s12276-022-00753-9. Epub 2022 Apr 19. PMID: 35440754; PMCID: PMC9076866.
Hollenberg AM, Maqsoodi N, Phan A,et al. Bone morphogenic protein-2 signaling in human disc degeneration and correlation to the Pfirrmann MRI grading system. Spine J. 2021 Jul;21(7):1205-1216. doi: 10.1016/j.spinee.2021.03.002. Epub 2021 Mar 5. PMID: 33677096; PMCID: PMC8356724.
Haschtmann D, Ferguson SJ, Stoyanov JV. BMP-2 and TGF-β3 do not prevent spontaneous degeneration in rabbit disc explants but induce ossification of the annulus fibrosus. Eur Spine J. 2012 Sep;21(9):1724-33. doi: 10.1007/s00586-012-2371-3. Epub 2012 May 26. PMID: 22639297; PMCID: PMC3459107.
Chan SC, Tekari A, Benneker LM, et al. Osteogenic differentiation of bone marrow stromal cells is hindered by the presence of intervertebral disc cells. Arthritis Res Ther. 2015 Dec 25;18:29. doi: 10.1186/s13075-015-0900-2. PMID: 26809343; PMCID: PMC4727301.
Tekari A, May RD, Frauchiger DA,et al. The BMP2 variant L51P restores the osteogenic differentiation of human mesenchymal stromal cells in the presence of intervertebral disc cells. Eur Cell Mater. 2017 Feb 23;33:197-210. doi: 10.22203/eCM.v033a15. PMID: 28266688.
Sung NJ, Kim NH, Surh YJ, et al. Gremlin-1 Promotes Metastasis of Breast Cancer Cells by Activating STAT3-MMP13 Signaling Pathway. Int J Mol Sci. 2020 Dec 3;21(23):9227. doi: 10.3390/ijms21239227. PMID: 33287358; PMCID: PMC7730512.
Wen T, Wang H, Li Y,et al. Bone mesenchymal stem cell-derived extracellular vesicles promote the repair of intervertebral disc degeneration by transferring microRNA-199a. Cell Cycle. 2021 Feb;20(3):256-270. doi: 10.1080/15384101.2020.1863682. Epub 2021 Jan 26. PMID: 33499725; PMCID: PMC7889239.
Aashaq S, Batool A, Mir SA, et al. TGF-β signaling: A recap of SMAD-independent and SMAD-dependent pathways. J Cell Physiol. 2022 Jan;237(1):59-85. doi: 10.1002/jcp.30529. Epub 2021 Jul 19. PMID: 34286853.
Schmitt S, Küry S, Giraud M,et al. An update on mutations of the SLC39A4 gene in acrodermatitis enteropathica. Hum Mutat. 2009 Jun;30(6):926-33. doi: 10.1002/humu.20988. PMID: 19370757.
Zhong W, Yang C, Zhu L, et al. Analysis of the relationship between the mutation site of the SLC39A4 gene and acrodermatitis enteropathica by reporting a rare Chinese twin: a case report and review of the literature. BMC Pediatr. 2020 Jan 27;20(1):34. doi: 10.1186/s12887-020-1942-4. PMID: 31987033; PMCID: PMC6983971.
Liang J, Huang G, Liu X, et al. The ZIP8/SIRT1 axis regulates alveolar progenitor cell renewal in aging and idiopathic pulmonary fibrosis. J Clin Invest. 2022 Jun 1;132(11):e157338. doi: 10.1172/JCI157338. PMID: 35389887; PMCID: PMC9151700.
Chen PH, Wu J, Xu Y,et al. Zinc transporter ZIP7 is a novel determinant of ferroptosis. Cell Death Dis. 2021 Feb 19;12(2):198. doi: 10.1038/s41419-021-03482-5. PMID: 33608508; PMCID: PMC7895949.
Staszkiewicz R, Bryś K, Gładysz D, et al. Changes in Elements and Relationships among Elements in Intervertebral Disc Degeneration. Int J Environ Res Public Health. 2022 Jul 25;19(15):9042. doi: 10.3390/ijerph19159042. PMID: 35897416; PMCID: PMC9332279.
Wang C, Cui L, Gu Q, et al. The Mechanism and Function of miRNA in Intervertebral Disc Degeneration. Orthop Surg. 2022 Mar;14(3):463-471. doi: 10.1111/os.13204. Epub 2022 Feb 9. PMID: 35142050; PMCID: PMC8926997.
Ji ML, Zhang XJ, Shi PL et al. Downregulation of microRNA-193a-3p is involved in invertebral disc degeneration by targeting MMP14. J Mol Med (Berl). 2016 Apr;94(4):457-68. doi: 10.1007/s00109-015-1371-2. Epub 2015 Dec 1. PMID: 26620678.
Ding Y, Wang L, Zhao Q, et al. MicroRNA‑93 inhibits chondrocyte apoptosis and inflammation in osteoarthritis by targeting the TLR4/NF‑κB signaling pathway. Int J Mol Med. 2019 Feb;43(2):779-790. doi: 10.3892/ijmm.2018.4033. Epub 2018 Dec 18. PMID: 30569118; PMCID: PMC6317687.
Xing H, Zhang Z, Mao Q, et al. Injectable exosome-functionalized extracellular matrix hydrogel for metabolism balance and pyroptosis regulation in intervertebral disc degeneration. J Nanobiotechnology. 2021 Sep 6;19(1):264. doi: 10.1186/s12951-021-00991-5. PMID: 34488795; PMCID: PMC8419940.
Lu K, Li HY, Yang K, et al. Exosomes as potential alternatives to stem cell therapy for intervertebral disc degeneration: in-vitro study on exosomes in interaction of nucleus pulposus cells and bone marrow mesenchymal stem cells. Stem Cell Res Ther. 2017 May 10;8(1):108. doi: 10.1186/s13287-017-0563-9. PMID: 28486958; PMCID: PMC5424403.
Zhang X, Cai Z, Wu M,et al. Adipose Stem Cell-Derived Exosomes Recover Impaired Matrix Metabolism of Torn Human Rotator Cuff Tendons by Maintaining Tissue Homeostasis. Am J Sports Med. 2021 Mar;49(4):899-908. doi: 10.1177/0363546521992469. PMID: 33719604.
Matsuyama H, Suzuki HI. Systems and Synthetic microRNA Biology: From Biogenesis to Disease Pathogenesis. Int J Mol Sci. 2019 Dec 24;21(1):132. doi: 10.3390/ijms21010132. PMID: 31878193; PMCID: PMC6981965.
Takada T, Nishida K, Doita M, et al. Fas ligand exists on intervertebral disc cells: a potential molecular mechanism for immune privilege of the disc. Spine (Phila Pa 1976). 2002 Jul 15;27(14):1526-30. doi: 10.1097/00007632-200207150-00009. PMID: 12131712.
Wang HQ, Samartzis D. Clarifying the nomenclature of intervertebral disc degeneration and displacement: from bench to bedside. Int J Clin Exp Pathol. 2014 Mar 15;7(4):1293-8. PMID: 24817926; PMCID: PMC4014210.
Risbud MV, Shapiro IM. Role of cytokines in intervertebral disc degeneration: pain and disc content. Nat Rev Rheumatol. 2014 Jan;10(1):44-56. doi: 10.1038/nrrheum.2013.160. Epub 2013 Oct 29. PMID: 24166242; PMCID: PMC4151534.
Shamji MF, Setton LA, Jarvis W, et al. Proinflammatory cytokine expression profile in degenerated and herniated human intervertebral disc tissues. Arthritis Rheum. 2010 Jul;62(7):1974-82. doi: 10.1002/art.27444. PMID: 20222111; PMCID: PMC2917579.
Gorth DJ, Shapiro IM, Risbud MV. Transgenic mice overexpressing human TNF-α experience early onset spontaneous intervertebral disc herniation in the absence of overt degeneration. Cell Death Dis. 2018 Dec 18;10(1):7. doi: 10.1038/s41419-018-1246-x. PMID: 30584238; PMCID: PMC6315044.
Worbs T, Hammerschmidt SI, Förster R. Dendritic cell migration in health and disease. Nat Rev Immunol. 2017 Jan;17(1):30-48. doi: 10.1038/nri.2016.116. Epub 2016 Nov 28. PMID: 27890914.

No competing interests reported.

Supplementarymaterial1.pdf
Supplementary material 1.pdf: Marker gene dataset for immune cells: “Metagene” represents the name of marker genes. “Cell type” represents different types of immune cells. ” Immunity” represents types of immune response.
Supplementarymaterial2.pdf
Supplementary material 2.pdf: Preprocessed AF(a)(b), CEP(c)(d) NP(e)(f) quality control violin plots: The horizontal coordinates indicate different samples, which are simplified here to be expressed as numbers;The vertical coordinates indicate different types of genes:nFeature_RNA and nCount_RNA representthe total number of genes and total number of gene expressions where there is a positive correlation between the two in regular conditions. Percent.mt and percent.ribo indicate the expression ratio of mitochondrial and ribosomal genes. The percentage of mitochondrial genes was controlled to less than 10%.
Supplementarymaterial3.pdf
Supplementary material 3.pdf: UMAP and cell clustering maps before and after the removal of the batch effect regarding AF(a)(b),CEP(c)(d) and NP(e)(f) :(a)(c)(e) denotes the clustering of samples within the group before removing the batch effect, and it can be seen that there are significant differences in the distribution of different samples; (b)(d)(f) indicates the clustering of cells in each group after removing the batch effect and and performs preliminary cell clustering.

Download PDF

Journal Publication

published 12 Sep, 2023

Read the published version in BMC Musculoskeletal Disorders →

Editorial decision: Major revision
27 Apr, 2023
Reviews received at journal
26 Mar, 2023
Reviews received at journal
11 Mar, 2023
Reviewers agreed at journal
03 Mar, 2023
Reviewers agreed at journal
03 Mar, 2023
Reviewers invited by journal
02 Mar, 2023
Editor assigned by journal
25 Feb, 2023
Editor invited by journal
25 Feb, 2023
Submission checks completed at journal
25 Feb, 2023
First submitted to journal
16 Feb, 2023

You are reading this latest preprint version

GREM1,LRPPRC and SLC39A4 as Potential Biomarkers of Intervertebral Disc Degeneration:A Bioinformatics Analysis based on Multiple Microarray and Single-cell Sequencing Data

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Methods

Acquisition and Preprocessing of Data

Identity DEGs and Functional Enrichment Analysis

Definition of Key Genes Set Using a Best Subset Selection Regression

Machine Learning Modeling

Immune Infiltration Analysis based on ssGSEA

Targeted miRNA Selection and Validation

Cellular Fractionation and Analysis in Different Tissues of IVDD

Results

1. Screening of DEGs

2. Functional Enrichment of Both DEGs

3. Results of Key Genes and Corresponding External Validation

4. Screening for Targeting miRNAs

5.ssGSEA Immune Infiltration Results

6. Spatial Distribution of Key Genes

Discussion

Conclusions

Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and materials

The datasets analysed during the current study are available in Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/): GSE34095, GSE147383, GSE70632, GSE23130, GSE63492, GSE160756.

Competing interests

Funding

Authors' contributions

Acknowledgements

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1