Conjoint research of WGCNA, single-cell transcriptome and structural biology reveals the potential targets of IDD development and treatment and JAK3 involvement

Abstract

Objectives: This study conducted integrated analysis of bulk RNA sequencing, single-cell RNA sequencing and Weighted Gene Co-expression Network Analysis (WGCNA), to comprehensively decode the most essential genes of intervertebral disc degeneration (IDD); then mainly focused on the JAK3 macromolecule to identify natural compounds to provide more candidate drug options in alleviating IDD.

Methods: In the first part, we performed single-cell transcriptome analysis and WGCNA workflow to delineate the most pivotal genes of IDD. Then series of structural biology approaches and high-throughput virtual screening techniques were performed to discover potential compounds targeting JAK-STAT signaling pathway, such as Libdock, ADMET, precise molecular docking algorithm and in-vivo drug stability assessment.

Results: Totally 4 hub genes were determined in the development of IDD, namely VEGFA, MMP3, TNFSF11, and TIMP3, respectively. Then, 3 novel natural materials, ZINC000014952116, ZINC000003938642 and ZINC000072131515, were determined as potential compounds, with less toxicities and moderate ADME characteristics. In-vivo drug stability assessment suggested that these drugs could interact with JAK3, and their ligand-JAK3 complexes maintained the homeostasis in-vivo, which acted as regulatory role to JAK3 protein. Among them, ZINC000072131515, also known as Menaquinone, demonstrated significant protective roles to alleviate the progression of IDD in vitro, which proved the nutritional therapy in alleviating IDD.

Conclusions: This study reported the essential genes in the development of IDD, and also the roles of Menaquinone to ameliorate IDD through inhibiting JAK3 protein. This study also provided more options and resources on JAK3 targeted screening, which may further expand the drug resources in the pharmaceutical market.

Keywords: Janus kinase family; high-throughput virtual screening; intervertebral disc degeneration; single-cell transcriptome; weighted gene co-expression network analysis.

Figure 1

(A) Sample clustering analysis to detect the outliers. (B) Determination of soft threshold power value. Left panel indicated scale-free model fit index; right panel indicated the mean connectivity of these values. (C) Dendrogram branch plot of genes based on dissimilarity measure and assignment modules. (D) Module-trait correlation heatmap between different clinical traits and modules. (E) Gene significance histogram plot of all clustered modules. (F) Eigenvalue correlation heatmap of the modules and clinical traits. (G) Correlation scatter plot between gene-significance and module membership.

Figure 2

(A) Volcano scatter plot of the genes analyzed by limma algorithm. (B) Biological functions about the up-regulated DEGs. (C, D) Aberrantly activated signaling pathways analyzed by KEGG, GSEA, respectively.

Figure 3

(A) tSNE visualization of all cells in NP tissues after sample integration, 9 clustered were generated. (B) Dot plot showing the NP cells marker genes ACAN, SOX9 and MIA within tSNE map; and violin plot showing the expression of macrophage marker genes CD14 and MRC1. (C) tSNE visualization of all cells in NP tissues between degeneration and non-degeneration group. (D) Selection of the most variable genes for monocle analysis. (E) PCA dimensional reduction by monocle algorithm between degeneration and non-degeneration group. (F) Volcano scatter plot of the genes analyzed by monocle method. (G) Hierarchical clustering heatmap of the analyzed DEGs in each cell.

Figure 4

(A) Venn diagram suggesting 33 genes were commonly expressed genes in three parts. (B) PPI construction of the commonly expressed genes to identify the most significant hub genes. (C) Featureplot visualization about the hub genes between degeneration and non-degeneration group at single-cell resolution. (D) Relative mRNA expression levels of the four hub genes in the control and degeneration groups.

Figure 5

The chemical structure of Janus kinase 3 (JAK3). (A) Initial crystal structure of JAK3 with active binding sphere addition. (B) Macromolecule with surface binding region surrounded, red indicated positive charge and blue indicated negative charge. (C) The Ramachandran chart of JAK3 macromolecule. (D) Docking structure model of FM-381-JAK3 complex.

Figure 6

(A) Residues interaction roles heatmap based on the top 50 compounds from Libdock module. (B–D) Detailed intermolecular interactions between ZINC000014952116-JAK3, ZINC000003938642-JAK3 and ZINC000072131515-JAK3, respectively. (E–G) Residues interaction roles heatmap of ZINC000014952116-JAK3, ZINC000003938642-JAK3, and ZINC000072131515-JAK3 complex.

Figure 7

Schematic drawing of interactions between ligands and JAK3, the surface of binding area was added, blue represented positive charge, red represented negative charge, and ligands were shown in sticks, the structure around the ligand-receptor junction was shown in thinner sticks. (A) ZINC000014952116-JAK3 complex, (B) ZINC000003938642-JAK3, and (C) ZINC000072131515-JAK3 complex.

Figure 8

Probability of interactive residues roles between the compounds and JAK3 based on (A) favorable count, (B) unfavorable count, (C) hydrogen bond count, (D) charge count, (E) hydrophobic count and (F) other count.

Figure 9

Results of molecular dynamics simulation of these two compounds. (A) Orthorhombic box with an explicit periodic boundary solvation water model. (B–D) Energy values of ZINC000014952116-JAK3, ZINC000003938642-JAK3 complex, and ZINC000072131515-JAK3 complex. (E) RMSD values of these three compounds.

Figure 10

(A) NP cells were treated with different concentrations (50, 100 μM) of TBHP after 24 h and imaged by phase-contrast microscopy. (B) Cellular viability of the NP cells with different concentration of TBHP for 24 h. (C) Cellular viability of the TBHP-induced (100 μM) NP cells with different concentration of ZINC000072131515. (D, E) The protein expression of JAK3 in NP cells treated with different concentrations of ZINC000072131515. (F) NP cells migration was determined by scratch experiments, the results were recorded at 0, 24 h. (G) Percentage of wound closure of NP cells at 24 h. *P < 0.05; **P < 0.01; ***P < 0.0001; ns, non-significance.