Upregulation of cGMP-dependent Protein Kinase (PRKG1) in the Development of Adolescent Idiopathic Scoliosis.

OBJECTIVE: To explore the molecular regulatory mechanisms underlying fibroblast differentiation and dysfunction in the development of adolescent idiopathic scoliosis (AIS) in an effort to identify candidate therapeutic targets for AIS.
METHODS: The GSE110359 dataset, obtained from the bone marrow stromal cells of 12 AIS patients and five healthy controls, was retrieved from the GEO database. The data were preprocessed and differentially expressed genes (DEGs) were identified. KEGG pathway and Gene Ontology (GO)-Biological Process (BP) enrichment analyses were performed to identify the function of the DEGs. A protein-protein interaction (PPI) and a microRNA-transcription factor (TF)-target co-regulatory network were constructed to identify hub genes in the development of AIS. In addition, hub DEGs were evaluated by quantitative PCR (qPCR) and immunohistochemical staining.
RESULTS: A total of 188 DEGs including 100 up-regulated and 88 down-regulated genes were obtained. The up-regulated DEGs were related to "p53 signaling pathway", "FoxO signaling pathway", and "cGMP-PKG signaling pathway" terms, while the down-regulated DEGs were significantly enriched in seven terms including "protein processing in endoplasmic reticulum". The key up-regulated genes, PRKG1, CCNG2, and KAT2B, and the key down-regulated genes, MAP2K1 and DUSP6, were identified by the PPI and miRNA-TF-Target regulatory network analyses. mRNA expression patterns for PRKG1, DUSP6, and KAT2B were successfully verified by qPCR. In addition, PRKG1 protein levels were found to be elevated during the immunohistochemical analysis.
CONCLUSION: Increased expression of PRKG1 in AIS patients might be an attractive therapeutic target for AIS. However, further gain or loss-of-function studies should be conducted.

adolescent idiopathic scoliosis

differentially expressed genes

Gene Ontology

Biological Process

protein–protein interaction

transcription factor

quantitative PCR

phosphoinositide 3‐kinase

marrow stromal cells

Gene Expression Omnibus

Overrepresentation Enrichment Analysis

Quantitative Reverse Transcription PCR

standard deviation

protein kinase cGMP‐dependent 1

Introduction

Adolescent idiopathic scoliosis (AIS) is the most common structural spinal deformity when the deformity exceeds 10 degrees . It affects approximately 2% to 4% of adolescents between 10 and 18 years of age , . Scoliosis surgery has a high rate of complication making it high‐risk and expensive, which can impose a heavy economic burden on families and society. Therefore, understanding the molecular mechanism underlying AIS development is of vital importance and could help develop novel medical interventions.

Currently, the molecular mechanism of AIS has not been determined. It is believed that the occurrence and development of AIS involve multiple cell types and genes. Studies suggest that joint laxity is associated with AIS , and the incidence of AIS in ballet dancers, who have a higher rate of hypermobility(70%), is significantly higher than that of non‐dancers of the same age . The distribution of collagen in adolescents with AIS is different from that of normal subjects . In addition, connective tissue abnormalities including Marfan syndrome and Ehlers‐Danlos syndrome are often associated with a higher incidence of spinal deformity , . These studies suggest that the pathogenesis of AIS may be closely related to collagen abnormalities. Collagen formation is mainly achieved by fibroblasts producing collagen‐I, α‐SMA, and fibronectin. Haller et al. found that variants in musculoskeletal collagen genes were significantly enriched in AIS compared to controls. They revealed that the fibrillar collagen gene COL11A2 was highly associated with AIS . In addition, many studies have suggested that TGF‐β/Smad4 signaling, MMP‐1, 2, 9, 13, Ras‐MAPK signaling, and phosphoinositide 3‐kinase (PI3K)‐Akt‐eNOS‐NO signaling are the main regulators of collagen synthesis and degradation , , , .

