ANGPTL4 May Regulate the Crosstalk Between Intervertebral Disc Degeneration and Type 2 Diabetes Mellitus: A Combined Analysis of Bioinformatics and Rat Models

Introduction

Low back pain (LBP) is one of the top ten injury causes affecting disability-adjusted life-years.¹ This has been reported to have an enormous socioeconomic impact.² Intervertebral disc degeneration (IVDD) is regarded as an important risk factor for LBP.^3–5 The intervertebral disc consists mainly of the nucleus pulposus (NP), cartilage endplates (CEP), and annulus fibrosus (AF). IVDD is characterized by progressive NP cell death, CEP calcification and sclerosis, AF fissure, breakdown of the extracellular matrix (ECM), immune cell infiltration and so on.^6–9 Based on the current lack of clear understanding of IVDD, it is necessary to further explore the underlying mechanisms.

Type 2 diabetes mellitus (T2DM) is a common metabolic disease characterized by insulin resistance, activation of inflammatory responses, accumulation of advanced glycation end-products (AGEs), and augmentation of oxidative stress. The activation of these mechanisms leads to alterations in connective tissues, including bone and cartilage, which also include the reduction of intervertebral height, degradation of CEP, and potential abnormalities in its diffusion performance.^10–13 It should also be noted that due to the regulation of bone metabolism by insulin, hyperglycemia induced by T2DM will lead to an increase in bone density and osteogenic activity, which leads to sclerosis of the bone endplate and may have detrimental effects on the nutrient supply of intervertebral discs, ultimately leading to IVDD.^14–16 The abnormal microenvironment in patients with T2DM also causes the accumulation of reactive oxygen species in IVDs, disorders of lipid metabolism and mitochondrial function, and stress-related premature senescence.^17–19 In addition, AGEs can cause calcification and sclerosis of the endplate, microstructural changes, intervertebral disc ECM degradation, and inflammatory cascade responses.^20–24 Therefore, the risk of IVDD in T2DM patients may inevitably increase, and it is necessary to understand the mechanisms of IVDD and T2DM.^25,26

Since there are various connections between IVDD and T2DM, our study aimed to find deeper links between them to provide new insights for the diagnosis and treatment of IVDD. In recent years, the development of RNA sequencing technology and bioinformatics has enabled the genetic investigation of the common etiologies of crosstalk in diseases genetically.^27–29 First, the GSE27494 and GSE20966 datasets were selected to identify differentially expressed genes (DEGs) that regulate IVDD and T2DM, respectively. Common DEGs were identified and protein interaction networks (PPIs) were constructed to select five hub genes. Hub genes were further verified and functional enrichment was performed to assess the biological processes involved. Functional correlation analysis, immune infiltration, transcription factors (TFs), miRNAs, and potential drugs estimated in the DSigDB database were performed for analysis based on hub genes. Finally, hub gene expression was verified by RNA sequencing and experiments using a rat model. The workflow of this study is illustrated in Figure 1.

Figure 1 The overall workflow of this study. Abbreviations: IVDD, intervertebral disc degeneration; DEGs, differentially expressed genes; GO, Gene Ontology; TFs, transcription factors.

Materials and Methods

Data Collection

All the included datasets were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). Keywords: “type 2 diabetes mellitus” or “intervertebral disc degeneration” were used to search for profiles of patients with T2DM and IVDD. The inclusion criteria were as follows:1) the arrays should contain both case and control groups; 2) the included datasets should contain at least three samples in each group for both IVDD and T2DM; 3) raw data could be easily acquired for deeper analysis; and 4) the included test specimens were from human sources. Finally, GSE27494, GSE20966, GSE25724, and GSE56081 were screened for further analyses (Table 1).

Table 1 Information of Included GEO Datasets Containing T2DM and IVDD Patients

Identification of Common DEGs Between IVDD and T2DM

DEGs from the GSE27494 and GSE20966 datasets were identified using the GEOquery and limma packages in R.^30,31 The cutoff criteria were set as P-value<0.05 and |fold change|≥1 transformed by “log2” algorithm to identify DEGs. The obtainment and exhibition of common DEGs were assisted by Venn diagrams using the interactive viewer Venny 2.1.³²

PPI Network Construction

