Introduction
Intervertebral disc degeneration (IVDD), a common aspect of the natural aging process, often leads to lower back pain (LBP), which is a major cause of disability worldwide and has a significant impact on public health.1 However, lumbar fusion surgery, which is currently used mainly to relieve end-stage symptoms, comes with several complications, including a high risk of nerve injury, increased radiation exposure, and accelerated degeneration of adjacent segments.2 Importantly, this procedure fails to address the underlying pathobiology of IVDD. The IVD, a non-self-renewing and avascular structure, consists of a gel-like nucleus pulposus (NP) encased within a fibrocartilaginous annulus fibrosus (AF) and bordered vertically by cartilaginous endplates (CEP), with the degeneration of the NP being the primary factor driving the onset and progression of IVDD.3 Therefore, identifying and mitigating the degeneration of nucleus pulposus cell (NPC) within the IVD during the early stages of IVDD is a crucial therapeutic strategy for addressing the condition.
NPC are generally described as a homogeneous population derived from embryonic notochord cells (NCs). However, data from multi-level sequencing techniques and lineage-tracing animal models suggest that NP cells may actually constitute a heterogeneous population.4 Significant changes in the composition and fate of cell types within the NP during degeneration have drawn increasing research attention to the role of NPC heterogeneity in the progression of IVDD. Single-cell transcriptome analysis provides an unbiased quantitative expression analysis at the single-cell level, which facilitates the identification of various cell types and the extent of immune cell infiltration, thereby enhancing our understanding of cell heterogeneity and microenvironment interactions during IVDD.5 Zhang et al used single-cell RNA sequencing (scRNA-seq) to analyze the transcriptional landscape of the NP and AF in healthy human IVD, identifying specific markers and their roles in maintaining IVD homeostasis.6 Gan et al classified IVD chondrocyte-like cells into three subpopulations based on their role in extracellular matrix (ECM) homeostasis and identified a new progenitor cell subpopulation characterized by specific genes and progenitor cell potential.7 Therefore, analyzing the development of different NPC subpopulation in IVDD using high-resolution single-cell transcriptomics helps identify potential targeted pathways to mitigate the progression of IVDD.8,9
Epidermal growth factor receptor (EGFR) is a tyrosine kinase receptor widely distributed across various tissues and organs in the human body.10 It recruits adaptor proteins and tyrosine kinase substrates through autophosphorylation of receptor tyrosine residues, thereby triggering a cascade of downstream signaling events that participate in the onset and progression of various degenerative diseases.10 Previous studies have demonstrated that molecules such as mitogen-inducible gene 6 (MIG6) can directly interact with the EGFR family, influencing both collagen synthesis and degradation, while EGFR signaling has been linked to the inhibition of cartilage differentiation.11,12 Caroline et al discovered that EGFR signaling involving transforming growth factor α (TGF-α) can inhibit both chondrogenesis and myogenic differentiation in limb mesenchyme.13,14 In addition, excessive activation of EGFR signaling can stimulate the activity of matrix metalloproteinases (MMPs), thereby promoting catabolic processes, which may be associated with abnormal degradation and remodeling of the ECM.11,15,16 Gefitinib, as an EGFR inhibitor, has been found to suppress the expression of MMP9, MMP13, and MMP14 during endochondral ossification, leading to an increase in the amount of collagen fibrils and a significant thickening of the epiphyseal growth plate.17 Additionally, Pan et al utilized an injectable thermosensitive hydrogel to achieve the sustained release of gefitinib, thereby delaying the progression of IVDD by promoting ECM synthesis, inhibiting ECM degradation, and activating autophagy.18
Machine learning, with its powerful classification capabilities, has recently been extensively used to learn the representation of high-dimensional features derived from high-throughput data. Additionally, when combined with high-throughput sequencing analysis, machine learning has been widely used to identify new diagnostic biomarkers. Therefore, we analyzed single-cell transcriptomic changes and integrated three machine learning algorithms—least absolute shrinkage and selection operator (Lasso), extreme gradient boosting (XGBoost), and Random Forest (RF)—to identify key EGFR-related regulatory factors in the progression of IVDD. We then created NP-specific EGFR inactivation mice using the Cre-loxP system and established an NP needle puncture model to validate the regulation of Janus Kinase 1 (JAK1) expression by EGFR during IVDD progression, proposing a new target for mitigating IVDD progression through the modulation of nucleus pulposus cell degeneration.
Materials and Methods
Single-Cell RNA Statistical Analysis
IVD gene expression datasets (GSE165722) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Seurat package (version: 4.0.3, https://satijalab.org/seurat/) was used for cell normalization and regression based on the expression table according to the UMI counts of each sample and percent of mitochondria rate to obtain the scaled data. Principal Component Analysis (PCA) was constructed based on the scaled data with the top 2000 highly variable genes, and the top 20 principals were used for TSNE construction. Utilizing graph-based cluster method, we acquired the unsupervised cell cluster result based on the top 20 principal-components of PCA. Then, after integrating the data, we applied the Harmony algorithm to correct for batch effects. Briefly, the merged data were first log-normalized, and the top 2,000 highly variable genes were selected. The effects of mitochondrial gene proportion and library size were regressed out. Harmony was run using sample origin as the grouping variable, and the effectiveness of batch correction was evaluated by comparing the distribution of cells in the low-dimensional embedding space before and after correction. Gene markers for each group were identified using the Wilcoxon rank-sum test in the FindMarkers function with Bonferroni correction for multiple testing.
Gene set variation analysis (GSVA): GSVA is a nonparametric and unsupervised method for evaluating enrichment of transcriptome gene set. By comprehensively scoring the gene set of interest, GSVA converts gene-level changes into pathway-level changes, and then judges the biological function of the sample. This study downloaded gene sets from the molecular signatures database, and the GSVA algorithm was used to comprehensively score each gene set, so as to evaluate potential biological functional changes between different NPC clusters.
Pseudotime analysis: The pseudotime analysis was performed using the Monocle2 R package (v2.14.0) to reveal cell state transitions in NPC clusters. The data were normalized by the estimateSizeFactors and estimateDispersions functions with the default parameters. Dimensional reduction and cell ordering were performed using the DDRTree method and the orderCells function. Genes exhibiting significant dynamic changes along Pseudotime were identified using a likelihood ratio test implemented in the differentialGeneTest function, with multiple testing correction via the Bonferroni method. Lastly, the plot_cell_trajectory, plot_genes_in_pseudotime, plot_pseudotime_heatmap and ggplot function was run for visualization.
RNA-Seq Analysis
RNA-seq datasets (GSE70362) were downloaded from the GEO database. Single-gene gene set enrichment analysis (GSEA) enrichment analysis was performed based on EGFR expression levels, with pathway identification using the Kolmogorov–Smirnov test implemented in the GSEA function and multiple testing correction via the Benjamini-Hochberg method. Then, we utilized the weighted gene co-expression network analysis (WGCNA) software to identify gene modules in the combined cohort that were significantly associated with severe degenerative disc (SDD). This enabled us to determine which genes were highly correlated with SDD. To determine the optimal value for the soft threshold in the module analysis, we examined the scale independence and average connectivity of the considered modules using various weighting parameters. Further clustering of genes was performed, and the different branches of this clustering represent different gene modules. Within each module, there is a high degree of gene co-expression, while between different modules, the degree of gene co-expression is low. Modules with higher correlation indicate a set of genes that are closely related to the target disease. Consequently, we were able to pinpoint the genes with a strong association with SDD.
Feature Selection of Machine Learning
To identify genes significantly associated with IVDD, we employed multiple machine learning algorithms as previously described. We utilized the randomForest package for robust ensemble-based prediction, capitalizing on its inherent resistance to overfitting. The algorithm was configured to randomly select 2 features at each split and build 1000 decision trees, ultimately quantifying gene importance through the mean decrease in Gini. Further, we applied Lasso regression, a classical method effective for both categorical and continuous variables with strong feature selection capabilities in high-dimensional data. Lasso enhances predictive accuracy and model interpretability via variable selection and regularization, thereby mitigating overfitting. We set the model family parameter to binomial for binary classification, specified the mixing parameter alpha as 1, and applied 10-fold cross-validation for model optimization, resulting in a final lambda value of 0.02966407. Finally, to prioritize features while addressing class imbalance, we implemented XGBoost. This method calculates feature importance scores based on each feature’s contribution to performance gain (derived from metric improvement after splits). Briefly, we set the learning rate to 0.3 and used repeated cross-validation to optimize the model, with feature importance calculated using the varImp function. Collectively, these approaches ensure robust identification of genes critically influencing IVDD prediction accuracy.
Human NP Tissue Collection and Processing
The Ethics Committee of Union Hospital in Tongji Medical College approved the collection of human NP tissues and the use of publicly available human data ([2025]0423). Human NP tissues were obtained from patients undergoing spinal surgery for conditions such as spinal fractures, adolescent idiopathic scoliosis, or lumbar disc herniation at our hospital. Briefly, human NP tissues were grossly isolated from fresh intervertebral disc samples with careful exclusion of the AF and surrounding tissues. The specimens were then fixed in 4% paraformaldehyde and paraffin-embedded for sectioning.
Rat Primary NPC Isolation and Treatment
All animal work performed in this report was approved by the Animal Experimentation Committee of Huazhong University of Science and Technology (IACUC number: 4085). All animal experiments were conducted in accordance with the national standard of Laboratory animal-Guideline for ethical review of animal welfare (GB/T 35892–2018). Sprague-Dawley rats (weighing 200–250g) were obtained from the Experimental Animal Center of Tongji Medical College, Huazhong University of Science and Technology (Wuhan, China). Following euthanasia via intraperitoneal injection of a lethal dose of pentobarbital, primary NPC wfreere isolated from the caudal intervertebral discs, as previously described.19 Briefly, NP tissues were excised, minced into small fragments, and enzymatically digested in 0.2% (w/v) type II collagenase (Sigma-Aldrich, St. Louis, MO, USA) at 37°C for 15 minutes. After digestion, the cell suspension was centrifuged, and the resulting pellet was resuspended and seeded into culture flasks containing DMEM/F-12 medium (Gibco, Grand Island, NY, USA) supplemented with 10% fetal bovine serum (Gibco, Grand Island, NY, USA) and 1% penicillin-streptomycin (Sigma-Aldrich, St. Louis, MO, USA). Cells were maintained in a humidified incubator at 37°C, and the culture medium was refreshed every three days. To evaluate the effect of EGFR signaling inhibition on NPC under oxidative stress, when the cell confluence reached 70%–80%, the cells were pretreated with Gefitinib (10μM, MCE, New Jersey, USA) for 24 hours, followed by treatment with tert-butyl hydroperoxide (TBHP; 100μM, Sigma-Aldrich, St. Louis, MO, USA) for 6 hours.
Western Blot Analysis
Following treatment, NPC were rinsed twice with phosphate-buffered saline (PBS) and subsequently lysed in RIPA lysis buffer (Beyotime, Shanghai, China) supplemented with protease and phosphatase inhibitors. The resulting lysates were subjected to 7.5% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE), followed by transfer onto polyvinylidene difluoride (PVDF) membranes (EMD Millipore, Billerica, MA, USA). Membranes were blocked with NcmBlot blocking buffer (NCM Biotech, Shanghai, China) for 10 minutes at room temperature and incubated overnight at 4 °C with primary antibodies targeting JAK1 (1:2000, A18323, Abclonal) and GAPDH (1:3000, 4970, Cell Signaling Technology) as a loading control. After washing with TBST buffer, membranes were incubated with species-specific horseradish peroxidase (HRP)-conjugated secondary antibodies for 1 hour at room temperature. Immunoreactive bands were visualized using an enhanced chemiluminescence (ECL) detection system (Biosharp, Hefei, China). Quantitative analysis was performed using ImageJ v1.52 software (NIH, Bethesda, MD, USA) by measuring the gray value of each band normalized to GAPDH.
In vivo Animal Experiments
In accordance with the standards for animal housing, mice were group housed at 23° to 25°C with a 12-hour light/dark cycle and allowed free access to water and standard laboratory pellets. To achieve NP-specific deletion of EGFR, we utilized Col2-Cre mice, in which Cre recombinase expression is driven by regulatory elements of the type II collagen gene.20–22 Col2-Cre mice were bred with EGFRflox/flox mice to generate Col2-Cre EGFRflox/+ mice. Then, these Col2-Cre EGFRflox/+ mice were bred with EGFRflox/flox mice to finally obtain Col2-Cre EGFRflox/flox (Egfr CKOCol2) mice (Figures S1 and S2).18 All mouse lines were purchased from the Cyagen Biosciences and Hubei Biont Biological Technology Co., Ltd.
To induce IVDD, male mice 3 months of age were subjected to NP acupuncture. Briefly, mice were anesthetized with 40 mg/kg pentobarbital. The coccygeal 6/7 and 7/8 IVDs were localized using X-ray imaging, followed by percutaneous puncture to a depth of 1–1.5 mm using 31-gauge needles (approximately 0.25 mm in diameter), inserted parallel to the endplates. X-ray imaging was used to confirm that the needle tip reached the center of the IVD. The mice were randomly separated into three groups, WT, WT+NP acupuncture, and Egfr CKOCol2+NP acupuncture group. The mice were sacrificed 4 weeks after the acupuncture, and the corresponding discs were fixed with 4% paraformaldehyde for 48 h, decalcified with 10% ethylenediaminetetraacetic acid, and embedded in paraffin.
Histology
For histology, 6-μm-thick sections were deparaffinized, rehydrated, and stained with hematoxylin and eosin (H&E) and safranin-O/fast green (SO/FG).23 For immunohistochemical (IHC) staining, the 6-μm-thick sections were incubated with rabbit anti-JAK1 (1:200, A18323, Abclonal) at 4 °C overnight, followed by incubation with horseradish peroxidase-conjugated secondary antibody at room temperature for 1 h and diaminobenzidine DAB color development, and stained with hematoxylin. The integrated optical density (IOD) from the IHC images was quantified using ImageJ, following established procedures as previously outlined.
Statistical Analysis
Receiver operating characteristic curve (ROC) analyses were performed on key genes using the plot.roc function in the R package pROC. The ROC curve was drawn to evaluate the diagnostic value of identified candidate biomarkers and to further screen key genes from candidates. The data of GSE244889 and GSE186542 were used for verification.
To ensure objectivity and reduce potential bias, the allocation of animals or samples, along with data collection and analysis for both in vivo and ex vivo experiments, was carried out under blinded conditions. The statistical analyses of the present study were based on R software (version 3.6.3) and GraphPad Prism (version 9.0; GraphPad Software, San Diego, CA, USA). Differences between groups were analyzed using the one-way analysis of variance, followed by the least significant difference test. P < 0.05 was regarded statistically significant.
Results
scRNA-Seq Atlas of NPC
The overall study flow is illustrated in Figure 1. Using the Pfirrmann grading system, we selected six samples with a Pfirrmann grade of II or higher, which previous researchers have defined as SDD.24 As shown in Figure S3, distinct separation of different samples was observed in PCA space prior to correction (Figure S3A). After Harmony correction, the samples were well mixed in the Harmony embedding space (Figure S3B), indicating that batch effects were effectively removed while preserving biologically relevant variation. Then, through unsupervised graph-based clustering of the integrated dataset, comprising 31,945 cells from the six samples, we identified 18 distinct clusters based on nearest neighbor approximations. NPC were sorted based on sex-determining region Y-box 9 (SOX9) and aggrecan (ACAN) expression for further heterogeneity analysis (Figure 2A and B).24 Consistent with highly expressed genes and findings from published single-cell studies, we categorized NPC into six subpopulations: each characterized by distinct gene expression profiles (Figure 2C and D; Table S1).25
|
Figure 2 scRNA-seq atlas of 6 SDD NP tissues. (A) t-SNE map of NP tissues. (B) t-SNE image of the NPC marker genes: ACAN and SOX9. (C) t-SNE map of the following 6 subpopulations: Eff-NPC, HTC-1/2, Fib-NPC, Reg-NPC, and Hom-NPC. (D) Heatmap of marker genes of the 6 subpopulations. (E) GSVA analysis of biological process.
|
To clarify the functions of each NPC subpopulations, we conducted GSVA analysis on key functions during the IVDD process (Figure 2E). The results showed that fibrous NPC (Fib-NPC) had higher enrichment scores in regulating collagen metabolism and ECM regulation compared to other subpopulations (Figure 2E). Given the close association between apoptosis and autophagy with NP degeneration in IVDD, we also performed enrichment scoring for these processes. The GSVA results indicated that Fib-NPC had significantly higher scores in regulating of apoptotic and autophagy signaling compared to all other subpopulations (Figure 2E).
Pseudo-Time Analysis of NPC
Previous analyses of expression levels suggested a continuum among NPC subpopulations. To explore this, we conducted a pseudotime analysis. The results showed that homeostatic NPC (Hom-NPC) and hypertrophic chondrocyte-like NPC 1(HTC-1) predominantly occupies the early part of the trajectory, while Fib-NPC is mainly distributed in the late stage (Figure 3A). While effector NPC (Eff-NPC) and regulatory NPC (Reg-NPC) present widespread distribution across all stages, suggesting a potential shared fate among these subpopulations (Figure 3A). The pseudotime cell density plot indicates that Hom-NPC and Fib-NPC are at opposite ends, with Eff-NPC closer to Fib-NPC (Figure 3B). Thus, according to the scRNA-seq pseudotime analysis and functional analysis results, the Fib-NPC was characterized as a late-stage subpopulation of NPC and closely related to IVDD. In addition, Eff-NPC may regulate or promote the emergence of Fib-NPC.
 |
