The gene expression profile of dataset GSE113212 was obtained from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/), National Center for Biotechnology Information (NCBI). This data set was based on the GPL17077 platform and contained a total of 8 samples, including 4 hypertrophic ligamentum flavum samples from the elderly individuals and 4 non-hypertrophic samples from the young individuals (Table 1). We performed standard procedures of quality control on the dataset via the R package array QualityMetrics (Kauffmann et al., 2009), and log2 normalization using the R package limma (Ritchie et al., 2015). In addition, principal component analysis (PCA) was performed using the R packages FactoMineR (Luo et al., 2014) and factoextra (Kassambara and Mundt, 2020). Differentially expressed genes (DEG) are derived from differential analysis by the R package limma and visualized using heatmaps and volcano plots. For identifying the significant DEGs, we set the |log (fold-change) | > 0 and p < 0.05 as the cutoff.

To analyze the expression difference of mitochondrial dysfunction-related genes in all samples, the list of 1,136 genes encoding proteins with strong support of mitochondrial localization, sub-mitochondrial compartment and pathway annotations were downloaded from the MitoCarta3.0 database (http://www.broadinstitute.org/mitocarta) (Rath et al., 2021). From them, 446 mitochondrial dysfunction-related genes were obtained by searching for the keyword “mitochondrial dysfunction” in the GeneCards database (https://www.genecards.org/) (Safran et al., 2010). Then the intersection of these selected mitochondrial dysfunction-related genes (MD genes) and the DEGs from dataset GSE113212 was identified as MDRDEGs.

To systematically explore the underlying biological functions of identified MDRDEGs in HLF, Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed via the R package clusterProfiler (Yu et al., 2012). GO analysis involved three categories: Biological process, cellular component, and molecular function. Gene Set Enrichment Analysis (GSEA) was performed using the R package clusterProfiler.

The R package RCircos (Zhang et al., 2013) is used to map the location of genes in chromosomes. The chromosome data is provided by the R package itself, and the location information of genes on chromosomes is downloaded from the ENSEMBL database (Yates et al., 2020). Pearson correlation between the expression of identified MDRDEGs in elderly individuals was performed.

The protein-protein interaction (PPI) network was constructed based on the Search Tool for the Retrieval of Interacting Genes (STRING) database with the MDRDEGs as input. Medium confidence (interaction score > 0.4) was chosen as a threshold and all positive interactions were included. The results were visualized using Cytoscape software (Shannon et al., 2003). The hub genes were filtered from the whole PPI network using CytoHubba (Chin et al., 2014), a plugin of Cytoscape software and calculated by the Degree method. Twenty genes in the maximum correlation criterion (MCC) were chosen by the CytoHubba plugin and sequentially ordered.

To further understand the regulation of these hub genes, we constructed genes-miRNAs and genes-transcription factors (TFs) regulatory network based on the miRNet database (Chang et al., 2020) predicted the miRNAs and TFs of the hub genes respectively.

To evaluate the infiltration level of immune cells, and the content of 22 different types of immune cells in each sample was calculated based on the LM22 background gene set provided by CIBERSORT (Chen et al., 2018). To further reveal the correlation between the hub genes and the infiltration level of immune cells, a correlation scatter plot for the gene-immune cell pair with a significant correlation with the fitting curve was drawn.

To further verify the selected hub MDRDEGs and related pathways, ligamentum flavum (LF) samples were collected from 30 patients (12 females and 18 males) who underwent posterior lumbar decompression surgery with removal of LF tissue from June 2021 to January 2022. The process of biosample collection followed the ethical approval of the Naval Medical University. For the HLF group, 15 LF specimens were collected from LSS patients with LF hypertrophy, and for the LDH group, 15 specimens were harvested from patients with uncomplicated lumbar disc herniation as control. The thickness of LF was measured on T2-weighted magnetic resonance imaging (MRI) and LF hypertrophy was defined as LF thickness > 4 mm according to the previous study (Zheng et al., 2020). Table 2 summarized the demographic variable of the cases included in different groups. Extensive or partial laminectomy with LF resection was performed in all patients during the operation. The resected ligamentum flavum was rinsed in 4°C physiological saline and then sent for examination immediately.

Ligamentum flavum cells were isolated as described previously (Specchia et al., 2001; Zheng et al., 2020). In brief, LF tissue was washed by PBS 3 times, cut into small pieces measuring around 0.5 mm³, and digested for one hour with 0.2% type I collagenase (Gibco), The digested fragments were then rinsed in DMEM (Gibco), supplemented with 10% FBS (Glpbio, United States), and 100 U/mL penicillin. Cells after the third passage were used for experiments.

LF samples were homogenised to isolate total DNA. real-time quantitative PCR was performed on the mitochondrial gene cytochrome c oxidase subunit 1 (COX1) and the nuclear DNA product (b-actin). The relative mtDNA copy number was calculated by subtracting the cycle of the nuclear gene from the cycle of the mitochondrial gene.

The level of ROS in LF cells was assessed by a C11-BODIPY probe assay kit (Invitrogen) according to the instructions. 1 × 10⁴ LF cells were seeded in 96-well plates and cultured for 30 min with a 2 μM C11-BODIPY probe, and the number of reactive oxygen species (ROS) was measured using a flow cytometer. The malondialdehyde (MDA) and glutathione (GSH) content in tissue homogenates and cell lysis were analyzed by a lipid peroxidation kit (Sigma, MAK085) and Glutathione Assay Kit (Sigma, CS0260) by the standard protocol.

About three cubic meters of ligamentum flavum tissue were homogenized in 800 μL of Trizol. Total RNA was isolated by Trizol reagent (Invitrogen, United States) according to the manufacturer’s instructions, and further reverse transcribed by using the iScript cDNA Synthesis kit (bio-rad). Real-time PCR and analysis were performed as previously described (Xu et al., 2016a). The fold changes of target genes were analyzed by the 2^−ΔΔCT method and 18s was used as an internal control.

The Wilcoxon rank sum test was used to compare the differences between groups. All statistical analyses and visualizations were performed via R software (version 4.0.2). p < 0.05 was considered to be significantly different, and all statistical tests were two-sided.

The gene expression profile of the GS113212 was normalized and shown in Figures 1A–D. A total of 3,742 genes were identified as DEGs, including1,457 downregulated ones and 2,285 upregulated ones (Figures 1E, F).

A total of 43 genes were derived from the intersection of mitochondrial dysfunction-related genes with the DEGs as MDRDEGs. There were 22 downregulated genes and 21 upregulated genes (Figure 2). The most significant enrichment terms for GO were organic acid catabolic process, carboxylic acid catabolic process, mitochondrial respiratory chain complex assembly (biological processes), mitochondrial matrix, mitochondrial inner membrane, mitochondrial protein complex (cellular components), coenzyme binding, oxidoreductase activity, and NADH dehydrogenase activity (molecular functions) (Figure 2; Table 1). KEGG pathway analysis demonstrated that the MDRDEGs were mainly enriched in valine, leucine, and isoleucine degradation, fatty acid metabolism, propanoate metabolism, and fatty acid degradation (Figure 2; Table 1).

The result of GSEA showed that most of the MDRDEGs were involved in cytokine-cytokine receptor interaction, focal adhesion, antigen processing and presentation, lysosome, ECM receptor interaction, ribosome, tight junction, fatty acid metabolism, retinol metabolism, PPAR signal pathway (Figure 3; Table 3).

The expression of MDRDEGs was shown with a heatmap (Figure 4A). The expression differences of all MDRDEGs in young and elderly samples were visualized with grouped box plots (Figure 4B). Compared with the young group, a higher expression of ATPAF2, CLPB, CPOX, LONP1, MRPS34, PREPL, SCO2, SHMT2, TK2, TOMM40 and TXNRD2, and a lower expression of ABCB7, AGK, DBT, IBA57, MFN2, PDE2A, POLG2 and TFAM were observed in the elderly group. The locations of these genes were presented in the genetic landscape (Figure 4C).

Correlations between the expression of MDRDEGs were calculated and presented as a heat map (Figure 4D). Among them, a total of 12 gene pairs met the threshold of statistical significance with an association greater than 0.9 and p < 0.01. They were ATPAF2-HSD17B10, CPOX-MFN2, MPV17-POLG2, SHMT2-ACAT1, SHMT2-DBT, SHMT2-DNAJC30, SLC25A24-NDUFS4, TOMM40-TRMU, TOMM40-LIPT1, TRMU-ACADM, ACADM-BCKDHB, LIPT1-NDUFAF4 (Figures 4E–P). There were 7 pairs with positive correlations and 5 pairs with negative correlations.

To determine the relationships among the mitochondrial dysfunction-related genes, a PPI network was constructed using the STRING database and then visualized using Cytoscape software (Figures 5A, B). According to the degree of the nodes, the top 20 genes including NDUFV1, SURF1, NDUFS4, POLG2, LONP1, ACADM, TK2, HADH, MPV17, NDUFAF4, TFAM, MFN2, NDUFB9, ACAT1, CRAT, HSD17B10, DBT, BCKDHB, SCO2, and TRMU were identified as hub genes (Figure 5C).

The transcription factors (TFs) and miRNAs of these hub genes were predicted using the miRNet database for a comprehensive analysis of their genetic background and regulatory network.

The gene-miRNA network consists of 20 hub genes and miRNAs (Figure 5D). The four interactive hub genes that most miRNA would target were MFN2, TFAM, TRMU, and DBT. Furthermore, the network showed that some of the miRNAs have multiple target hub genes.

A total of 8 pairs of gene-TF interactions among 5 hub genes including NDUFV1, MFN2, SURF1, TFAM, and ACAT1, and 6 TFs including SP1, PPARGC1A, YY1, MYC, STAT1, and PPARG were identified (Figure 5E).

From the enrichment analysis results, we found that there were significant differences in the immune process between the two groups of samples. Therefore, immune infiltration levels of 22 types of immune cells in the two groups of samples were analyzed. Significant differences between groups in CD8⁺ T cells and M0 macrophages were observed (Figures 6A–C).

Correlations between immune cells were analyzed separately in the two groups of samples. The results showed that there were strong correlations between the infiltration levels of immune cells in the elderly group (Figure 6D).

We further analyzed the correlation between hub genes and immune cells, and the results showed that ACADM, BCKDHB, HSD17B10, MFN2, and NDUFAF4 were positively correlated with B cells naïve, resting memory CD4⁺ T cells, regulatory T cells (Tregs) and M1 Macrophages, respectively, while the LONP1 gene was negatively correlated with resting memory CD4⁺ T.

The relative mtDNA copy number was significant lower in patients with HLF. The MDA content and ROS level were significantly increased in the HLF group, whereas the GSH content and SOD activity were markedly decreased in the HLF group (Figure 7).

The expression of the above-mentioned 20 hub genes was verified by qPCR experiment. The relative mRNA expression level of LONP1, TK2, SCO2, TRMU, and MPV17 were significantly higher in HLF samples than in control samples, whereas the expression level of DBT, TFAM, MFN2, POLG2, SURF1, ACADM, NDUFS4 were significantly lower in HLF samples (Figure 8).