Clinical-transcriptomic classification of lumbar disc degeneration enhanced by machine learning | Military Medical Research

Clinical characteristics of the discovery cohort

The discovery cohort comprised 122 NP tissue samples from 108 patients. In all, 94 patients each contributed a single NP sample, whereas the remaining 14 patients each provided 2 distinct samples from different levels, for a total of 28 samples (Additional file 2: Tables S1, S2). The clinical profiles of this cohort are comprehensively described in Table 1. Briefly, the mean age of this cohort was 59.2 years. Among these samples, 59.8% were from male patients, and 77.8% were collected from patients over 50 years of age. The Pfirrmann grade distribution was as follows: 1.6% grade II, 21.3% grade III, 68.9% grade IV, and 8.2% grade V. Regarding disc morphology, 15.6% of samples were classified as normal, 41.8% as bulge, 26.2% as contained herniation, and 16.4% as uncontained herniation. Additionally, 90.1% of samples were located at the L4-5 and L5-S1 levels.

Transcriptome-based molecular classification of 4 distinct LDD subtypes

To identify potential molecular subtypes of degenerative discs from a high-dimensional matrix, SC3 [17] and other clustering solutions were employed. Based on the top 8000 HVGs, 122 NP samples were classified into 4 robust subtypes: 52 (42.6%) into cluster 1 (C1), 32 (26.2%) into cluster 2 (C2), 19 (15.6%) into cluster 3 (C3), and 19 (15.6%) into cluster 4 (C4) (Fig. 1b; Additional file 1: Fig. S1a, b; Additional file 2: Table S1). Surprisingly, of the 28 samples from the 14 patients, both samples from a single patient at different spinal levels were clustered into the same molecular subtype for 16 samples from 8 patients; for the remaining 12 samples from 6 patients, both samples from a single patient exhibited 2 divergent molecular subtypes (Additional file 2: Table S2). These results suggest that discs from different levels of the same LDD patient can present subtype heterogeneity.

To investigate the biological processes of the 4 LDD subtypes, limma [18] was used to identify differentially expressed genes (DEGs) among C1, C2, C3, and C4. Overall, C1, C2, C3, and C4 exhibited 235, 111, 165, and 1,538 uniquely upregulated DEGs and 463, 144, 212, and 1,287 downregulated DEGs, respectively, with cutoffs of log₂ fold change greater than 0.5 and P < 0.05 (Fig. 1c). Gene ontology (GO) analysis revealed that the biological processes of the identified clusters were as follows: C1 was related to axon development and glycosaminoglycan metabolism; C2 was related to chondrocyte development; C3 was related to aerobic metabolism and energy consumption; and C4 was related to inflammation, the immune response and collagen metabolism (Fig. 1d). Subtype-specific upregulated DEGs were visualized via a heatmap (Additional file 1: Fig. S1c). Notably, neuropathic pain mediators (NPY1R and NPY5R) [19] and factors involved in axonal growth (SEMA3A, NGF, and GFRA1) [20, 21] were enriched in C1, whereas signaling molecules involved in bone and cartilage development (DLX5, BMP2, SMAD2, SMAD3, and SMAD5) [22, 23] were specifically upregulated in C2 (Additional file 1: Fig. S1c). Genes encoding the identity of the mitochondrial respiratory chain (MRC) or C-respirasome (COX7A1 and COX7A2) [24] and genes encoding subunits of mitochondrial respiratory chain complex I (NDUFA1, NDUFA2, NDUFA3, NDUFB1, and NDUFC1) were highly abundant in C3 (Additional file 1: Fig. S1c). In C4, genes encoding chemokines (CCL2, CCL3, and CXCL2), proinflammatory cytokines [IL1A, IL1B, and TNF (TNFA)] and the inflammasome (NLRP3) showed significantly high expression (Additional file 1: Fig. S1c). Representative signature genes of each subtype are visualized in a violin diagram (Fig. 1e). Overall, these results indicate that LDD can be categorized into 4 subtypes with different biological processes.