Figure 3 Pseudo-time analysis of NPC. (A) cell trajectory plot of NPC. (B) Density plots of NPC along the pseudo-time axis. (C) Expression dynamics of the following genes plotted across pseudo-time: EGFR, FBLN1, COL1A1, COL3A1, IGFBP6, MMP2, PRSS23, COL6A1. (D) Heatmap of significantly gene expression dynamics across pseudotime (P<0.05).
|
Gene expression trends along the temporal trajectory reveal an overall increase in COL1A1, COL3A1, and COL6A1 expression (Figure 3C). Additionally, we plotted the expression of Fib-NPC subpopulation marker genes mentioned in related studies, which showed an increasing trend for FBLN1, IGFBP6, PRSS23, and MMP2 (Figure 3C). Previous studies have shown that EGFR signaling is linked to NP degeneration. Therefore, we examined changes in EGFR expression within Fib-NPC, alongside the previously mentioned Fib-NPC subpopulation-related genes. The results show that EGFR has an increasing trend in the pseudotime trajectory (Figure 3C). To further analyze the impact of pseudotime on gene expression, we used natural spline functions to model this relationship, capturing nonlinear expression patterns. We identified the 2,000 genes with the most significant changes and presented the expression of certain key genes (Figure 3D; Table S2). The results demonstrated a significant increase in EGFR expression, while COL1A1, COL3A1, COL6A1, and other Fib-NPC subpopulation marker genes significantly increased (P<0.05) (Figure 3D).
Single-Gene GSEA and WGCNA of RNA-Seq
The GSE70362 dataset comprises 22 NP tissue samples. Based on related studies, samples with Pfirrmann grades I and II are classified as mild degenerative disc (MDD), while those with grades above II are defined as SDD. Herein, we conducted GSEA on EGFR to explore potential KEGG pathways involved in IVDD progression and 121 pathways are enriched (P.adjust < 0.05; Table S3). The pathways that exhibited the most notable upregulation include Malaria, Autoimmune thyroid disease and Hematopoietic cell lineage. Conversely, the pathways that showed the most significant downregulation consist of Ribosome biogenesis, Spliceosome and Aminoacyl-Trna biosynthesis (Figure 4A).
 |