Fibroblasts were differentiated from bone marrow stromal cells (MSC). The role of fibroblast differentiation and dysfunction in the pathogenesis of scoliosis has received extensive attention. For example, several genes including TGF‐β1 and MMP3 were differentially expressed in the dense connective tissues of patients with AIS and could participate in the development of scoliosis , . However, current research on the role of fibroblast differentiation and dysfunction in the development of scoliosis remains unclear; this prevents them from being used as targets for clinical intervention.

Microarray‐based studies have been widely conducted to investigate the molecular mechanism of AIS development. For example, Fendri et al. identified 145 differentially expressed genes (DEGs) between AIS patients and healthy control, including FAM101A, ZIC2, PITX1, and COMP . This study aimed to: (i) explore the molecular regulatory mechanisms underlying fibroblast differentiation and dysfunction in the development of AIS; and (ii) identify candidate therapeutic targets for AIS. DEGs involved in collagen‐related biological processes could be used as candidate therapeutic targets for the treatment of AIS in the future.

Methods

Microarray Data

In this study, we systematically screened DEGs in fibroblasts from normal and AIS patients, and explore the molecular mechanism of AIS pathogenesis using bioinformatic analyses (Fig. 1). The GSE110359 dataset was downloaded from the NCBI GEO database (Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo/) , (species: Homo sapiens). The data were obtained from MSCs from 12 AIS patients and five healthy controls. All subjects were less than 18 years of age. This data was generated on a GPL17586[HTA‐2_0] platform using the Affymetrix Human Transcriptome Array 2.0 [transcript (gene) version].

Fig 1

Schematic representation of the study design