To externally validate the identified subtypes of LDD, we optimized gene selection by evaluating the mean decrease in accuracy and Gini coefficient of DEGs. The optimal 248-gene set comprising subtype-specific DEGs (44 for C1, 15 for C2, 19 for C3, and 170 for C4) showed the minimal error in 10-fold cross-validation (Additional file 1: Fig. S1d; Additional file 3), and was therefore selected to build an RF-based gene classifier, ensuring robust subtype discrimination. The variable importance results showed high similarity among the top 30 genes in the 248-gene classifier (Additional file 1: Fig. S1e, f). In the discovery cohort, the 248-gene classifier allocated 122 NP samples to each of the 4 subtypes with a global precision of 94.3% (115/122), resulting in AUROCs of 0.945, 0.952, 1, and 1 for C1, C2, C3, and C4, respectively, and excellent APs of 0.917, 0.935, 1, and 1 for C1, C2, C3, and C4, respectively (Additional file 1: Fig. S1g). Subsequently, to assess the generalizability of this gene classifier, we applied it to 3 independent microarray datasets from Gene Expression Omnibus (GEO) database representing 62 IVD samples from 62 donors (GSE70362, GSE23130, and GSE15227). We validated the classifier and confirmed that 24 samples from GSE70362 were assigned to C1 (79%, 19/24), C2 (17%, 4/24), and C3 (4%, 1/24), 23 samples from GSE23130 were assigned to C3 (100%, 23/23), and 15 samples from GSE15227 were assigned to C4 (100%, 15/15) (Additional file 1: Fig. S1h). These results supported the predictive potential of the 248-gene classifier and the generalizability of the molecular classification. To further validate the transcriptomic differences revealed by bulk RNA-seq analysis among clusters, we performed quantitative immunohistochemistry (IHC) for selected proteins, including GFRA1, DLX5, COX7A1, and IL-1β. Consistent with the bulk RNA-seq data, semi-quantitative IHC confirmed significantly elevated protein expression of GFRA1 in C1, DLX5 in C2, COX7A1 in C3, and IL-1β in C4 (P < 0.0001) (Fig. 1f). Taken together, these results indicate that 4 molecular subtypes of LDD with different biological processes were identified and that this molecular classification of LDD might be generalizable.

Matrisome characteristics of subtype-specific ECM dysregulation traits

To better understand the specific characteristics of degenerated discs, we evaluated 4 clusters using several gene sets illustrating biological processes related to disc degeneration (Fig. 2a). The results revealed that the core matrisome significantly differed among the clusters, highlighting the existence of different ECM dysregulation patterns. Next, the core ECM structure scores for proteoglycans and collagens revealed that proteoglycans dominated in C1 and that collagens were enriched in C4 (Fig. 2b), further indicating distinct ECM composition and structure among the clusters.

Matrisome was categorized into the core matrisome (collagens, proteoglycans, and ECM glycoproteins) and the associated matrisome (ECM-affiliated proteins, ECM regulators, and secreted factors) [25]. Given that ECM composition and structure govern IVD homeostasis, we speculted that subtype-specific ECM protein aggregation and degradation might contribute to distinct mechanical and chemical microenvironments. To investgate this, we first performed comprehensive profiling of matrisome gene expression across all samples (Additional file 1: Fig. S2a), and then identified subtype-specific differentially expressed matrisome genes (DEMGs) by intersecting subtype DEGs with matrisome genes (Additional file 1: Fig. S2b, c). Subsequently, we mapped the interaction network of upregulated DEMGs (Fig. 2c). The results revealed that the DEMGs of the 6 matrisome subcategories and their potential interactions were distinct among the clusters, reflecting a unique regulatory pattern per subtype. These patterns might contribute to the subtype-specific mechanical and biochemical microenvironments of the IVD and the different stages of disc degeneration.

In C1, DEMGs encoding proteoglycans and collagens were predicted to directly interact with DEMGs encoding other matrisome subcategories (Fig. 2c). Gene set enrichment analysis (GSEA) revealed that proteoglycan biosynthesis-related terms were enriched (Fig. 2d). These results suggest that C1 represents mild degeneration. Additionally, genes encoding ECM core components (ACAN and COL2A1) were highly expressed, and the cell vitality promoter (SERPINA1) [26] and growth factors promoting chondrocyte proliferation (EGF and FGF2) [27] presented similar high expression patterns (Fig. 2c, e). The expression of mechanosensors (PIEZO1 and TRPV4) was also strongly increased, which could enhance ECM synthesis (mainly collagen II) to sustain mechanical loading [28, 29] (Fig. 2e). This collagen II-dominated ECM remodeling in response to mechanical loading was considered to indicate collagenesis, suggesting potential mechanical adaptive repair of LDD.