Figure 4 Single gene GSEA and WGCNA of RNA-seq. (A) Bar graph showing the NES of significantly enriched pathways (significance threshold: P.adjust < 0.05; ranking criteria: |NES|; Table S3). (B) Determination of the optimal soft threshold power. (C) Genes with similar expression patterns were clustered, with different colors representing different gene clusters, and gray modules indicating genes not assigned to any clusters. (D) Heat map of module-trait correlations. (E) Heat map of module-module correlations.
|
We conducted WGCNA on NP samples for identifying gene module associated with SDD, which involved the construction of co-expression networks and the identification of co-expression modules. For the remaining 22 samples, the top 50% of genes with the highest expression variance were selected for the subsequent WGCNA. To determine the appropriate soft threshold power for WGCNA, assessments of scale independence and mean connectivity were conducted. A power value of 14 was selected based on a correlation coefficient greater than 0.9, which ensures optimal scalability of the network (Figure 4B). The was constructed with the soft threshold power of 14 and established the minimum number of genes required for each module in the gene network to be 30 (Figure 4C). Modules that exhibit a strong association with clinical features are frequently discovered to possess significant and specific biological implications. Therefore, our next research will focus on the genes contained in the lightcyan, cyan, green, and darkred modules (Figure 4D; Table S4). In addition, cyan, lightcyan, and green modules have a positive correlation with each other (Figure 4E). Darkred module has a negative correlation with cyan module, but a weak positive correlation with lightcyan and green modules (Figure 4E).
Identification of Candidate Key Genes via Multi-Machine Learning Methods
To investigate key regulatory genes associated with EGFR signaling in the degeneration of NPC, we selected genes from the 60 KEGG pathways most significantly correlated with EGFR expression, as well as the 2000 genes showing the most notable changes in pseudotime trajectories (Table S5). To further identify significant genes, the aforementioned genes were intersected with those in the four modules showing significant correlation in WGCNA, resulting in the identification of five hub genes (Figure 5A).
 |