The original CEL data were preprocessed and normalized using the Oligo program in R (version 1.36.1, http://bioconductor.org/packages/release/bioc/html/oligo.html). The probe was annotated in the platform annotation file and was used to remove probes that did not match the Gene symbol; for different probes that mapped to the same gene, the mean value for these probes was taken as the final expression value of the gene.

Differential expression analysis was conducted using the classical Bayesian method from the limma package (Version 3.10.3, http://www.bioconductor.org/packages/2.9/bioc/html/limma.html). Genes with a P value <0.01 were identified as DEGs.

Enrichment Analysis of

DAVID (Version6.8, https://david-d.ncifcrf.gov/) was used to perform KEGG pathway, Gene Ontology (GO)‐Biological Process (BP) enrichment analysis for upregulated and downregulated genes. The functional terms with enriched gene counts ≥2 and P values < 0.05 were considered significant.

Protein–Protein Interaction (

The interactions between the proteins encoded by the DEGs were analyzed using the STRING (Version: 10.0, http://www.string-db.org/) database . The DEGs acted as the input gene set and the species was set as human. To obtain the highest number of interactions, the PPI score was set to 0.15 (low confidence). The PPI network was built using Cytoscape (version: 3.2.0, http://www.cytoscape.org/).

Proteins in the same module of the PPI network often have the same or similar functions. We used Cytoscape's MCODE (Version 1.4.2, http://apps.cytoscape.org/apps/MCODE) plugin to identify the most significant clustering modules in the PPI network. The threshold value was a score > 5. KEGG pathway enrichment analysis of the genes in significant clustering modules was then performed.

Transcription Factor (

The miRNA‐target and TF‐target enrichment prediction was performed using the Overrepresentation Enrichment Analysis (ORA) enrichment method from the WEB‐based Gene SeT AnaLysis Toolkit (Web Gestalt, http://www.webgestalt.org/) . Only genes from significant modules were used for this analysis. The thresholds were set at a count ≥2 and P value < 0.05. The TF‐target gene regulatory network and the miRNA‐target gene regulatory network were combined to obtain them iRNA‐TF‐Target regulatory network.

Quantitative Reverse Transcription

Ten female AIS patients (mean age = 16.3) and 10 female subjects with vertebral fractures or lumbar disc herniation (mean age = 16.5) were enrolled in this study, and gave informed consent. Total RNA was extracted by TRIzol Reagent (TaKaRa, Dalian, China) from the ligament tissues of AIS patients and healthy subjects (n = 3 for each) as per the commercial protocol. Total RNA was reverse transcribed to cDNA using prime Script RT Master MIX (TaKaRa) and expression levels forMAP2K1, PRKG1, KAT2B, DUSP6, and CCNG2 were detected with a Power SYBR Green PCR Master (Thermo Fisher) on an ABI 7500 (Applied Biosystems, Carlsbad, CA, USA). The relevant primer sequences can be found in Table 1. GAPDH was used as the housekeeping gene.

TABLE 1

The forward and reverse primers of genes and GAPDH normalizer for RT‐PCR

Primer	Direction	Sequence (5′‐3′)
MAP2K1	Forward	CAATGGCGGTGTGGTGTTC
	Reverse	GATTGCGGGTTTGATCTCCAG
PRKG1	Forward	GACAACGATGAACCACCAC
	Reverse	GCTTTGCTTCAGGACCAC
DUSP6	Forward	GAAATGGCGATCAGCAAGACG
	Reverse	CGACGACTCGTATAGCTCCTG
KAT2B	Forward	CGAATCGCCGTGAAGAAAGC
	Reverse	CTTGCAGGCGGAGTACACT
CCNG2	Forward	TAGCAGGGACAGTGGATA
	Reverse	AATACTTGGGCAATAGGA
GAPDH	Forward	TGACAACTTTGGTATCGTGGAAGG
	Reverse	AGGCAGGGATGATGTTCTGGAGAG

Immunohistochemical Localization of

Immunohistochemical analysis was performed as described previously . Briefly, ligament tissues were embedded and cut into 8 μm thick frozen sections. Antigen recovery was performed by incubating the sections in 0.01 M citric acid buffer and the activity of endogenous peroxidases were blocked with 3% hydrogen peroxide for 15 min, followed by incubation with 10% goat serum for 1 h at 37°C to block non‐specific binding reactions. The sections were incubated with a primary antibody against PRKG1 (1:100, Cat.no.: ab38007, Abcam) at 4°C overnight. We then added the goat anti‐Rabbit IgG‐HRP (H + L, 1: 400, Jackson ImmunoResearch, Cat.no.: 111–035‐003). Sections were developed in diaminobenzidine solution for 1 min and counterstained with hematoxylin for 30 s, dehydrated, and sealed with neutral gum.

Statistical Analysis

Data were expressed as mean ± standard deviation (SD). Experiments were conducted in triplicate. Gene expression differences detected by qRT‐PCR between the two groups were compared using the Student's t‐test in GraphPad prism 5.0 (GraphPad Software Inc., San Diego, CA, USA). P values < 0.05 were regarded as statistically significant.

Results

Differential Expression Analysis

The distribution of the normalized expression values are represented in Fig. 2A and show the median values on the same horizontal line. A total of 188 DEGs with P value <0.01 were identified. Of these, 100 were up‐regulated and 88 were down‐regulated (Fig. 2B and Table S1). The gene expression heatmap (Fig. 2C) and PCA plots (Fig. 2D) showed that DEGs can be clearly distinguished in the disease samples when compared to the control.

Fig 2

Box diagram after data preprocessing (A), volcano plot (B) heatmap diagram (C), and principal component analysis (PCA) plots (D) of gene expression profiles. Green arethe AIS samples, and yellow are the control samples

The KEGG pathway enrichment analysis showed that up‐regulated DEGs were significantly enriched in three pathways, including “p53 signaling pathway” (n = 3; P = 3.04 × 10−3), “FoxO signaling pathway” (n = 3; P = 1.84 × 10−2), and “cGMP‐PKG signaling pathway” (n = 3; P = 3.64 × 10−2). The down‐regulated DEGs were significantly enriched into seven pathways, including “Protein processing in endoplasmic reticulum” (n = 4; P = 2.78× 10−3) and “Complement and coagulation cascades” (n = 2; P = 3.25 × 10−2) (Table 2).

TABLE 2

KEGG pathway and Gene Ontology (GO)‐Biological Process (BP) analysis results of differentially expressed genes (DEGs)

	ID	Name	Count	P Value	Gene
Up‐regulated genes
Pathway	hsa04115	p53 signaling pathway	3	3.04E‐03	DDB2;STEAP3;CCNG2
	hsa04068	FoxO signaling pathway	3	1.84E‐02	FBXO32;RBL2;CCNG2
	hsa04022	cGMP‐PKG signaling pathway	3	3.64E‐02	GTF2I;PRKG1;PDE5A
GO‐BP	GO:0006357	regulation of transcription from RNA polymerase II promoter	8	4.06E‐03	CAMTA2, TSHZ1, ANKRA2, RBL2, GTF2I, TFDP2, VEZF1, BRD8
	GO:0006633	fatty acid biosynthetic process	3	2.37E‐02	ELOVL4, FAXDC2, CBR4
	GO:0006091	generation of precursor metabolites and energy	3	2.46E‐02	FECH, SLC25A27, DHTKD1
	GO:0007264	small GTPase mediated signal transduction	5	2.65E‐02	PLD2, RASL11A, ARL17A, RALBP1, DOCK11
	GO:0007165	signal transduction	11	3.88E‐02	RAP2B, RALBP1, GTF2I, PPP2R5D, PKIG, PDE5A, MCC, PRKG1, BRD8, DAPK1, DTNA
Down‐regulated genes
Pathway	hsa04141	Protein processing in endoplasmic reticulum	4	2.78E‐03	ERO1A;STT3A;DERL1;TXNDC5
	hsa05020	Prion diseases	2	7.98E‐03	LAMC1;MAP2K1
	hsa05110	Vibrio cholerae infection	2	1.58E‐02	KDELR2;ERO1A
	hsa00510	N‐Glycan biosynthesis	2	1.64E‐02	STT3A;MGAT5
	hsa00590	Arachidonic acid metabolism	2	2.20E‐02	PTGS1;CBR1
	hsa04610	Complement and coagulation cascades	2	3.24E‐02	A2M;MBL2
	hsa05222	Small cell lung cancer	2	4.25E‐02	LAMC1;TRAF3
GO‐BP	GO:0006457	protein folding	4	2.65E‐02	GRPEL1, ERO1A, TXNDC5, PDIA6
	GO:0034976	response to endoplasmic reticulum stress	3	2.95E‐02	ERO1A, TXNDC5, PDIA6
	GO:0022617	extracellular matrix disassembly	3	3.02E‐02	A2M, LAMC1, CTSG
	GO:0045454	cell redox homeostasis	3	3.10E‐02	ERO1A, TXNDC5, PDIA6
	GO:0019371	cyclooxygenase pathway	2	3.52E‐02	CBR1, PTGS1
	GO:0034975	protein folding in endoplasmic reticulum	2	4.55E‐02	ERO1A, EMC1

The GO‐BP enrichment analysis revealed that the up‐regulated DEGs were significantly enriched in five GO‐BP terms, and the down‐regulated DEGs were significantly enriched in six GO‐BP terms (Table 2).

The

To investigate the associations between the DEGs, we built a PPI network. As shown in Fig. 3A, the PPI network included 70 upregulated genes, 58 downregulated genes, and 483 pairs of interactions. Topological property of each node is displayed in Table S2. A total of 14 genes had a degree that was equal or larger than 20, including RAP2B (Upregulated gene, degree = 32), KAT2B (Upregulated gene, degree =27), PRKG1 (Upregulated gene, degree = 26), PDIA6 (Downregulated gene, degree = 25), ARL17A (Upregulated gene, degree =25), and ANKRA2 (Upregulated gene, degree =25). Since this network was too large and complex for more detailed analysis, we created a further mined sub‐network of modules. As a result, we created a sub‐network module (score = 6.471) which contained 18 nodes and 55 interaction pairs (Fig. 3B). The important upregulated (PRKG1, CCNG2, KAT2B, and RBL2) and down‐regulated (MAP2K1 and DUSP6) genes were included in this module. KEGG pathway analysis was performed on the genes in the significant clustering modules, and a total of nine KEGG pathway terms were shown to be enriched, including “FoxO signaling pathway” (MAP2K1, RBL2, and CCNG2), “vascular smooth muscle contraction” (PRKG1 and MAP2K1), “thyroid hormone signaling pathway” (MAP2K1 and KAT2B), “cGMP‐PKG signaling pathway” (PRKG1 and MAP2K1), “MAPK signaling pathway” (DUSP6 and MAP2K1), and “PI3K‐Akt signaling pathway” (MAP2K1 and RBL2) (Table 3).

Fig 3

Theprotein‐protein interaction (PPI)networkof differentially expressed genes in AIS (A).Thesignificant clustering modules of the PPI network (B). The circle nodes indicate the up‐regulated genes, and the diamond nodes represent the down‐regulated genes.

TABLE 3

KEGG pathway enriched by differentially expressed genes (DEGs) in the sub‐network module

ID	Name	Count	P Value	Genes
hsa04068	FoxO signaling pathway	3	2.17E‐04	MAP2K1;RBL2;CCNG2
hsa04730	Long‐term depression	2	1.53E‐03	PRKG1;MAP2K1
hsa04540	Gap junction	2	3.45E‐03	PRKG1;MAP2K1
hsa04270	Vascular smooth muscle contraction	2	5.98E‐03	PRKG1;MAP2K1
hsa04919	Thyroid hormone signaling pathway	2	6.09E‐03	MAP2K1;KAT2B
hsa04022	cGMP‐PKG signaling pathway	2	1.16E‐02	PRKG1;MAP2K1
hsa05203	Viral carcinogenesis	2	1.50E‐02	RBL2;KAT2B
hsa04010	MAPK signaling pathway	2	2.53E‐02	DUSP6;MAP2K1
hsa04151	PI3K‐Akt signaling pathway	2	4.45E‐02	MAP2K1;RBL2

A total of 21 miRNAs and six TFs that target the genes in the significant modules were predicted using WebGestalt, and a miRNA‐TF‐Target regulatory network involving these 21 miRNAs, six TFs, three down‐regulated DEGs, and 10 up‐regulated DEGs was then constructed (Fig. 4). For the DEGs related to vascular smooth muscle contraction, PRKG1 was targeted by miR‐199a, miR‐522 and OCT1/MAP2K1 were targeted by miR‐324‐3p, miR‐181a, miR‐181b, miR‐181c, miR‐181d, miR‐330, miR‐199a, and PAX4. Other upregulated genes including TNRC6, HECA, DOCK1, RASL11A, KLHL24, CCNG2, KAT2B, and RBL2, and the down‐regulated genes KLF10 and DUSP6 were also included in the network.

Fig 4

The microRNA‐transcription factor‐Target regulatory network. The triangle nodesaremiRNAs, the hexagon nodesare transcription factors, and edges with an arrow indicate regulatory relationships.

Five critical genes, MAP2K1, PRKG1, DUSP6, KAT2B, and CCNG2 were selected for qRT‐PCR validation. As shown in Fig. 5, qRT‐PCR showed that the mRNA expression levels of PRKG1, DUSP6, and KAT2B were the same as those predicted by the DEG analysis from the microarray data (P < 0.05). Since PRKG1 is involved in “Vascular smooth muscle contraction”, we further characterized its expression by immunohistochemistry. The results showed that the expression levels of PRKG1 protein in the AIS group were also elevated when compared with the control group (Fig. 6).

Fig 5

The qRT‐PCR analysis of five crucial genes, including MAP2K1, PRKG1, KAT2B, CCNG2 and DUP6.Data were expressed as mean ± standard deviation. Comparisons between groups were calculated by student's t test. *, P < 0.05; **, P < 0.01***, P < 0.001.

Fig 6

The protein expression of PRKG1 detected by immunohistochemicalanalysis (×200). The arrow indicatesPRKG1 expression.

Discussion

AIS is the most common structural spinal deformity. However, its pathological mechanism is still largely unknown. In this study, we retrieved an AIS microarray dataset and used it to investigate the potential molecular mediators of its pathological mechanism. As a result, we identified a total of 188 DEGs between AIS and control samples. Further analysis suggested that the “p53 signaling pathway”, “FoxO signaling pathway”, and “cGMP‐PKG signaling pathway” were activated while “Protein processing in endoplasmic reticulum” and “Complement and coagulation cascades” pathways were suppressed in AIS patients. Five crucial nodes in the PPI and miRNA‐TF‐target networks were identified. These included MAP2K1, PRKG1, DUSP6, KAT2B, and CCNG2 which were further verified using qRT‐PCR and we were able to confirm the expression profiles of PRKG1, DUSP6, and KAT2B. In addition, immunohistochemical analysis was carried out to evaluate the protein expression of PRKG1 in AIS patients.

PRKG1 (protein kinase cGMP‐dependent 1) was significantly enriched in “cGMP‐PKG signaling”, “Gap junction”, and “Vascular smooth muscle contraction” and may be important in the complex pathological mechanisms of AIS. PRKG1 acts as key mediator of the nitric oxide/cGMP signaling pathway and is an important component in many signal transduction processes in diverse cell types. PRKG1 has been shown to be associated with RhoA activity and Akt signaling and was identified as a novel negative regulator of osteoblast differentiation in human skeletal (mesenchymal) stem cell . In our study, PRKG1 was shown to be upregulated at both the transcript and protein levels, and was linked to the following KEGG pathways “cGMP‐PKG signaling”, “Long‐term depression”, “Gap junction”, and “Vascular smooth muscle contraction”, as well as the GO‐BP term “signal transduction”. Hewitson et al. indicated that treatment with dipyridamole decreased fibroblast mitogenesis and collagen synthesis via the cGMP‐PKG pathway. Phosphodiesterase inhibitors, which increase cellular cGMP concentrations, induce the synthesis of collagen in keloid fibroblasts . Collagen‐induced platelet secretion and aggregation is also mediated through cAMP and cGMP‐independent signaling .

In addition, MAP2K1, which was shown to be down regulated in the expression study, interacts directly with PRKG1, and is involved in cGMP‐PKG signaling, as well as MAPK and PI3K‐Akt signaling. Mutant MAPK7‐induced AIS has been linked to impaired osteogenesis . Previous study has indicated that p38 and MAPK play important roles in the production of collagen in fibroblasts . Moreover, the activation of MAPK signaling and the PI3K/Akt pathway promote fibroblast migration . The activation ofPI3K/Akt pathway is associated with the normal regulation of MSC functions and has applications in tissue regeneration , . Degeneration of the intervertebral discs, which are primarily composed of water, aggrecan, and fibrillar collagens, could be prevented by the activation of the PI3K/Akt pathway as it increases extracellular matrix production, preventing apoptosis and oxidative damage, and facilitating cell proliferation . Therefore, MAP2K1 interactions with PRKG1 may indicate that PRKG1 is also an important candidate in the pathogenesis of AIS.

However, this study inevitably has some limitations. First, the sample size of this study is relatively small. Only 12 AIS patients and five healthy controls were analyzed from microarray study and key genes were validated in 10 AIS patients and 10 controls. Second, the molecular mechanism of PRKG1 in the development of AIS still needs further investigation. Therefore, further in vitro and in vivo studies should be conducted.

Conclusion

In conclusion, we identified a total of 188 DEGs following a comparison of expression data from AIS and healthy control samples. Further analysis suggested that the “p53 signaling pathway”, “FoxO signaling pathway”, and “cGMP‐PKG signaling pathway” were activated, while “Protein processing in endoplasmic reticulum” and “Complement and coagulation cascades” pathways were suppressed in AIS patients. Increased expression of PRKG1 in AIS patients was confirmed by qRT‐PCR and immunohistochemical staining, indicating that it may be an attractive therapeutic target for novel intervention strategies in the treatment of AIS. However, further gain and loss of function studies should be conducted to evaluate this hypothesis.

Table S1 The list of 188 differentially expressed genes.

Click here for additional data file.

Table S2 Topological property of each node of the protein–protein interaction network.

Click here for additional data file.