In C2, the ECM interactome exhibited among DEMGs encoding ECM glycoproteins, secreted factors, and ECM regulators (Fig. 2c). However, DEMGs encoding proteoglycans and collagens were not abundant. These results suggest that ECM degradation might lead to structural disc failure. The GSEA scores also indirectly illustrated the potential for structural repair of discs via bone morphogenesis and chondrocyte differentiation (Fig. 2d). In addition, ADAMTS5 upregulation accelerated proteoglycan degradation, and THBS1 upregulation triggered the TGF-β1/Smad3 signaling pathway and promoted ECM remodeling in response to mechanical stress [30], whereas ASPN upregulation promoted collagen mineralization [31] (Fig. 2e). In particular, TGFB2, BMP2, BMP4, IGF1, and COL10A1 were upregulated to promote bone morphogenesis and chondrocyte hypertrophy [23] (Fig. 2c, e). Overall, the C2 subtype might present with structural failure of the ECM, and endochondral ossification might be activated to repair the disc and stabilize the functional spinal unit.

In C3, the abundance of interacting DEMGs was minimal, with ECM glycoproteins, ECM regulators and secreted factors mainly serving as the ECM interactome, whereas the abundance of collagens and proteoglycans was significantly decreased (Fig. 2c). These results suggest that the ECM structure of C3 becomes dysfunctional. However, the significant enrichment of oxidative energy metabolism-related terms (Fig. 2d) correlated with C-respirasome activation (Fig. 1e), potentially enhancing ECM integrity and mechanostability in cartilage [32]. These results suggest that the C-respirasome was activated to adapt to oxidative stress, because the C-respirasome was more bioenergetically efficient under oxidative conditions [24]. The downregulation of HIF1A (Fig. 2e), a transcription factor that regulates cell behavior and viability in disc degeneration [33], also supported the oxidative microenvironment. Additionally, genes mediating inflammation and immunity (S100A8, S100A9, and S100B) (Fig. 2e), especially S100B, were found to play fundamental roles in the spatial and temporal regulation of chondrogenesis [34, 35]. Unfortunately, a fewer chondrogenic genes were highly expressed in C3 (Fig. 2e). C3-subtype discs might exhibit an adverse oxidative microenvironment that inhibits chondrogenesis, resulting in severe structural damage.

In C4, collagens rather than proteoglycans were significantly enriched and interacted with other modules, forming abundant potential interactions (Fig. 2c). These results suggest that in C4, the ECM was remodeled via fibrogenesis. Inflammation and immune response-related biological processes (Fig. 2d) indicated that LDD patients experienced a severe inflammatory cascade. In addition, fibrosis markers (COL1A1, COL1A2, TNF, IL1B, MMP2, and MMP9), especially COL1A1 and COL1A2, contributed to the fibrotic phenotype (Fig. 2e). ADAMTS4, MMP2, and MMP9 are involved in ECM degradation, a process regulated by IL1B [36]. Notably, TNF and IL1B increase ADAMTS4 expression to degrade aggrecan and upregulate the expression of COL1A1 [36, 37], potentially promoting fibrogenesis. These results suggest that the inflammatory environment might drive fibrogenesis in C4-subtype discs.

Immunostaining confirmed the representative ECM protein expression among the 4 subtypes: C1 showed a significantly higher aggrecan (ACAN, P < 0.0001) and collagen II (COL2, P < 0.01) levels, whereas C4 exhibited elevated collagen I (COL1, P < 0.0001) (Fig. 2f). Taken together, these results indicate that the 4 subtypes exhibit distinct ECM dysregulation patterns, suggesting subtype-specific ECM remodeling in LDD.

Subtype-specific cell subpopulations dominated the ECM phenotype per subtype