Figure 5 Selection of key genes for IVDD using multi-machine learning algorithms. (A) The Venn diagram shows the intersection of significant genes from pseudotime analysis, single gene GSEA analysis, and WGCNA. (B) Four genes with importance greater than 1.0 were selected based on random forest algorithm. (C) Features screening in the LASSO regression model based on 10-fold cross validation, three genes corresponding to the lowest point were considered the most superior biomarkers. (D) Three genes with importance scores greater than 50 were selected via the XGBoost algorithm. (E) The intersection of RF, lasso and XGBoost was performed and visualized with a Venn plot, resulting in JAK1 and IDI1 being selected for further ROC evaluation.
|
We employed Lasso, XGBoost, and RF machine learning strategies to screen candidate key genes for further analysis (Table S6). RF in combination with feature selection was used to determine the connection between the error rate, the number of classification trees, and ranking genes in descending order based on their importance scores (Figure 5B). The Lasso algorithm successfully identified three feature variables when the lambda value was minimized (Figure 5C). The XGBoost algorithm selected three genes with importance scores greater than 50 (Figure 5D). Ultimately, intersection on the results of three machine learning algorithms and visualized through Venn diagram. Two candidate key genes (JAK1and IDI1) were obtained for further ROC assessment (Figure 5E).
JAK1 Association with Nucleus Pulposus Degeneration and Regulation by EGFR
To evaluate the diagnostic potential of JAK1 and IDI1 for predicting IVDD, we constructed ROC curves based on individual characteristics and gene expression. The Area Under Curve (AUC) was determined for both the training and validation datasets (GSE244889 and GSE186542). In the training dataset, the AUC for JAK1 and IDI1 exceeded 70% (Figure 6A). But in the validation dataset, only JAK1 (AUC=83.3%) exhibited favorable efficiency (Figure 6A). We further collected human NP tissues and found that JAK1 expression was significantly elevated in degenerated NP tissues (Figure 6B and C). Then, the JAK1 was selected as the core gene in this study for subsequent experiments.
 |