Common DEGs of T2DM and IVDD were identified and imported into the STRING.³³ The visualization and integration of PPIs networks were further assisted with the Cytoscape platform (https://cytoscape.org/) for the analysis of protein and genetic interactions.³⁴

Establishment of the Topological Algorithm and Identification of Hub Nodes

In this study, highly interconnected hub nodes were identified using cytoHubba, a Cytoscape plugin.³⁵ Nodes were calculated using 11 topological algorithms with cytoHubba (https://apps.cytoscape.org/apps/cytohubba), and MCC, which has a better performance in predicting essential proteins, was used to determine hub genes.³⁵

Functional Analysis

Gene Ontology (GO) based on hub genes was analyzed using |fold change|≥1 and adjusted P-value <0.05 and visualized using clusterProfiler.³⁶ P-values were adjusted using the Benjamini–Hochberg method, and the top 10 pathways for each with cutoff criteria adjusted P-value<0.05 were exhibited.³⁰ CIBERSORT is a novel tool for analyzing biomarkers based on the expression of 22 immune cell subtypes (http://cibersort.stanford.edu/).³⁷ The calculation of the cell composition of the 22 immune cell subtype distributions was performed according to the codes downloaded from CIBERSORT. The scores were further analyzed in patients with disease and hub genes with high or low expression levels. The results were visualized using the R package ggplot.

Identification of the Correlation Between Hub Genes and Multiple Function-Related Genes

1491 genes associated with mitochondrial function were acquired from the Integrated Mitochondrial Protein Index (IMPI; https://mitominer.mrc-mbu.cam.ac.uk/release-4.0/impi.do). A total of 828 genes associated with lipid metabolism were acquired from Reactome (https://reactome.org/) and KEGG (https://www.kegg.jp/, Kyoto Encyclopedia of Genes and Genomes) databases. A total of 507 genes associated with ECM were obtained from a previous report and the KEGG database.³⁸ Intersection genes of DEGs and mitochondrial function-related genes (MFRGs), lipid metabolism-related genes (LMRGs) and ECM-related genes (ECMRGs) were identified as differentially expressed mitochondrial function-related genes (DE-MFRGs), lipid metabolism-related genes (DE-LMRGs) and ECM-related genes (DE-ECMRGs). Pearson’s correlation between hub genes and DEGs was determined using the R statistical package. The results were displayed in the form of a heatmap based on scale bar of z-score.

Recognition of Related TFs and miRNAs

TFs and miRNAs, which were identified by TRRUST (http://www.grnpedia.org/trrust/) and miRnet (https://www.mirnet.ca/miRNet/home.xhtml), constitute two major regulatory modes of gene expression, including transcription and post-transcription.^39,40 To select TFs and miRNAs, a coregulatory network was visualized using Cytoscape.

Potential Therapeutic Drug Prediction

Protein-drug interaction data from DSigDB (http://tanlab.ucdenver.edu/DSigDB) were used to predict potential drugs for treating IVDD and T2DM based on hub genes, with the cutoff criteria of adjusted P-value < 0.05 and composite score > 1000.⁴¹

NP Primary Cells Isolation and Culture

Primary NP tissues were achieved from six-week-old male Sprague-Dawley (S-D) rats following cervical dislocation, rinsed in 75% ethanol, and digested with 1% collagenase II solution for 2 h (Beyotime Biotechnology, China). Specifically, the gelatinous NP tissues inside IVD were isolated from annular AF and hyaline CEP using a sterile scalpel. After centrifugation and suspension, the monolayer culture of cells was according to the conditions of DMEM medium with 10% fetal bovine serum (Thermo Fisher Scientific, USA) and 1% penicillin-streptomycin (Thermo Fisher Scientific, USA), room temperature of 37°C and 5% CO₂.

RNA Sequencing

NP cells were stimulated with or without TNF-α (10 ng/mL; R&D Systems, Minneapolis, MN, USA) for 24 h at 37°C and 5% CO2. Then, TRIzol reagent (Thermo Fisher Scientific, USA) was used to isolate RNA by protocol and total RNAs were subjected to further analysis via transcriptome sequencing in Wuhan Huada Gene Technology Co., Ltd. (China): KEGG pathways and volcano plot (|fold change| ≥ 0, normalized to control, Q-value <0.001) were applied for the determination of changes under stimulation of TNF-α using transcripts per kilobase million TPM for analysis, which were documented in the Mybgi platform (WuhanHuada Gene Technology, https://mybgi.bgi.com/tech/login).

RNA Extraction and Real-Time Quantitative PCR Analyses

Cells were treated with TNF-α (10ng/mL) or exogenous addition of glucose and palmitic acid (HGHF, 25.5 mmol·L⁻¹ glucose + 300 μmolL⁻¹ palmitic acid, Thermo Fisher Scientific, USA) for 24 h to mimic the microenvironment of T2DM.⁴² RNA was extracted using the TRIzol reagent. Following with the reversed transcription to first strand complementary DNAs (cDNAs) from RNAs using A cDNA Synthesis Kit (Takara Bio, Otsu, Japan) was used to reverse-transcribe cDNA from RNA. Real-time qPCR was performed using an Applied Biosystems QuantStudio 6 Flex Real-Time PCR System (Thermo Fisher Scientific, USA). Primers were designed using NCBI BLAST, and primer information is provided in Table 2. The 2^−ΔΔCT method was used to quantify the relative expression levels of target genes after normalization to β-actin.

Table 2 Primers Information and siRNA Sequences

siRNA Transfection

si-ANGPTL4 was purchased from Integrated Biotech Solutions Co. Ltd. (China). Six-well plates were used for cell inoculation and culture after 24h when the cells achieved 60–70% density. Then, 50nM of control (siNC) or siANGPTL4 was added using Lipofectamine 3000 siRNA transfection protocols. The siNC and siANGPTL4 sequences were generated as previous study.⁴³ The siRNA sequences are listed in Table 2. The efficacy of siANGPTL4 was verified by RT-qPCR.

Protein Extraction and Western Blot Analysis

To analyze the protein expression and function of ANGPTL4, NP cells were subjected to previous treatments. Total proteins were obtained by washing twice with 1 × phosphate-buffered saline (PBS) using a mixture of RIPA and inhibitors of phosphatase, and protease inhibitors. The lysate was quantified by centrifugation (12,000×g for 15 min) using BCA quantification kit (Beyotime Biotechnology, China). Samples were separated on a 4–20% SDS-PAGE gel, electroblotted, and transferred onto 0.22-μm PVDF membranes after dissolving in SDS-sample loading buffer. The PVDF membranes were blocked using 5% bovine serum albumin (BSA) at room temperature for 2 h and then incubated overnight at 4 °C with primary antibodies (1:1000 in 5% BSA–TBST). Primary antibodies against ANGPTL4 (ab196746), SOX9 (ab185966), MMP13 (ab39012), and β-actin (D6A8) were purchased from Abcam and Cell Signaling Technology. Membranes were washed thrice with TBST for 10 min and subsequently incubated with anti-rabbit or anti-mouse IgG (Cell Signaling Technology, Danvers, MA, United States) secondary antibodies (1:10,000) for 2 h at room temperature to avoid exposure to light. Immunoreactivity was detected using the LI-COR Odyssey fluorescence imaging system after washing thrice with TBST for 10 min. The OD values were measured using the ImageJ software (National Institutes of Health, United States).

Immunofluorescence

Slides were placed in a six-well plate following disinfection with 75% ethanol to culture NP cells at a density of 3.5×10⁵ per well. The NP cells were then treated with TNF-α (10ng/mL) or HGHF for 24 h after attachment and proliferation. After fixation with paraformaldehyde (4%) for 30 min, cell membranes were permeabilized with 0.25% Triton X-100 (Beyotime Biotechnology, China) for 10 min, blocked with 5% BSA for 2 h at room temperature, NP cells then incubated with the primary antibody ANGPTL4 (1:250, Abcam) overnight at 4°C. Cells were washed three times with PBS and then incubated with an Alexa Fluor 594 conjugate secondary antibody (anti-rabbit, 1:500; Cell Signaling Technology) for 1 h. Then, the cells were incubated with DAPI for 20 min. Eventually, fluorescence images were determined with a fluorescence microscope (Olympus, Inc., Tokyo, Japan). Semi-quantitative analysis of integrated optical density (IOD) was performed using Image-Pro Plus (version 6.0; Media Cybernetics, Silver Spring, MD, USA).

Surgical Procedure

Animal experiments were performed in accordance with Shanghai Ninth People’s Hospital (Approval #SH9H-2021-TK326-1). Twenty male S-D rats aged 12 weeks were kept under the appropriate preoperative and postoperative conditions. A needle puncture was used to construct an IVDD rat model.⁴⁴ Briefly, anesthesia was administered to the rats by intraperitoneal injection of pentobarbital sodium (40 mg/kg of body weight), and 75% ethanol was used to disinfect the tails. A sterile needle of 20G was utilized to penetrate the Co7/8 discs from the dorsal skin into the NP center. The IVDD model was completed by a subsequent rotation of 360° and holding for 30s in position. Eventually, re-sterilization was performed to prevent infection. 5μL saline was injected weekly into the Co7/8 discs for IVDD. Co6/7 discs were not used as controls.

Immunohistochemistry

After rat tails were transected and fixed by paraformaldehyde, coccyxes were exposed by skinning. Intervertebral discs were cut into 5 μm sections and attached to the thickness of a glass sheet by neutral resin to seal them.⁴⁵ Deparaffinization and rehydration were first performed on the sections, and 3% hydrogen peroxide (H₂O₂) was used for 30 min to block endogenous peroxidase. The primary antibodies were incubated overnight at 4°C, followed by blocking with 10% goat serum for 2h at room temperature. Next, the secondary antibody conjugated with HRP was incubated the next day for 2h at RT after three washes in PBS. The microscopic data were obtained using a Leica DM4000 B microscope (Leica Microsystems). Semi-quantitative analysis of IOD was performed using Image-Pro Plus (version 6.0; Media Cybernetics, Silver Spring, MD, USA).

Statistical Analysis

R software (Version 4.2.2) and GraphPad Prism 8.3 were used to analyze the data and construct the figures. Continuous variables were analyzed using the Mann–Whitney U-test. All experiments presented in this study were conducted three times to obtain data. Data are presented as the mean (M) ± standard deviation (SD). Significant differences between two groups were analyzed using t-test, whereas among multiple groups, one-way analysis of variance was used. Statistical significance was set at p < 0.05.

Results

Hub Genes Screening and Functional Enrichment Analysis in IVDD and T2DM Datasets

The expression matrices of the two datasets GSE27494 and GSE20966 were normalized, and the distribution trends of the box plots were straight lines (Figure 2A and B). To assess the repeatability of the data within the group, Principal component analysis (PCA) of the two datasets was performed, and the results showed good repeatability (Figure 2C and D).

Figure 2 Expression levels in GSE27494 and GSE20966 via data wrangling. (A and B) Normalized expression matrices of GSE27494 and GSE20966 demonstrated by boxplot. (C and D) PCA diagrams. (E and F) DEGs of GSE27494 and GSE20966 shown by volcano plot. Abbreviations: PCA, principal component analysis; DEG, differentially expressed genes.

With a screening threshold of |fold change|≥1 and P-value<0.05, 648 DEGs (343 upregulated and 305 downregulated) were identified in the GSE27494 dataset, and 171 DEGs (123 upregulated and 48 downregulated) were identified in the GSE20966 dataset. Volcano plots of DEGs in the above two datasets are shown in Figure 2E and F. Venn plots were also created (Figure 3A), which showed that 13 genes were differentially expressed between the two datasets (Figure 3B). Based on the common DEGs, five hub genes were identified using the Cytohubba plugin in Cytoscape (Figure 3C). To further assess the regulation of functions based on hub genes, GO was performed, which showed the involvement of immune regulation (negative regulation of immune system processes, negative regulation of leukocyte chemotaxis, negative regulation of mononuclear cell migration, negative regulation of leukocyte migration, and monocyte chemotaxis), SMAD protein signal transduction, and ECM-related structural and functional interactions (collagen-containing extracellular matrix, basement membrane, and glycosaminoglycan binding).

Figure 3 Determination of the hub genes and functional enrichment analysis. (A) Common DEGs in GSE27494 and GSE20966 demonstrated using venn plot. (B) PPI network of common DEGs. (C) Top 5 hub genes screened and analyzed by Cytoscape (Threshold as following, |fold change|≥1 and adjusted P-value <0.05). (D) GO enrichment analysis based on 5 hub genes. Arrows represent signaling pathways related to immune filtration; Box represents pathways related to protein signal transduction; Checkmarks represent pathways related to anabolism and catabolism of ECM. Abbreviations: PPI, protein interaction network; DEGs, differentially expressed genes; GO, Gene Ontology; ECM, extracellular matrix.

Identification of the Relationship Between Hub Genes and ECM, Lipid Metabolism and Mitochondrial Function Respectively

Since the progression of IVDD and T2DM is accompanied by mitochondrial dysfunction, disturbance of lipid disorders, and ECM degradation, the relationship between hub genes and related genes was analyzed. Furthermore, based on GO enrichment analysis (Figure 3D), the relationship between hub genes and ECM was analyzed.

In GSE27494, among 1419 MFRGs and 828 LMRGs, 29 and 38 genes were categorized as DE-MFRG and DE-LMRGs, respectively. 21 DE-MFRGs were upregulated, and eight were downregulated. 24 DE-LMRGs were upregulated and 14 were downregulated. Among the 507 ECMRGs, 40 were categorized as DE ECMRGs. Seventeen up-regulated and 23 down-regulated genes were identified. The results of the Pearson correlation analysis indicated that the majority of DE-MFRG, DE-LMRG, and DE-ECMRG were highly correlated with hub genes (|r| ≥ 0.5, p < 0.05), which indicated the close relationship between hub genes and mitochondria function, lipid metabolism and ECM in GSE27494 (Figure 4A–C, the scale bar including 1, 0.5, 0, −0.5 and −1 represents z score).

Figure 4 The relationship between 5 hub genes and lipid metabolism, mitochondrial function and ECM degradation respectively in dataset GSE27494. (A) Correlation between the hub genes and DE-MFRGs. (B) Correlation between the hub genes and DE-LMRGs. (C) Correlation between the hub genes and DE-ECMRGs. The scale bar including 1, 0.5, 0, −0.5 and −1 represents z score. Abbreviations: ECM, extracellular matrix; DE-MFRGs, differentially expressed- mitochondrial function related genes; DE-LMRGs, differentially expressed-lipid metabolism related genes; DE-ECMRGs, differentially expressed-ECM related genes.

In GSE20966, 1419 MFRGs, and 828 LMRGs, six and nine genes were categorized as DE-MFRGs and DE-LMRGs, respectively. Five DE-MFRGs were upregulated, and one was downregulated. Seven DE-LMRGs were upregulated, and two were downregulated. Of the 507 ECMRGs, 4 genes were categorized as DE-ECMRGs. 3 upregulated and 1 downregulated. The correlation results indicated that partial DE-MFRG, DE-LMRG, and DE-ECMRG were highly correlated with hub genes (|r| ≥ 0.5, p < 0.05) which indicated the close relationship between hub genes and mitochondria function, lipid metabolism and ECM in GSE20966 (Figure 5A–C, the scale bar including 1.5, 1, 0.5, 0, −0.5, −1 and −1.5 represents z score).

Figure 5 The relationship between 5 hub genes and lipid metabolism, mitochondrial function and ECM degradation respectively in dataset GSE20966. (A) Correlation between the hub genes and DE-MFRGs. (B) Correlation between the hub genes and DE-LMRGs. (C) Correlation between the hub genes and DE-ECMRGs. The scale bar including 1.5, 1, 0.5, 0, −0.5, −1 and −1.5 represents z score. Abbreviations: ECM, extracellular matrix; DE-MFRGs, differentially expressed- mitochondrial function related genes; DE-LMRGs, differentially expressed-lipid metabolism related genes; DE-ECMRGs, differentially expressed-ECM related genes.

Differences in Immune Characteristics Between Groups

The distribution of immune cells among the groups was assessed. The CIBERSORT results demonstrated that patients with IVDD had significantly higher levels of activated dendritic cells than those in the non-degenerate group (Figure 6A and B), and patients in the T2DM group showed no significant changes in immune filtration patterns compared with the control group (Figure 6C and D).

Figure 6 The immune filtration pattern in GSE27494 and GSE20966 datasets. (A and B) The immune filtration pattern grouping by disease in GSE27494 using boxplot and stacked bar chart. (C and D) The immune filtration pattern grouping by disease in GSE20966 using boxplot and stacked bar chart. *P < 0.05.

Furthermore, the immune characteristics between different expression of the hub genes in the datasets were analyzed. In GSE27494, patients with high expression of ANGPTL4 showed significantly higher levels of resting NK cells, patients with higher expression of CCL2 showed significantly lower levels of activated NK cells and higher levels of resting NK cells, patients with higher expression of CCN3 and THBS2 showed significantly lower expression of activated dendritic cells, and patients with higher expression of INHBA showed no significant changes in immune filtration compared to those with lower expression levels (Figure 7A–E). In addition, in GSE20966, patients with higher expression of ANGPTL4 showed higher expression of CD8⁺ T cells; patients with higher expression of CCL2 showed significantly lower expression of M0 macrophages and resting mast cells, higher expression of M2 macrophages and activated mast cells, and higher expression of INHBA showed significantly lower expression of resting memory CD4⁺ T cells; and patients with higher expression of CCN3 and THBS2 showed no significant changes in immune filtration patterns (Figure 8A–E).

Figure 7 The immune filtration of 5 hub genes in GSE27494. (A) Immune cell composition based on the expression level of ANGPTL4. (B) Immune cell composition based on the expression level of CCL2. (C) Immune cell composition based on the expression level of CCN3. (D) Immune cell composition based on the expression level of INHBA. (E) Immune cell composition based on the expression level of THBS2. *P < 0.05.

Figure 8 The immune filtration of 5 hub genes in GSE20966. (A) Immune cell composition based on the expression level of ANGPTL4. (B) Immune cell composition based on the expression level of CCL2. (C) Immune cell composition based on the expression level of CCN3. (D) Immune cell composition based on the expression level of INHBA. (E) Immune cell composition based on the expression level of THBS2. *P < 0.05, **P < 0.01.

Construction of TF-miRNA Coregulatory Network

The MiRnet and TRRUST tools were used to predict the target miRNAs and TFs of hub genes. A total of 281 target miRNAs of the five hub genes were identified.

We used the miRNet tool to predict the target miRNAs of hub genes. Finally, we obtained 295 target miRNAs and 32 TFs of five hub genes and determined 254 mRNA-miRNA pairs. Notably, SMAD3 was identified as a TF associated with ANGPTL4, which is in accordance with the GO enrichment results (Figure 3D, SMAD protein signal transduction). Based on the prediction results, a co-expression network of mRNAs, miRNAs, and TFs, comprising 282 nodes and 319 edges, was constructed using Cytoscape (Figure 9). There were 18 miRNAs regulating ANGPTL4, 42 miRNAs regulating CCL2, 0 miRNAs regulating CCN3, 131 miRNAs regulating THBS2, and 103 miRNAs regulating INHBA. There were 5 TFs regulating ANGPTL4, 22 TFs regulating CCL2, 1 TF regulating THBS2, and 4 TFs regulating CCN3 (NOV). Partial miRNAs were also coregulated by multiple hub genes. After further searching the PubMed database (https://pubmed.ncbi.nlm.nih.gov/) for literature related to miRNAs with both T2DM and IVDD to screen for those overlapping with the above 295 miRNAs, 23 miRNAs were identified in Table 3.

Table 3 Reported miRNAs Related to Both IVDD and T2DM Among 295 miRNAs

Figure 9 TF-mRNA-miRNA coregulatory network. The blue square nodes represent TFs. Yellow round nodes represent hub genes. Green oval nodes represent miRNAs.

Validation of Hub Genes in T2DM and IVDD

To validate the expression levels of hub genes in other datasets of patients with T2DM and IVDD, two microarray datasets (GSE56081 and GSE25724) were selected for analysis. The two datasets were normalized and the distribution trends of the box plots were straight lines (Figure 10A–E). To assess the repeatability of the data within the groups, PCA of GSE56081 and GSE25724 was performed, and the results showed that the repeatability of the data was good (Figure 10B–F). The volcano plots in the above two datasets are shown in Figure 10C and G, with a threshold of P-value<0.05 and |fold change|≥1. The plots show that there was a large difference in the expression levels of the control group compared with those in IVDD. The expression levels of ANGPTL4 and INHBA differed between GSE56081 and GSE25724. However, INHBA expression level of INHBA up-regulated in GSE56081 and down-regulated in GSE25724. In addition, there was no statistically significant difference in CCL2 expression levels, and we did not determine the expression of CCN3 (or NOV) in either dataset, potentially due to limitations of gene probes and sequencing techniques (Figure 10D–H).

Figure 10 Validation of 5 hub genes in datasets GSE56081 and GSE25724. (A) Normalized expression in each sample of GSE56081 and GSE25724 using boxplot. (B) PCA analysis of GSE56081. (C) Volcano plot of DEGs of GSE56081. (D) Validation of the expression level of ANGPTL4 using boxplot in GSE56081. (E) Normalized expression in each sample of GSE25724 using boxplot. (F) PCA analysis of GSE25724. (G) Volcano plot of DEGs of GSE25724. (H) Validation of the expression level of ANGPTL4 using boxplot in GSE25724. *P < 0.05, **P < 0.01. Abbreviations: PCA, principal component analysis; DEG, differentially expressed genes.

Targeted Drug Prediction

The DSigDB database was used to predict potential target drugs associated with hub genes, which may be useful for the treatment of both IVDD and T2DM. A total of 67 target drugs were identified. The top 10 predicted target drugs according to the combined scores are shown in Table 4. A study has shown that chromium chloride could reduce lipid accumulation and improve insulin sensitivity.^46,47 The most notable among these compounds was pepstatin, for which scientific evidence supports its effectiveness as a novel predictive drug for treating both IVDD and T2DM. Among these predictions, pepstatin, for which there is now scientific evidence, supports its effectiveness as a novel therapeutic compound for treating both IVDD and T2DM.^48–53

Table 4 Top 10 Predicted Drugs Based on 5 Hub Genes

The Expression and Function of ANGPTL4 Under the Environment of IVDD and T2DM

To confirm the expression and function of ANGPTL4 in IVDD and T2DM, we mimicked the in vitro disease microenvironment of IVDD and T2DM stimulated by TNF-α and HGHF using RNA sequencing, RT-qPCR, Western blotting, and immunofluorescence. After stimulation with TNF-α or HGHF, there was a significant upregulation of ANGPTL4, based on volcano plot analysis (Figure 11A). KEGG results indicated that the DEGs that participated in cell fate pathways varied significantly, including cell cycle signaling pathway, cellular senescence, apoptosis, and DNA replication pathway, and pathways associated with IVDD and diabetic microenvironments, such as focal adhesion and AGE-RAGE signaling pathway in diabetic complications (Figure 11B). RT-qPCR showed upregulation of ANGPTL4 after stimulation with TNF-α and HGHF (Figure 11C). In addition, the Western blotting and immunofluorescence results were consistent with the results of this study (Figure 11D, E, I and 11J). Knockdown of ANGPTL4 using siRNA showed high efficacy (Figure 11F). Under the stimulation of TNF-α and HGHF, in vitro knockdown of ANGPTL4 showed significant rescuing effects on ECM metabolism, as demonstrated by the decrease in MMP13 and increase in SOX9 (Figure 11G and H). Immunohistochemical staining showed that ANGPTL4 expression in the puncture-induced rat IVDD model was in line with the results of our in vitro experiments (Figure 11K and L).

Figure 11 Experimental validation of ANGPTL4 under T2DM and IVDD. (A) Volcano plot of the expression level of NP cells stimulated by TNF-α via RNA sequencing. (Threshold as following, |fold change| ≥ 0, normalized to control, Q-value <0.001) (B) KEGG enrichment analysis based on the expression level of genes via RNA sequencing. (C–E) Quantification of ANGPTL4 in the RNA and protein level by RT-qPCR and Western blotting stimulated by HGHF and TNF-α. (F) Quantification of the efficacy of siANGPTL4 assessed by RT-qPCR. (G and H) Quantification of ECM degradation-related genes assessed by Western blotting treated by siANGPTL4 under the stimulation of TNF-α and HGHF in vitro. (I and J) Representative images and quantification of ANGPTL4 detected by immunofluorescence stimulated by HGHF and TNF-α. (K and L) Representative images and quantification of ANGPTL4 by immunohistochemistry in puncture-induced rat IVDD model, n=6. All data are presented as mean ± SD. *P < 0.05, **P < 0.01, ***P<0.001 and ****P<0.0001, n = 3. Arrows represent signaling pathways related to ECM degradation and T2DM in IVDD; T2DM, type 2 diabetes mellitus. Abbreviations: IVDD, intervertebral disc degeneration; NP, nucleus pulposus; HGHF, exogenous addition of glucose and palmitic acid; ECM, extracellular matrix.

Discussion

A close relationship between IVD and DM was reported by Hirsh et al in 1984.⁵⁴ As research on T2DM and IVDD continues, numerous studies have demonstrated a mutual correlation between T2DM and IVDD. A meta-analysis in 2016 reported that diabetes (OR 1.19, 95% CI 1.06–1.32) served as a significant risk factor in recurrent lumbar disc herniation in a systematic review including 17 studies.⁵⁵ In 1998, Robinson et al confirmed that human diabetic IVD tissues showed a reduction in matrix sulfated GAG content and glycosylation rate.⁵⁶ Researches based on Experiments have revealed that some genes such as CBP, PGC1α, Runx2, ChREBP, and P300 participate in the progression of IVDD and T2DM. However, this only partially explains the genetic involvement and lacks thorough analysis.^57,58 Furthermore, diabetes-related products, such as amyloid protein, leptin, and glucagon-like peptide-1 (GLP-1, a safe polypeptide hormone for treating T2DM), are also related to IVDD. However, the exact molecular mechanisms underlying crosstalk between IVDD and T2DM remain unclear. Our study is the first to explore the commonly shared gene characteristics of IVDD and T2DM to optimize the strategies of detection, diagnosis, and treatment.

A deeper investigation of the roles of T2DM has made it possible to better understand the potentially interactive pathogenesis of T2DM and IVDD. In addition to obesity and secondary mechanical overloading, which are characterized as classic types of T2DM, inflammation, hyperglycemia, and advanced glycation end products (AGEs) induced by T2DM could accelerate the progression of IVDD. The IVD consists of the central NP, CEP, and surrounding AF in a relatively hypoxic, acidic, hyperosmotic, and avascular microenvironment. During the progression of degeneration, changes in the microenvironment can significantly influence IVD cell state. For instance, it is possible for NPCs in the physiological state to be transformed into catabolic activities in certain environments, such as stimulation by AGEs.²⁴ Previous studies have demonstrated that the dysregulation of the IVD microenvironment could have direct disadvantages.

Multiple pathways participate in coordinating the various steps in T2DM progression and IVD. The PI3K/AKT pathway is reported to induce insulin impairment, and its dysregulation accelerates T2DM.⁵⁹ Activation of the PI3K/AKT pathway can alleviate intervertebral disc degeneration.^60,61 TXNIP signaling is another newly discovered target that is a novel bridge between IVDD and DM.^62–64 We propose that the common pathogenesis of IVDD and DM may be due to similar responses of genes and regulators to changes in the microenvironment.

To determine the hub genes, we examined two cohorts by differential gene expression analysis and identified 13 common genes between T2DM and IVDD. Interestingly, we found that the five hub genes might play a major role in the progression of IVDD and T2DM. Among the five hub genes, ANGPTL4 was finally chosen as the hub gene after validation at both the bioinformatics and experimental levels. ANGPTL4 encodes a glycosylated secreted protein containing a C-terminal fibrinogen domain that functions as a serum hormone to regulate glucose homeostasis, lipid metabolism, and insulin sensitivity. Other angiopoietin-like proteins (ANGPTLs), such as ANGPTL3, ANGPTL4, and ANGPTL8, play important roles in regulating the activity of lipoprotein lipase (LPL), a key enzyme involved in the hydrolysis of plasma lipoproteins.^65–68 Kersten et al provided an overview of ANGPTL4 and suggested a link between ANGPTL4 and lipid metabolism.⁶⁹ Furthermore, it has been found to be involved in many types of musculoskeletal diseases, such as rheumatoid arthritis, osteoarthritis, and osteoporosis, especially in the progression of IVDD.^70–74 Another study confirmed that genetic inactivation of ANGPTL4 in ANGPTL4-deficient mice improved insulin sensitivity and glucose homeostasis, thereby reducing the risk of T2DM.⁷⁵ Moreover, it was also reported that ANGPTL4 overexpression could alter the function of mitochondrial activities and modulate methionine metabolic cycles in liver tissues.⁷⁶ Similar to previous reports, our results also showed an interaction between ANGPTL4 and mitochondrial activity, lipid metabolism, and ECM degradation based on GO results and heatmap correlation analysis.^69,71

The bioinformatic analysis results were further confirmed by experiments. According to a previous study, palmitic acid (PA) is the most abundant saturated free fatty acid.⁷⁷ High glucose and free fatty acid levels are the main hallmarks of the diabetic microenvironment. In the following experiments, we aimed to mimic the diabetic microenvironment using exogenous supplementation with glucose and PA. NP cells were subjected to the treatment with or without TNF-α (10 ng/mL; R&D Systems, Minneapolis, MN, USA) for 24 h at 37°C and 5% CO2, after which the cells were isolated using TRIzol and subjected to RNA sequencing. The sequencing results demonstrated that ANGPTL4 expression was significantly upregulated in IVDD. Focal adhesion (ECM-related changes) and the AGE-RAGE signaling pathway in diabetic complications, based on KEGG functional enrichment, also confirmed their involvement. In addition, using RT-PCR, Western blotting, and immunofluorescence, our results further confirmed the upregulation of ANGPTL4 in degenerate or diabetic microenvironments by mimicking the conditions of degeneration and diabetes using TNF-α (10ng/mL) and HGHF (25.5 mmol·L⁻¹ glucose + 300 μmol·L⁻¹ palmitic acid, PA) respectively.⁴² These results show that ANGPTL4 could serve as a common mediator in the IVDD and T2DM pathways.

Cibersort was introduced in 2015 as a powerful tool better than others for characterizing the cell composition in tissues from expression profiles.³⁷ In this study, the pattern of immune cell distribution in the control, degenerate, and T2DM samples did not show multiple differences. We observed the involvement of dendritic cell activation in IVDD progression. Dendritic cells are multifunctional cells characterized by their ability to present antigens, thereby initiating inflammation-related dyshomeostasis.⁷⁸ Previous studies also reported that the proportion of dendritic cells was significantly upregulated in degenerate or herniated IVD tissues than in normal tissues, possibly due to the involvement of inflammation and adaptive immune.^79–81 Further analysis based on different expression levels of ANGPTL4 showed downregulation of NK cell activation and upregulation of resting NK cells in the low expression ANGPTL4 group. This phenomenon has also been validated by immune-related studies associated with IVDD, demonstrating a close relationship between ANGPTL4 and immune filtration in the progression of IVDD.^79,82,83 In the development of T2DM, there was no significant change in the immune filtration pattern in our results, which may be due to the sample size or metabolic mechanisms of T2DM rather than the immune disorders of T1DM.⁸⁴ Further analysis by grouping ANGPTL4 by expression levels demonstrated that CD8⁺ T cells were significantly upregulated in T2DM patients. However, the distribution of CD8⁺ T cells varies among different tissues in T2DM based on previous reports. Okano et al reported that CD8+ T cells did not change significantly in T2DM patients.^85,86 The analysis of peripheral blood samples from patients with T2DM showed a significant decrease in CD8⁺ cells.^87–90 Stentz et al reported the activation of CD8⁺ lymphocytes in muscle tissues of T2DM patients.⁹¹ Our results showed that upregulation of CD8+ cells grouped by ANGPTL4 in diabetic islets might be due to the migration of CD8⁺ T cells and involvement of islet inflammation.⁹² Overall, these results demonstrated that ANGPTL4 plays an important role in the distribution of immune cells.

Further searches of the PubMed database for studies related to miRNAs and T2DM/IVDD identified 35 miRNAs. Among these, miR-495-3p, for example, plays a pathogenic role in diabetic nephropathy and targeted NUP160.⁹³ For IVDD, miR-495-3p attenuates apoptosis and inflammation by targeting IL5RA.⁹⁴ There was significant downregulation of miR-335-3p and secondary abnormal cell subtypes in both T2DM and IVDD, which might shed light on the therapeutics for diabetes and degeneration by targeting mir-335-3p.^93,95 Although many studies related to T2DM or IVDD miRNAs have been performed, there is a gap in research related to the crosstalk between them, which requires more meaningful studies.

In step therapy for IVDD, drug treatment in the early stage is of great significance for relieving LBP and avoiding the risk of surgery. Considering the potentially huge socioeconomic burden caused by surgery, medical administration combined with health education could greatly improve the symptoms of patients to a certain extent and greatly reduce the risk of surgery, thus returning to normal life. For T2DM, drug therapy such as metformin has been greatly promoted as a classic method. In our study, we predicted several target drugs for T2DM and IVDD based on hub genes, among which pepstatin has been proven to have therapeutic effects in the improvement of insulin resistance and inhibition of IVD apoptosis.^48,53 Pepstatin is a highly efficient cell-permeable aspartyl protease inhibitor.⁹⁶ In addition, pepstatin could improve the biosynthesis of collagen and restore the binding capacity of plasma proteins in treating various complications of T2DM, such as wounded skin and deterioration of plasma composition and function.^50,52,97 Taken together, these studies indicate that pepstatin has huge potential for the treatment of T2DM and IVDD. Based on the results of our study and previous reports, we introduced a novel application of pepstatin.

Our study had some limitations. First, the limitation of clinical information documented in GEO and the possible contamination that induces bias in the final outcomes could not be neglected. Second, further studies, especially in humans, and more functional studies are needed for validation. Finally, clinical trials with a large sample size are needed to verify whether ANGPTL4 is a good indicator of T2DM and IVDD.

Conclusion

In conclusion, our research proposed common gene signatures to demonstrate the potential underlying mechanisms of IVDD and T2DM and identified ANGPTL4 as an immune-associated biomarker and a possible therapeutic target for both IVDD and T2DM.