Given that bulk RNA-seq can provide an overview of NP biological differences in the process of disc degeneration, it obscures the intricacies of changes within and across cell types and cannot reflect the functions of individual cell subsets. To determine the cell subpopulation composition of each subtype and explore their biological differences, two publicly available single-cell RNA-seq (scRNA-seq) datasets of IVD samples were retrieved from GEO database (GSE160756 and GSE165722) for data integration and analysis (Fig. 3a). A single-cell reference matrix of degenerated NP samples was identified by referencing our published NPC markers [38] and other published immune cell markers [39], which revealed relatively comprehensive cellular populations in the degenerated NPs (Fig. 3a, b; Additional file 1: Fig. S3a, b). Subsequently, BayesPrism [40] was employed to integrate a single-cell reference matrix of degenerated NP and the bulk RNA matrix. The results revealed that Chond2, Chond3, and NP progenitor cells (NPPCs) were highly abundant in C3, C1, and C2, respectively. Notably, C4 had a significantly high macrophage content (Fig. 3c). A relatively higher fraction of macrophages was also found in C4 by the CIBERSORT algorithm with the leukocyte signature matrix (LM22) gene signature [41] (Additional file 1: Fig. S3c). To investigate whether a high cell content is related to the ECM phenotype per subtype, we visualized the aforementioned DEMGs expression landscape of each subtype in all the cell subpopulations and found that Chond1-3 presented a relatively high abundance of ACAN and COL2A1, which contribute to the ECM phenotype in C1. In parallel, TRPV4 and PIEZO1 were upregulated in Chond2, further suggesting that Chond2 enhances ECM synthesis to adapt to biomechanical changes (Additional file 1: Fig. S3d). Notably, Chond3, with high SERPINE1 expression, was the dominant contributor to the antiapoptotic effect and cell viability (Additional file 1: Fig. S3d). Bone morphogenesis-related genes (TGFB2, BMP4, and THBS1) were enriched in NPPCs (Additional file 1: Fig. S3d), suggesting that NPPCs undergo osteogenic differentiation during the progression of disc degeneration in C2. C-respirasome signature molecules COX7A1 and COX7A2 [24] were expressed primarily in Chond1, with downregulation of HIF1A (Additional file 1: Fig. S3d), further indicating that Chond1 is involved in the hypoxic microenvironment dysregulation in C3. Surprisingly, COL1A1 and COL1A2 were especially expressed in stromal cells, which might contribute to the fibrotic remodeling in C4 (Additional file 1: Fig. S3d). Notably, NPPCs also expressed relatively high levels of COL1A1 and COL1A2, suggesting that NPPCs in C4 presented a fibrotic phenotype.

By utilizing Scissor [42], we further revealed a high percentage of macrophages (17.2%) in C4 (Fig. 3d, e). The high expression of IL1B and TNF in macrophages (Fig. 3f) suggested inflammation and an immune response in C4. To elucidate the potential role of macrophages, we performed CellChat analysis, and the results revealed that macrophages were the dominant influencer in the TNF signaling pathway and interacted with NPPCs and Chond3 (Fig. 3g). Moreover, the C-X-C motif chemokine ligand (CXCL) signaling pathway network also revealed interactions between NPPCs and Chond3 (Fig. 3h). These intercellular interactions likely promote fibrogenesis in the NP tissues of C4. Similarly, Scissor-based prediction for C1, C2, and C3 revealed that the dominant cell subpopulation of each subtype contributed to the subtype-specific ECM dysregulation patterns (Additional file 1: Fig. S3e), which was in line with the above results (Fig. 2d). Taken together, each subtype presented distinct and unique cell subpopulations contributing to the subtype-specific ECM dysregulation phenotypes.

Interpretable machine learning-based diagnostic prediction models that stratify molecular subtypes on the basis of clinical features

Diagnosing the molecular subtypes of LDD at the bedside remains challenging because of the potential impairment caused by disc puncture. Thus, we investigated whether a machine learning-based diagnostic prediction model could be developed on the basis of clinical features to predict molecular subtypes in clinical scenarios. First, the LASSO regression analysis was conducted to reduce the aforementioned 23 independent variables to 12 (Additional file 1: Fig. S4a), including age, body mass index (BMI), course of disease, numerical rating scale (NRS) score of low back pain, NRS score of sciatica, neurogenic claudication, straight-leg-raising test, spondylolisthesis, intervertebral disc height (IDH), Pfirrmann grade, disc morphology, and Modic changes. Next, 122 samples were randomly split into a training set (70%) and a testing set (30%) to avoid overfitting. Comparisons of clinical data between the model training and testing cohorts revealed no statistically significant differences (Table 2), indicating comparability. These potential variables were assessed among subtypes (Fig. 4a; Table 3). These 4 subtypes exhibited significant differences in age, course of disease, neurogenic claudication, spondylolisthesis, IDH, Pfirmann grade, and disc morphology (P < 0.05). Specifically, samples in C1 and C2 were predominantly from patients aged 50–70 years. Samples in C3 were mainly from patients older than 70 years. Most samples in C4 presented with uncontained disc herniation. Five different machine learning models, RF, SVM, XGBoost, MLR, and NNet, were evaluated using 10-fold cross-validation on the basis of the selected features. The hyperparameters of each model were optimized through a grid search within the cross-validation framework. Comprehensive analysis of the classified multimodel demonstrated that the RF model was considered the optimal model, because it achieved excellent discrimination with the highest AUROC and AP (Fig. 4b, c).