Figure 6 ROC analysis and validation of JAK1 expression. (A) The ROC curves of JAK1 and IDI1 in the training dataset (AUC: 72.9% and 85.4%) and integrated validation dataset (AUC: 83.3% and 50.0%). (B, C) Representative IHC images of JAK1 in NP and corresponding statistical analysis (n = 7/group). (D, E) Representative Western blot plots of JAK1 and corresponding statistical analysis (n = 3/group). (F, G) Representative images of H&E, SO/FG, and IHC staining in IVD tissue, with corresponding quantification of JAK1 expression (n = 7/group). (**P < 0.01, ****P < 0.0001).
|
To preliminarily investigate whether EGFR regulates JAK1 expression in NPC, we inhibited EGFR signaling using Gefitinib and established an oxidative stress model via TBHP treatment. Western blot demonstrated significantly elevated JAK1 expression in TBHP-treated NPC, whereas Gefitinib pretreatment abolished this upregulation (Figure 6D and E). We then performed H&E staining to assess morphology, fibrous tissue, and NP boundaries, while SO/FG staining was employed to detect proteoglycans and collagen within the IVD (Figure 6F). To validate the diagnostic role of JAK1 and its regulation by EGFR in vivo, we developed an IVDD model using Egfr CKOCol2 mice. Immunohistochemistry was then used to compare the expression of JAK1 in NP tissues. The results showed an increase in JAK1 expression in the IVDD model mice, while this increase was significantly reduced in Egfr CKOCol2 mice (Figure 6F and G).
Discussion
IVDD is a complex process characterized by the degeneration of NPC, activation of matrix-degrading enzymes, imbalances in catabolic and anabolic activities, and heightened inflammatory responses. The NP is a vital component of the IVD, providing an avascular and hypoxic environment that supports chondrocyte-like cells within an ECM rich in proteoglycans and type II collagen. This environment is crucial for regulating disc function and deformation. In this study, we began our analysis at the single-cell transcriptional level, examining gene expression changes during NPC degeneration. We explored the interactions among different NPC subpopulations and integrated bulk RNA sequencing data with various machine learning methods. Through this approach, we identified key molecules that could inform potential treatment strategies to delay NPC degeneration and mitigate IVDD.
Previous studies have confirmed that changes in NPC significantly contribute to IVDD, and scRNA-seq offers the potential to analyze NPC subpopulation. Fib-NPC, predominantly found in degenerative NP tissues, are classified as late-stage NPC and are positively correlated with the severity of IVDD.25 Using ACAN and SOX9 markers, we identified NPC and divided them into six distinct subpopulations. The Fib-NPC subpopulation was the most prevalent, consistent with our sample selection results. Previous research has shown that Fib-NPC exacerbate local inflammation in IVDD by increasing the secretion of inflammatory cytokines. This promotes macrophage infiltration and M1 polarization, ultimately leading to highly disorganized lamellar fibers with rim lesions, delamination, and radial fissures.24 After categorizing NPC subpopulations based on characteristic genes, we analyzed their cellular functions using GSVA. Given their close association with NP fibrosis, we focused on the functional status of Fib-NPC in relation to collagen and ECM.26 Our results indicate that the Fib-NPC subpopulation has a higher functional enrichment score for the regulation of collagen metabolic and cell matrix adhesion compared to other subpopulations.
The enrichment of Fib-NPC in apoptosis and autophagy related signaling pathways is generally higher than other subpopulations. But interestingly, the HTC-2 subpopulation has similar enrichment scores to Fib-NPC in terms of positive regulation of autophagy and positive regulation of macroautophagy, and even higher enrichment scores in mitochondrion related autophagy signaling. In addition, HTC-2 is similar to Fib-NPC in terms of ECM assembly and organization functions. These findings suggest that Fib-NPC and HTC-2 NPC share similar functional states. Previous studies have shown that during cartilage development, chondrocytes undergo terminal differentiation and become hypertrophic, a process known to contribute to osteoarthritis (OA).26 Similarly, during the progression of IVDD, NP tissue undergoes hypertrophic differentiation, which may involve apoptosis and autophagy-related cellular dysfunction events.27 Therefore, HTC-2 NPC may represent a cellular functional state that precedes Fib-NPC in the process of NPC degeneration.
To trace cellular differentiation and developmental states, and to identify key cell types and gene expression changes during disease progression, we conducted a pseudotime analysis on the various NPC subpopulations. The results indicate that Hom-NPCs are located at the beginning of the trajectory, while Fib-NPC are positioned at the end. Consistent with the GSVA results, HTC-1 was primarily located in the early stage of the trajectory, whereas HTC-2 was positioned closer to the Fib-NPC subpopulation. This finding further suggests that HTC-2 may represent a terminally hypertrophic subpopulation that precedes Fib-NPC. Additionally, the results of gene expression along the pseudotime trajectory indicate that the expression of other Fib-NPC subpopulation marker genes (FBLN1, IGFBP6, MMP2, PRSS23) significantly increases in the later stages of the trajectory. Additionally, the collagen components COL1A1, COL3A1, and COL6A1 show a pronounced upward trend in the trajectory. Therefore, our analysis indicates that Fib-NPC a late-stage degenerative NPC subpopulation closely associated with IVDD. Additionally, identifying key regulatory genes involved in NPC degeneration by screening for genes with significant expression changes along this pseudotime trajectory might be an effective approach.24
Interestingly, although Eff-NPC are distributed throughout the entire trajectory, their peak distribution appears to be closer to Fib-NPC. Meanwhile, GSVA analysis also shows that Eff-NPC have high enrichment scores in both collagen and ECM metabolism, as well as in apoptosis and autophagy signaling. Previous studies have further performed subpopulation analysis of Eff-NPC, revealing that a subset of Eff-NPC is primarily present in normal and mildly degenerated IVD, where they are involved in maintaining cartilage-like properties and synthesizing glycosaminoglycans.8,25 In contrast, another subset of Eff-NPC predominates in severely degenerated discs, exhibiting functional shifts toward matrix mineralization and the regulation of senescence-related signaling pathways, which are closely associated with the pathological features of IVDD.8,25 Considering that the scRNA-seq data included in this study were all derived from patients with IVDD, it is plausible that Eff-NPC may also participate in extracellular matrix remodeling and apoptotic signaling during IVDD progression.
Single-gene GSEA enrichment analysis helps identify the role of target genes in biological pathways and assess changes in cellular functions under different phenotypic states.28 Using data from GSE70362, we performed GSEA enrichment analysis to identify pathways significantly associated with EGFR. Consistent with previous studies, our results indicate that the PI3K/Akt signaling pathway is associated with the progression of IVDD, and that EGFR and PI3K/Akt signaling are subject to cross-regulatory interactions.29,30 In normal IVD, moderate activation of the PI3K/Akt/mTOR pathway promotes protective autophagy, facilitating the clearance of damaged organelles and maintenance of cellular homeostasis.31 However, sustained activation may inhibit autophagy and accelerate cellular senescence.31 Moreover, PI3K/Akt signaling can suppress the expression of matrix-degrading enzymes such as MMPs and a disintegrin and metalloproteinase with thrombospondin motifs-4/5 (ADAMTS-4/5), thereby reducing collagen and proteoglycan degradation. Under inflammatory conditions, persistent phosphorylation of Akt can activate the IKK/NF-κB pathway, leading to the release of pro-inflammatory cytokines and the formation of a positive feedback loop that exacerbates IVDD progression.31,32
The JAK/STAT signaling pathway is involved in many critical biological processes, such as immune regulation, inflammatory responses, and cell apoptosis.33 Previous studies have shown that cytokines like TNF-α, IL-1β, and IL-6 promote autophagy, senescence, or apoptosis in IVD by triggering a series of inflammatory responses. The JAK/STAT pathway, as a downstream pathway of IL-6, has been reported to have a positive feedback loop that frequently and significantly promotes IL-6 autocrine secretion in various cancer cell lines and clinical samples.34 Wu et al discovered that treatment with various concentrations of Resveratrol reduced the expression levels of IL-6 and the phosphorylation levels of JAK1/STAT3 in NPCs, suggesting that Resveratrol interrupts the IL-6/JAK1/STAT3 positive feedback loop.35 Additionally, activation of the JAK/STAT3 pathway can directly increase the expression levels of cyclooxygenase-2 (COX-2) and MMP13, leading to IVDD, without the mediation of classical cytokines IL-1β or TNF-α.35 Integrated analysis and clinical sample validation revealed significantly elevated JAK1 expression during IVDD progression. This finding indicates that activation of the JAK/STAT signaling pathway may contribute to the pathogenesis and progression of diverse diseases, including IVDD.36
WGCNA has the capability to identify gene sets that exhibit highly coordinated changes and to calculate the associations between these gene sets and phenotypes based on their internal connectivity, thereby delineating relationship patterns across numerous different samples.28 Consequently, by identifying significantly changing genes related to the EGFR-related KEGG pathways and the NPC degeneration process, and by incorporating susceptible modules uncovered by WGCNA, we preliminarily screened out five hub genes. To uncover patterns and associations within high-throughput sequencing data, we employed three machine learning techniques and ROC curve to identify key molecular features from complex datasets, identifying JAK1 as a potential key biomarker in NP degeneration, regulated by EGFR. Previous studies have reported that the JAK/STAT signaling pathway is activated in lumbar disc herniation and is involved in the abnormal degradation of collagen matrix, as well as the disruption of cell adhesion between NPCs.37–39 For instance, miR-22-3p can activate the expression of JAK1, promote the phosphorylation of STAT3, and facilitate its nuclear translocation, ultimately leading to the transcription of a series of degradative enzymes such as MMP13 and ADAMTS5.40,41 Additionally, the JAK1/STAT3 signaling pathway can regulate NLRP3-mediated pyroptosis, thereby influencing ECM degradation and oxidative stress in IVDD.33 Therefore, Hu et al42 found that Ruxolitinib, as an effective and selective oral inhibitor of JAK1 and JAK2 protein kinases, can improve cell adhesion in degenerative NP, reduce ECM degradation, and maintain the basic cytoskeletal morphology. Previous studies also demonstrated that EGFR directly phosphorylates JAK1 via its kinase activity, thereby initiating the JAK1/STAT3 signaling cascade.43 Furthermore, EGFR sustains JAK1 activation through a STAT3-mediated positive feedback loop.43,44 Consistently, our experiments using rat NPC revealed that EGFR signal inhibition reduces JAK1 expression. To validate this in vivo, we employed using Egfr CKOCol2 mice to construct an IVDD model, found that increased JAK1 expression in degenerated NP regions decreases with the inactivation of EGFR. Therefore, EGFR may influence the degeneration process of NPC by regulating JAK1.
Given that the EGFR-JAK1 axis may play a critical role in the progression of IVDD, this pathway could serve as a potential therapeutic target for drug development. However, it is important to note that certain members of the JAK family, SRC family, and MAPK pathway possess kinase domains that are highly homologous to that of JAK1.45,46 Therefore, to reduce off-target effects caused by the non-specific inhibition of JAK family members or other kinases, JAK1 inhibitors can be designed by leveraging JAK1-specific binding pockets and introducing steric hindrance at key residues.47 Furthermore, JAK1 is involved in critical immune signaling pathways mediated by interferons (IFN-α/β/γ) and IL-6, and excessive inhibition may increase susceptibility to infections and impair NK cell-mediated antitumor activity.46,48,49 Systemic risks related to lipid metabolism and cardiovascular function also warrant attention, as the JAK1/STAT3 pathway regulates lipid metabolism and coagulation, and its suppression may elevate low-density lipoprotein (LDL) levels and thrombosis risk.48 Therefore, in addition to addressing off-target concerns, the development of targeted delivery systems that direct therapeutics specifically to NPC and their EGFR signaling could enhance treatment precision and reduce systemic side effects.
There are several minor shortcomings in our study. Firstly, while we have used high-throughput sequencing data and various machine learning methods to enhance the reliability of our findings, validation in human-derived NPC is still required. Secondly, given the close association between IVDD and spatial architecture as well as biomechanics, future studies are warranted to further explore spatial transcriptomics and subpopulation analyses during the degeneration process.
Conclusion
In summary, we found that during the process of IVDD, some NPC undergo gradual degeneration, which is linked to cellular autophagy and apoptosis. From the perspective of cellular function and molecular characteristics, JAK1 expression increases during the degeneration of NPC and is regulated by EGFR. Therefore, our research provides a solid foundation for considering JAK1 as a key molecule in delaying NPC degeneration and offers new directions for further studies aimed at mitigating IVDD.
Data Sharing Statement
The datasets supporting the conclusions of this article are available in the GEO repository (https://www.ncbi.nlm.nih.gov/geo/), with accession numbers: GSE165722; GSE70362; GSE244889; GSE186542.
Ethics Approval and Consent to Participate
The collection of human NP tissues and use of publicly available human data were approved by the Ethics Committee of Union Hospital in Tongji Medical College ([2025]0423). All animal work performed in this report was approved by the Animal Experimentation Committee of Huazhong University of Science and Technology (IACUC number: 4085). Studies involving animal experimentation adhered to all relevant ethical regulations.
Funding
This work was supported by the National Key R&D Program of China (2021YFA1102600), the National Natural Science Foundation of China (82372407, 82002315, 82472435), the Fundamental Research Funds for the Central Universities (2024BRB021), and Wuhan Knowledge Innovation Project (2023020201020228).
Disclosure
The authors declare that they have no competing interests.
References
1. Kang L, Zhang H, Jia C, Zhang R, Shen C. Epigenetic modifications of inflammation in intervertebral disc degeneration. Ageing Res Rev. 2023;87:101902. doi:10.1016/j.arr.2023.101902
2. Dowdell J, Erwin M, Choma T, Vaccaro A, Iatridis J, Cho SK. Intervertebral disk degeneration and repair. Neurosurgery. 2017;80(3S):S46–S54. doi:10.1093/neuros/nyw078
3. Mwale F, Wang HT, Roughley P, Antoniou J, Haglund L. Link N and mesenchymal stem cells can induce regeneration of the early degenerate intervertebral disc. Tissue Eng Part A. 2014;20(21–22):2942–2949. doi:10.1089/ten.TEA.2013.0749
4. Gao B, Jiang B, Xing W, Xie Z, Luo Z, Zou W. Discovery and application of postnatal nucleus pulposus progenitors essential for intervertebral disc homeostasis and degeneration. Adv Sci. 2022;9(13):e2104888. doi:10.1002/advs.202104888
5. Ling Z, Liu Y, Wang Z, et al. Single-cell RNA-seq analysis reveals macrophage involved in the progression of human intervertebral disc degeneration. Front Cell Dev Biol. 2021;9:833420. doi:10.3389/fcell.2021.833420
6. Zhang Y, Han S, Kong M, Tu Q, Zhang L, Ma X. Single-cell RNA-seq analysis identifies unique chondrocyte subsets and reveals involvement of ferroptosis in human intervertebral disc degeneration. Osteoarthritis Cartilage. 2021;29(9):1324–1334. doi:10.1016/j.joca.2021.06.010
7. 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;9(1):37. doi:10.1038/s41413-021-00163-z
8. Wang D, Li Z, Huang W, et al. Single-cell transcriptomics reveals heterogeneity and intercellular crosstalk in human intervertebral disc degeneration. iScience. 2023;26(5):106692. doi:10.1016/j.isci.2023.106692
9. Xu C, Luo S, Wei L, et al. Integrated transcriptome and proteome analyses identify novel regulatory network of nucleus pulposus cells in intervertebral disc degeneration. BMC Med Genomics. 2021;14(1):40. doi:10.1186/s12920-021-00889-z
10. Sun X, Liu Y. Matrix metalloproteinase-10 in kidney injury repair and disease. Int J Mol Sci. 2022;23(4). doi:10.3390/ijms23042131
11. Shepard JB, Jeong JW, Maihle NJ, O’Brien S, Dealy CN. Transient anabolic effects accompany epidermal growth factor receptor signal activation in articular cartilage in vivo. Arthritis Res Ther. 2013;15(3):R60. doi:10.1186/ar4233
12. Zhang YW, Su Y, Lanning N, et al. Targeted disruption of Mig-6 in the mouse genome leads to early onset degenerative joint disease. Proc Natl Acad Sci USA. 2005;102(33):11740–11745. doi:10.1073/pnas.0505171102
13. Dealy CN, Scranton V, Cheng HC. Roles of transforming growth factor-alpha and epidermal growth factor in chick limb development. Dev Biol. 1998;202(1):43–55. doi:10.1006/dbio.1998.8988
14. Yoon YM, Oh CD, Kim DY, et al. Epidermal growth factor negatively regulates chondrogenesis of mesenchymal cells by modulating the protein kinase C-alpha, Erk-1, and p38 MAPK signaling pathways. J Biol Chem. 2000;275(16):12353–12359. doi:10.1074/jbc.275.16.12353
15. Appleton CT, Usmani SE, Bernier SM, Aigner T, Beier F. Transforming growth factor α suppression of articular chondrocyte phenotype and Sox9 expression in a rat model of osteoarthritis. Arthritis Rheum. 2007;56(11):3693–3705. doi:10.1002/art.22968
16. Usmani SE, Appleton CT, Beier F. Transforming growth factor-alpha induces endothelin receptor A expression in osteoarthritis. J Orthop Res. 2012;30(9):1391–1397. doi:10.1002/jor.22099
17. Zhang X, Siclari VA, Lan S, et al. The critical role of the epidermal growth factor receptor in endochondral ossification. J Bone Miner Res. 2011;26(11):2622–2633. doi:10.1002/jbmr.502
18. Pan Z, Sun H, Xie B, et al. Therapeutic effects of gefitinib-encapsulated thermosensitive injectable hydrogel in intervertebral disc degeneration. Biomaterials. 2018;160:56–68. doi:10.1016/j.biomaterials.2018.01.016
19. Chen S, Lv X, Hu B, et al. Critical contribution of RIPK1 mediated mitochondrial dysfunction and oxidative stress to compression-induced rat nucleus pulposus cells necroptosis and apoptosis. Apoptosis. 2018;23(5–6):299–313. doi:10.1007/s10495-018-1455-x
20. Ovchinnikov DA, Deng JM, Ogunrinu G, Behringer RR. Col2a1-directed expression of Cre recombinase in differentiating chondrocytes in transgenic mice. Genesis. 2000;26(2):145–146. doi:10.1002/(SICI)1526-968X(200002)26:2<145::AID-GENE14>3.0.CO;2-C
21. Xiang Z, Zhang P, Jia C, et al. Piezo1 channel exaggerates ferroptosis of nucleus pulposus cells by mediating mechanical stress-induced iron influx. Bone Res. 2024;12(1):20. doi:10.1038/s41413-024-00317-9
22. Wang D, Peng P, Dudek M, et al. Restoring the dampened expression of the core clock molecule BMAL1 protects against compression-induced intervertebral disc degeneration. Bone Res. 2022;10(1):20. doi:10.1038/s41413-022-00187-z
23. Rutges JP, Duit RA, Kummer JA, et al. A validated new histological classification for intervertebral disc degeneration. Osteoarthritis Cartilage. 2013;21(12):2039–2047. doi:10.1016/j.joca.2013.10.001
24. Chen F, Lei L, Chen S, et al. Serglycin secreted by late-stage nucleus pulposus cells is a biomarker of intervertebral disc degeneration. Nat Commun. 2024;15(1):47. doi:10.1038/s41467-023-44313-9
25. Tu J, Li W, Yang S, et al. Single-cell transcriptome profiling reveals multicellular ecosystem of nucleus pulposus during degeneration progression. Adv Sci. 2022;9(3):e2103631. doi:10.1002/advs.202103631
26. Rudolphi K, Gerwin N, Verzijl N, van der Kraan P, van den Berg W. Pralnacasan, an inhibitor of interleukin-1β converting enzyme, reduces joint damage in two murine models of osteoarthritis. Osteoarthritis Cartilage. 2003;11(10):738–746. doi:10.1016/s1063-4584(03)00153-5
27. Lakstins K, Yeater T, Arnold L, Khan S, Hoyland JA, Purmessur D. Investigating the role of culture conditions on hypertrophic differentiation in human cartilage endplate cells. J Orthop Res. 2021;39(6):1204–1216. doi:10.1002/jor.24692
28. Ding X, Qin J, Huang F, Feng F, Luo L. The combination of machine learning and untargeted metabolomics identifies the lipid metabolism -related gene CH25H as a potential biomarker in asthma. Inflamm Res. 2023;72(5):1099–1119. doi:10.1007/s00011-023-01732-0
29. Zhang W, Yang H, Zhu L, Luo Y, Nie L, Li G. Role of EGFR/ErbB2 and PI(3)K/AKT/e-NOS in lycium barbarum polysaccharides ameliorating endothelial dysfunction induced by oxidative stress. Am J Chin Med. 2019;47(7):1523–1539. doi:10.1142/S0192415X19500782
30. Du R, Xiao N, Han L, et al. Dexrazoxane inhibits the growth of esophageal squamous cell carcinoma by attenuating SDCBP/MDA-9/syntenin-mediated EGFR-PI3K-Akt pathway activation. Sci Rep. 2024;14(1):9167. doi:10.1038/s41598-024-59665-5
31. Wang Z, Ma J, Sun Y, et al. Isorhapontigenin delays senescence and matrix degradation of nucleus pulposus cells via PI3K/AKT/mTOR-mediated autophagy pathway in vitro and alleviates intervertebral disc degeneration in vivo. Int Immunopharmacol. 2024;139:112717. doi:10.1016/j.intimp.2024.112717
32. Liang ZH, Song J, Shangguan WJ, Zhang QQ, Shao J, Zhang YH. Melatonin mitigates matrix stiffness-induced intervertebral disk degeneration by inhibiting reactive oxygen species and melatonin receptors mediated PI3K/AKT/NF-kappaB pathway. Am J Physiol Cell Physiol. 2024;327(5):C1236–C1248. doi:10.1152/ajpcell.00630.2023
33. Yu P, Ma Z, Jiang H, Liu J, Li H. Sinapine thiocyanate alleviates intervertebral disc degeneration by not regulating JAK1/STAT3/NLRP3 signal pathway. Adv Clin Exp Med. 2024;33(9):965–977. doi:10.17219/acem/174508
34. Huang W-L, Yeh -H-H, Lin -C-C, et al. Signal transducer and activator of transcription 3 activation up-regulates interleukin-6 autocrine production: a biochemical and genetic study of established cancer cell lines and clinical isolated human cancer cells. Mol Cancer. 2010;9(1):309. doi:10.1186/1476-4598-9-309
35. Wu C, Ge J, Yang M, et al. Resveratrol protects human nucleus pulposus cells from degeneration by blocking IL-6/JAK/STAT3 pathway. Euro J Med Res. 2021;26(1):81. doi:10.1186/s40001-021-00555-1
36. Xin P, Xu X, Deng C, et al. The role of JAK/STAT signaling pathway and its inhibitors in diseases. Int Immunopharmacol. 2020;80:106210. doi:10.1016/j.intimp.2020.106210
37. Kazezian Z, Li Z, Alini M, Grad S, Pandit A. Injectable hyaluronic acid down-regulates interferon signaling molecules, IGFBP3 and IFIT3 in the bovine intervertebral disc. Acta Biomater. 2017;52:118–129. doi:10.1016/j.actbio.2016.12.029
38. Osuka K, Usuda N, Aoyama M, et al. Expression of the JAK/STAT3/SOCS3 signaling pathway in herniated lumbar discs. Neurosci Lett. 2014;569:55–58. doi:10.1016/j.neulet.2014.03.045
39. Shouda T, Yoshida T, Hanada T, et al. Induction of the cytokine signal regulator SOCS3/CIS3 as a therapeutic strategy for treating inflammatory arthritis. J Clin Invest. 2001;108(12):1781–1788. doi:10.1172/JCI13568
40. Heinrich PC, Behrmann I, Haan S, Hermanns HM, Muller-Newen G, Schaper F. Principles of interleukin (IL)-6-type cytokine signalling and its regulation. Biochem J. 2003;374(Pt 1):1–20. doi:10.1042/BJ20030407
41. Chen Z, Liao Z, Liu M, et al. Nucleus pulposus-targeting nanocarriers facilitate mirna-based therapeutics for intervertebral disc degeneration. Adv Healthc Mater. 2023;12(28):e2301337. doi:10.1002/adhm.202301337
42. Hu F, Pan Z, Liu C, et al. Identification of inflammatory regulation roles of thalidomide/ruxolitinib in nucleus pulposus and construction of polyelectrolyte nanocomplexes-impregnated injectable hydrogels for synergistic intervertebral disk degeneration treatment. Nano Today. 2022;44:101462. doi:10.1016/j.nantod.2022.101462
43. Kang JH, Uddin N, Kim S, et al. Tumor-intrinsic role of ICAM-1 in driving metastatic progression of triple-negative breast cancer through direct interaction with EGFR. Mol Cancer. 2024;23(1):230. doi:10.1186/s12943-024-02150-4
44. Palmer MA, Kirchhoff R, Buerger C, Benatzy Y, Schebb NH, Brune B. RNAi-based ALOX15B silencing augments keratinocyte inflammation in vitro via EGFR/STAT1/JAK1 signalling. Cell Death Dis. 2025;16(1):39. doi:10.1038/s41419-025-07357-x
45. Song Y, Yoon DH, Yang H, et al. Phase I dose escalation and expansion study of golidocitinib, a highly selective JAK1 inhibitor, in relapsed or refractory peripheral T-cell lymphomas. Ann Oncol. 2023;34(11):1055–1063. doi:10.1016/j.annonc.2023.08.013
46. Song Y, Malpica L, Cai Q, et al. Golidocitinib, a selective JAK1 tyrosine-kinase inhibitor, in patients with refractory or relapsed peripheral T-cell lymphoma (JACKPOT8 Part B): a single-arm, multinational, Phase 2 study. Lancet Oncol. 2024;25(1):117–125. doi:10.1016/S1470-2045(23)00589-2
47. Zhou S, Mao W, Su Y, et al. Identification of TUL01101: a novel potent and selective JAK1 inhibitor for the treatment of rheumatoid arthritis. J Med Chem. 2022;65(24):16716–16740. doi:10.1021/acs.jmedchem.2c01550
48. Winthrop KL. The emerging safety profile of JAK inhibitors in rheumatic disease. Nat Rev Rheumatol. 2017;13(4):234–243. doi:10.1038/nrrheum.2017.23
49. Xu Q, He L, Yin Y. Risk of herpes zoster associated with JAK inhibitors in immune-mediated inflammatory diseases: a systematic review and network meta-analysis. Front Pharmacol. 2023;14:1241954. doi:10.3389/fphar.2023.1241954