To identify whether the RF model was clinically useful for this 122-sample cohort, the original 4-class molecular subtype was restructured into a set of binary classification problems, each corresponding to a different class against all others. DCA suggested that this RF model provided significant clinical benefit for identifying the C1–C4 subtypes, with probability thresholds of 5–64%, 0–87%, 9–80%, and 0–78% for each subtype, respectively (Additional file 1: Fig. S4b). To visually explain the selected features in the RF model, SHAP was utilized to illustrate how these features stratify discs into different molecular subtypes. Attributes of the most important features for each subtype in the RF model are shown in Fig. 4d and Additional file 1: Fig. S4c. These results revealed that a high IDH, bulge morphology, short course of disease, and high BMI positively contribute to C1. Interestingly, the high IDH in C1 aligned with the high expression of ACAN and COL2A1 (Fig. 2e). A low IDH positively contributed to C2, followed by a low BMI and Pfirrmann grade V. A long course of disease, advanced age and contained herniation positively contributed to C3. Uncontained disc herniation, which triggers an inflammatory foreign body response [3], was a positive contributor to C4, consistent with the high expression of TNF and IL1B.

Additionally, intersubtype comparisons of the 12 clinical features provided significant information for further understanding the associations between molecular subtypes and clinical features (Table 3). Patients in the C2 and C3 groups were significantly older than those in the C1 and C4 groups. A long course of disease (median: 120 months) was observed for C3-subtype patients, highlighting that the hypoxic microenvironmental dysregulation was not achieved overnight but rather was the result of prolonged disease progression. Patients in the C1 group had excellent IDH retention while those in the C2, C3, and C4 groups presented with IDH loss, which aligned with the low expression of ACAN and COL2A1, supporting ECM collagenesis remodeling in C1. Thus, we speculate that extended conservative treatment may potentially benefit patients in the C1 group. Nearly half of the C2 patients presented with spondylolisthesis, significant IDH loss and Pfirrmann grade V degeneration (Table 3). Given that spondylolisthesis is an indicator of spinal instability [43], IDH recovery and spinal stability reconstruction may be the top priorities of surgical intervention for these patients. Most C3 patients presented neurogenic claudication (Table 3), the typical clinical presentation of spinal stenosis, indicating ischemic injury or mechanical compression of nerve roots [44]. This means that surgical intervention may be necessary for C3 patients rather than biotherapy as an initial approach. The C4 subtype is characterized by uncontained herniation triggering the inflammatory cascade, which is a double-edged sword that promotes herniated material resorption and causes pain [3]. Thus, managing inflammation may be the core challenge of treating patients in the C4 group. Moreover, to illustrate the interpretability of the RF model, a waterfall plot was used to visualize the impact of features in a typical example from each subtype on the model output. The SHAP values for typical patients from the C1, C2, C3, and C4 groups were 0.906, 0.972, 0.976, and 0.996, respectively (Additional file 1: Fig. S4d). Collectively, the attributes of the clinical features could explain the stratification of each patient.

To examine the accuracy of the RF model in stratifying LDD molecular subtypes, 25 NP samples from 25 LDD patients were collected at Army Medical Center of PLA and used as an independent validation cohort (Additional file 2: Table S3). On the basis of the bulk RNA matrix, the 248-gene classifier, the gold standard for LDD molecular classification, assigned 12, 5, 4, and 4 samples to the C1, C2, C3, and C4 subtypes, respectively. The RF model leveraging 12 clinical features stratified 13, 5, 3, and 4 samples to the C1, C2, C3, and C4 subtypes, respectively (Fig. 4e). Confusion matrix and ROC analyses revealed that the RF model identified LDD molecular subtypes with an accuracy of 0.84 and an AUROC of 0.9312 (Fig. 4e, f). The RF model performed robustly in an independent cohort, indicating that the RF model is a clinically applicable tool for LDD patient stratification and diagnosis. Taken together, the RF model-based association between molecular classification and clinical features holds significant potential for guiding the treatment of LDD patients in the clinical setting.

Macrophages dominate inflammation-induced fibrosis in C4-subtype discs

Considering that the macrophage-rich C4 subtype could be predicted using the proposed RF model, ten NP samples were collected at Army Medical Center of PLA and used as another independent cohort. The RF model was used to predict the molecular subtypes of these samples based on the clinical features of the patients (Additional file 4), followed by single-cell extraction. CD235a⁻CD31⁻CD68⁺ macrophages were examined by fluorescence-activated cell sorting, and CD235a⁻CD31⁻CD68⁻CD45⁻ NPCs were sorted for passage culture (Fig. 5a). As expected, the proportion of CD235a^–CD31^–CD68⁺ macrophages in C4 was significantly higher than that in other subtypes (Fig. 5b), which aligns with prior studies indicating a 3-fold increase of CD68⁺ macrophages in uncontained disc herniation [45].

To validate the ECM phenotype of M1 macrophage-NPC interactions, THP-1 (ATCC TIB-202, Manassas, VA)-derived M1 macrophages were indirectly cocultured with P2 CD235a⁻CD31⁻CD68⁻CD45⁻ NPCs from C4 samples (Fig. 5a). The results revealed that TNF-α (10 ng/ml)-treated NPCs upregulated collagen I expression, whereas the expression of aggrecan and collagen II was significantly decreased. In parallel, NPCs cocultured with M1 macrophages showed similar expression patterns of collagen I, collagen II and aggrecan; moreover, these interaction effects were partially reversed by the TNF-α inhibitor etanercept (10 ng/ml, HY-108847, MCE) (Fig. 5c; Additional file 1: Fig. S5a). Collectively, these findings indicate that macrophage-secreted TNF-α contributes to COL1A1 upregulation.

We further investigated the mechanism of TNF-α regulation of COL1A1. Immortalized NP cells (INPCs) (iCELL-0028a, iCell Bioscience Inc., Shanghai) were exposed to varying doses of TNF-α. The results revealed a significant increase in COL1A1 expression at both the mRNA and protein levels with 10 ng/ml TNF-α treatment (P < 0.001; Additional file 1: Fig. S5b-d), which aligns with previous findings [37]. Interestingly, NFKB1 was significantly upregulated in C4 (Fig. 5d), whereas other fibrosis-related transcription factors (JUN, SP1, TGFB1, and SMAD3) did not show significant changes in C4 (Additional file 1: Fig. S5e). Furthermore, a positive correlation was confirmed between COL1A1 and NFKB1 (r = 0.61, P < 0.001; Fig. 5e) but not between COL1A1 and JUN, SP1, TGFB1, or SMAD3 (Additional file 1: Fig. S5f). These results suggested that targeting inflammation might have potential benefits for ameliorating fibrosis of C4. We hypothesized that TNF-α might induce NF-κB1 (p50) activation to increase COL1A1 expression. Subsequent treatment of INPCs with 10 ng/ml TNF-α and 10 μmol/L JSH-23 (a selective inhibitor of NF-κB signaling) significantly reduced COL1A1 expression in the JSH-23-pretreated group (10 μmol/L) (P < 0.01; Fig. 5f), demonstrating the impact of NF-κB activation on COL1A1 expression.

In unstimulated cells, NF-κB remains inactive when bound to inhibitors of κB (IκB, IκBα, and IκBβ) in the cytoplasm. Treatment of INPCs with TNF-α resulted in rapid degradation of IκBα within 30 min and a concomitant increase in cytosolic p-IκBα and nuclear NF-κB1 (p50) levels (Fig. 5g, h; Additional file 1: Fig. S5g, h). Immunostaining confirmed NF-κB1 (p50) sequestration in the nucleus after 30 min of TNF-α treatment (Fig. 5i). Thus, TNF-α induces IκBα degradation in INPCs and facilitates rapid translocation of NF-κB1 (p50) to the nucleus, a hallmark of NF-κB signaling activation. Subsequently, luciferase reporter vectors containing wild-type and mutant (MUT) NF-κB1 (p50) binding sequences of COL1A1 were constructed. Cotransfection of the COL1A1 promoter fragment with the NF-κB1 (p50)-expressing vector significantly increased the luciferase activity compared with that of the control (the COL1A1 promoter fragment cotransfected with the empty vector) (P < 0.001; Fig. 5j). Conversely, cotransfection of the NF-κB1 (p50)-expressing vector with the COL1A1-MUT construct led to a significant decrease in luciferase activity over the former (P < 0.001) but did not affect luciferase activity compared with the control (Fig. 5j). These findings suggest that NF-κB1 (p50) positively modulates COL1A1 promoter activity (Fig. 5k). Overall, TNF-α-induced activation of NF-κB1 (p50) leads to the upregulation of COL1A1 expression.