RNA-seq based integrative analysis of potential crucial genes and pathways associated with patellar instability.

Patellar instability (PI) is a common knee injury in adolescents, but the crucial biomarkers and molecular mechanisms associated with it remain unclear. We established a PI mouse model and investigated PI-related changes in gene expression by RNA sequencing (RNA-seq). Differentially expressed gene (DEG) analysis and enrichment analysis were performed to identify crucial genes and pathways associated with PI. Subsequently, a protein-protein interaction, DEG-miRNA, DEG-transcription factors, and DEG-drug interaction networks were constructed to reveal hub genes, molecular mechanism, and potential drugs for PI. Finally, the reliability of the sequencing results was confirmed by real-time quantitative polymerase chain reaction (RT-qPCR) and immunohistochemistry. Upon comparison with the control group, 69 genes were differently expressed in PI, including 17 upregulated and 52 downregulated ones. The DEGs were significantly enriched in Janus kinase (JAK)/signal transducer and activator of transcription (STAT) signaling pathway and immune responses. The protein-protein interaction network identified ten PI-related hub genes, all of which are involved in the JAK/STAT signaling pathway or inflammation-related pathways. DEG-miRNA and DEG-transcription factor networks offered new insights for regulating DEGs post-transcriptionally. We also determined potential therapeutic drugs or molecular compounds that could restore dysregulated expression of DEGs via the DGIdb database. RT-qPCR results were consistent with the RNA-seq, confirming the reliability of the sequencing data. Immunohistochemistry results suggested that JAK1 and STAT3 expression was increased in PI. Our study explored the potential molecular mechanisms in PI, provided promising biomarkers and suggested a molecular basis for therapeutic targets for this condition.

Highlights

Our study explored the potential molecular mechanisms in PI for the first time.

We provided the promising biomarkers of PI.

We suggested a molecular basis for therapeutic targets.

Introduction

Patellar instability (PI) is a common knee injury in adolescents [1]. Its incidence was estimated to be 29–104/100,000 individuals per year, with females experiencing the highest incidence [2,3]. Owing to abnormal loading in the femoral trochlea, primary PI can cause bone loss and trochlear dysplasia [4,5]. It has also been found that PI can aggravate cartilage degeneration, and the condition would generally deteriorate with time [6,7]. While the treatment of PI has constantly progressed, the above complications are ultimately irreversible. Hence, it is necessary to find candidate genes and molecular mechanisms involved in PI to obtain a better understanding of this condition.

Identifying biomarkers and gene networks associated with PI can elucidate the pathophysiological mechanisms of PI while simultaneously benefiting research and the development of new therapies. One study found that the level of transient receptor potential vanilloid 4 (TRPV4) decreased significantly in PI, suggesting that mechanical loading plays a key role [8]. In addition, Lin et al. highlighted the importance of the NF-κB signaling pathway in PI development and cartilage degeneration [9]. The activation of the PI3K/AKT signaling pathway was also considered to be associated with PI and trochlear dysplasia [10]. Besides, increased matrix metallopeptidase (MMP)-2 and MMP-13 expression, and decreased tissue inhibitor of metalloproteinase (TIMP)-2 expression were observed in PI, indicating that proteinases may be involved in PI [9,11,12]. The number of osteoclasts was also found to be increased in PI which may contribute to the significant subchondral bone loss [8,13]. However, no studies performed to date have been able to comprehensively available to elucidate the molecular mechanism behind PI. Understanding the etiological basis of PI will lead to better treatment of it.

Our main purpose of this study is to identify the potential crucial genes and molecular mechanisms in PI. We identified genes differentially expressed in PI by RNA-seq for the first time. Subsequently, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to identify the biological functions and mechanisms of action of the DEGs. Then, we established a protein-protein interaction (PPI) network and selected the top 10 potential hub genes with the highest calculated degrees. We also identified the possible miRNAs and transcription factors (TFs) involved in the regulation of DEGs using NetworkAnalyst and the ENCODE database, respectively. Finally, we explored the potential therapeutic drugs or molecular compounds for PI in the DGIdb database. Our study was the first to reveal molecular mechanisms underlying PI. Taken together, our findings provide promising biomarkers of PI and suggest a molecular basis for therapeutic targets of this condition.

Methods

Animal models

Ethical approval for this study was obtained from the ethics committee of Hebei Medical University Third Affiliated Hospital (H2019-021-1). Three-week-old male C57BL/6 mice were purchased from Vitalriver (Beijing, China). After acclimation for 1 week, the mice were randomly divided into two groups (control group and PI group, n = 7 for each group). The protocols of modeling were performed with reference to the previous reports [9,14]. The mice were anesthetized with sevoflurane. To establish the PI model, we cut the skin and subcutaneous tissues of the right knees sequentially to expose the patella and its retinaculum. Then, a 3-µm incision was made on the medial retinaculum to induce lateral patellar dislocation. In the control group mice, only the skin and subcutaneous tissues were cut. Finally, the incisions were sutured layer by layer with 6–0 absorbable sutures. Postoperative antibiotics were administered in the first 3 days after surgery. The distal femoral samples were collected 4 weeks after the surgery [13]. The samples were stored in RNAlater stabilization solution to prevent RNA degradation.

Gross observation

Immediately after sample collection, we assessed and recorded the morphology of the femoral trochlea to evaluate the rate of successful modeling.

RNA isolation

Total RNA was extracted using Eastep Super Total RNA Extraction Kit (Promega, Shanghai, China). The concentration, purity, and integrity of the total RNA were checked using a NanoPhotometer spectrophotometer (IMPLEN, CA, USA) and Agilent Bioanalyzer 2100 (Agilent, CA, USA). The RNA with RNA integrity number (RIN)≥ 8.0 was subjected to RNA-seq.

RNA-seq analysis

RNA-seq analysis was performed by Genechem Co., Ltd (Shanghai, China). Sequencing libraries were generated using NEBNext Ultra RNA Library Prep Kit (NEB, MA, USA). Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads and then fragmented into short pieces. The fragmented mRNA was used as the template to synthesize the cDNA. The cDNA was purified with the AMPure XP system (Beckman Coulter, Beverly, CA, USA) to preferentially select 250 ~ 300-bp-long cDNA fragments. Then PCR amplification was performed to construct a cDNA library. The quality of the library was checked using Agilent Bioanalyzer 2100 (Agilent, CA, USA). Finally, the library preparations were sequenced on an Illumina NovaSeq platform (Illumina, CA, USA).

Clean data (clean reads) were obtained by removing reads containing adapters, reads containing poly-N, and low-quality reads from the raw data. Then, Q20, Q30, and GC content of the clean data were calculated. All downstream analyses were based on clean data of high quality.

DEG identification

The qualified clean reads were aligned to the mouse reference genome using Hisat2 (version 2.0.5). The fragments per kilobase per million reads (FPKM) values were calculated for gene quantification. Differential expression analysis of the control and PI groups was performed using DESeq2 (version 1.16.1), and genes with padj< 0.05 and |log2(fold change)| > 1 were identified as DEGs. A volcano plot and a heatmap were generated with the ‘Ggpolt2’ and ‘ComplexHeatmap’ R packages [15].

GO and KEGG pathway analysis of DEGs

To characterize their biological functions, all DEGs were assessed using the GO and KEGG databases. GO terms were divided into three subgroups: biological process (BP), molecular function (MF), and cellular component (CC). GO terms and KEGG pathways with corrected p-value < 0.05 were considered significantly enriched. The GO function and KEGG pathway enrichment analysis were performed with ‘clusterProfiler’ R package [16], and the results were visualized with ‘Ggpolt2’ R package.

PPI network construction and hub gene identification

The PPI networks were predicted with the STRING database (http://string‐db.org; version 11.5) [17]. The PPI networks of DEGs were constructed with a combined score >0.4, and the network was visualized with Cytoscape (version 3.8.0). CytoHubba, a Cytoscape plugin [18], was employed to identify the hub genes. The top 10 genes with the highest calculated degrees were considered hub genes.

Target gene-miRNA network and the DEG-TF network analysis

The target gene-miRNA network was constructed with NetworkAnalyst (https://www.net workanalyst.ca/) [19], miRTarBase (http://mirtarbase.mbc.nctu.edu. tw/php/download.php), and TarBase (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=tarbase/ index), and the DEG-TF network was constructed with ENCODE (http://cistrome.org/BETA/). The target gene-miRNA networks and DEG-TF networks were visualized with Cytoscape (version 3.8.0).

Identification of potential drugs

The DEG-drug interaction network was constructed with the Drug Gene Interaction Database (DGIdb, version 3.0.2) (https://www.dgidb.org) [20]. We identified the potential drugs or molecular compounds interacting with the DEGs. The DEG-drug interaction network was visualized with Cytoscape (version 3.8.0).

RT-qPCR

The six most differentially expressed genes associated with PI were confirmed by real-time quantitative polymerase chain reaction (RT-qPCR). The sequences of the primers are listed in Table 1. Total RNA was reverse transcribed to cDNA using the GoScript Reverse Transcription System (Promega, Shanghai, China). The RT-qPCR was carried out in the Bio-Rad CFX96 system (Bio-Rad, CA, USA). The amplification procedure was as follows: one cycle of 2 min at 95°C, followed by 40 cycles of 15s at 95°C and 1 min at 60°C. The dissociation curve was obtained after the last cycle. GAPDH was used as an endogenous reference, and the 2−∆∆Ct method was used for evaluating gene expression. The significance of differences between the two groups was analyzed by Student’s t-test in SPSS 21.0 (SPSS Inc., IL, USA). All RT-qPCR experiments were conducted in triplicate.

Table 1.

Primer sequences used for RT-qPCR

Gene	Forward primer (5’-3’)	Reverse primer (5’-3’)
Il2ra	CAAGAACGGCACCATCCTAAA	TCCTAAGCAACGCATATAGACCA
Bach2	GAGGAAGGAGTTCCGAGCC	CAAGTCATCTTTCGTCTGTCCA
Slc12a3	GCCTTTGATGGACGGCAAG	GGATCACTCCCCAGATGTTGA
Fcrla	GATGATGGCGATATGACCCAAT	GCAGAACCAATGTGTCTCCTTC
Lgr5	GGACCAGATGCGATACCGC	CAGAGGCGATGTAGGAGACTG
Spib	AGGAGTCTTCTACGACCTGGA	GAAGGCTTCATAGGGAGCGAT
Gapdh	AGGTCGGTGTGAACGGATTTG	GGGGTCGTTGATGGCAACA

Immunohistochemistry (IHC)

Samples were decalcified in decalcifying solution for 2 weeks at 4°C. After decalcification, the samples were embedded in paraffin and then 4-μm sections were cut. Sections were de-waxed and rehydrated. Then, 3% hydrogen peroxide was applied to block endogenous peroxidase and protease K was applied to repair the antigen. Antibodies to Janus kinase (JAK)-1 (ab68153, 1:250; Abcam, UK) and signal transducer and activator of transcription (STAT)-3 (PA5-105,265, 1:200; Thermol Biotech Inc., USA) were added and incubated overnight at 4°C. The secondary antibodies were added and incubated for 30 min at 37°C, followed by immunochemical staining with 3, 3′ -diaminobenzidine. All sections were imaged with a microscope for analysis.

Results

The main aim of this study is to identify the potential crucial genes and molecular mechanisms involved in PI. RNA-seq was performed to assess gene expression. Integrative bioinformatic analysis was further conducted to provide the promising biomarkers of PI and suggest a molecular basis for therapeutic targets.

Four weeks after modeling, compared with the control group, the PI group exhibited a wider and shallower femoral trochlear (Figure 1). The results suggested that modeling was successful in six mice, with typical trochlear dysplasia, while all mice in control group had a normal femoral trochlear. One mouse in the PI group was excluded from further RNA-seq analysis because of unsuccessful modeling. Thus, a total of 13 samples (control group, n = 7; PI group, n = 6) were subjected to RNA-seq.

Figure 1.

Gross observation of the distal femoral samples in the two groups. (a) The control group; (b) The PI group.

Quality control of sequencing data

The results of quality control of sequencing data are shown in Table 2. After removal of data with poor quality, a total of 140.95 G clean data (with each sample over 10 G) were obtained for subsequent bioinformatics analyses. Both Q20 and Q30 scores of each sample were not less than 92.92%, with a sequencing error rate of less than 0.03%, indicating the high quality of our data. There was no GC bias.

Table 2.

Sequencing data quality

Sample	Raw reads	Clean reads	Clean bases	Error rate	Q20	Q30	GC content
Con 1	79,345,142	77,526,744	11.63 G	0.03%	97.93	94.19	52.34%
Con 2	73,958,028	72,354,576	10.85 G	0.02%	98.11	94.58	51.47%
Con 3	77,056,798	74,722,612	11.21 G	0.02%	98.03	94.50	52.76%
Con 4	71,612,042	69,983,152	10.50 G	0.03%	97.97	94.23	51.97%
Con 5	73,180,524	71,819,996	10.77 G	0.02%	98.01	94.33	52.77%
Con 6	72,080,330	70,472,942	10.57 G	0.03%	97.77	93.88	53.15%
Con 7	78,749,130	74,878,522	11.23 G	0.02%	98.06	94.57	53.56%
PI 1	68,682,320	66,894,460	10.03 G	0.02%	98.05	94.42	52.88%
PI 2	76,387,772	73,764,176	11.06 G	0.02%	98.02	94.50	53.88%
PI 3	73,639,656	69,062,156	10.36 G	0.02%	98.08	94.77	53.67%
PI 4	77,106,690	73,830,052	11.07 G	0.03%	97.93	94.23	52.87%
PI 5	70,332,966	67,591,882	10.14 G	0.02%	98.38	95.32	52.71%
PI 6	80,039,618	76,894,804	11.53 G	0.03%	97.42	92.92	52.04%

To reveal the molecular mechanism of PI, we conducted a DEG analysis to characterize the gene expression changes. Compared with levels in the control group, a total of 69 genes were differentially expressed in PI, of which 17 were upregulated and 52 downregulated (Figure 2 and Table 3). This result suggested that the great changes in gene expression in PI.

Figure 2.

Identifying the differentially expressed genes (DEGs) in our RNA-seq results. (a) The heatmap of DEGs; (b) The Volcano plot of DEGs.

Table 3.

DEGs in patellar instability

Gene name	log2FoldChange	p-value	padj
Up-regulated
Ighv14-4	1.230411708	2.21E-08	1.75E-05
Il11	3.15581519	3.18E-08	2.23E-05
Igkv4-68	1.156193746	1.10E-07	6.20E-05
Igkv6-14	1.70074259	1.96E-07	9.45E-05
Ighv5-9	1.382471754	4.40E-07	0.000190273
Khdc3	1.120919152	6.13E-07	0.000224567
Igkv6-25	1.181494972	9.57E-07	0.000316368
Igkv12-46	1.500363062	1.34E-06	0.000423206
Ighv1-53	1.033289659	1.36E-06	0.000423206
Ighv1-42	1.241004374	2.63E-06	0.000764885
Bpifb1	1.30171014	8.49E-06	0.002056468
Ighv4-1	1.479502639	3.26E-05	0.005384786
Ighg3	1.349365682	7.25E-05	0.010271636
Gpr55	1.126567421	0.000200518	0.021529524
Igkv4-57-1	1.399956803	0.000343995	0.032037175
Crlf1	1.173779701	0.000401287	0.036360478
Gm43442	1.422967796	0.000518718	0.043891742
Down-regulated
Il2ra	−1.581367181	4.37E-15	7.37E-11
Bach2	−1.092881649	3.91E-14	3.29E-10
Slc12a3	−1.185185393	1.64E-10	6.90E-07
Fcrla	−1.189134462	2.57E-10	7.38E-07
Lgr5	−1.233461068	2.63E-10	7.38E-07
Spib	−1.099262086	1.21E-09	2.55E-06
Casr	−2.293439462	1.60E-09	2.70E-06
Prg4	−1.331600891	5.66E-09	7.75E-06
Myl4	−1.219919296	5.98E-09	7.75E-06
Il7r	−1.107849771	1.65E-08	1.63E-05
Fam129c	−1.039459885	1.73E-08	1.63E-05
Gm30211	−1.191591246	1.83E-08	1.63E-05
Hist1h1a	−2.723512734	2.28E-08	1.75E-05
Hist1h1e	−1.543447586	1.21E-07	6.57E-05
4930426D05Rik	−1.109105307	1.42E-07	7.46E-05
Egfl6	−1.587267739	2.13E-07	9.97E-05
Hist1h2bj	−1.967690895	4.91E-07	0.000196914
Gm6525	−1.424661873	6.92E-07	0.000248155
Srpk3	−1.003529329	7.86E-07	0.00027023
Gm37065	−1.141423284	9.95E-07	0.000322469
Gm34095	−1.279253348	2.42E-06	0.000716262
Tmem45b	−2.611299235	2.86E-06	0.000817499
Pde4c	−1.247782838	4.47E-06	0.00125534
Hist1h1b	−1.576664829	4.95E-06	0.001344644
Gm18724	−1.150426594	5.91E-06	0.001556446
Gm38043	−2.770043626	7.76E-06	0.001926103
Capsl	−1.416787652	7.77E-06	0.001926103
Cytl1	−1.281203189	1.02E-05	0.002290215
Gpr12	−1.252992389	2.38E-05	0.004357949
Bach2os	−1.575327943	2.71E-05	0.004706088
Hist1h2ak	−2.07268382	3.04E-05	0.005175726
Gm1604a	−3.807290415	5.02E-05	0.007561664
Gm44891	−2.235769942	5.48E-05	0.008167655
5830487J09Rik	−1.294588163	7.93E-05	0.011051243
Plin1	−1.313514007	0.000102095	0.013010704
Adrb3	−2.562673414	0.00010659	0.013309555
Hist1h1d	−1.158848745	0.000121713	0.014760491
9630013D21Rik	−1.082347616	0.000154682	0.01749983
A530030E21Rik	−1.014603041	0.000186908	0.020459164
Hist1h4n	−1.273688948	0.000189121	0.020567866
Mcpt4	−1.70730339	0.000228382	0.024061472
Gm20743	−1.242422266	0.00023075	0.024160005
Igkj4	−1.431817534	0.000238191	0.024633019
Cidec	−1.493687242	0.00026973	0.027556571
Lrrc10b	−1.511437959	0.000312523	0.03010402
Gm45745	−1.105342772	0.000323989	0.030855793
2310008N11Rik	−1.990365136	0.000367681	0.033868862
Gm47730	−1.166361777	0.000405515	0.036360478
Cyp3a13	−1.416891958	0.000515551	0.043891742
Gria1	−1.237426726	0.000569684	0.046315443
Cilp	−1.053794767	0.00057149	0.046315443
Gm38120	−1.810672744	0.000628689	0.049292116

GO and KEGG pathway enrichment analyses were performed to characterize the biological functions of DEGs (Table 4). The GO terms were selected based upon p-value rankings when more than three terms were enriched for a given category. The GO analysis results showed that 8 GO terms were significantly enriched, with 3 related to BP, 2 related to CC, and 3 related to MF (Figure 3 and Table 5). BP terms included recognition of phagocytosis, complement activation, classical pathway, and positive regulation of lymphocyte activation. CC terms included circulating immunoglobulin complex and immunoglobulin complex. MF terms included immunoglobulin receptor binding, antigen binding and adrenergic receptor binding. Three pathways that were notably enriched among the DEGs were hematopoietic cell lineage, JAK-STAT signaling pathway and regulation of lipolysis in adipocytes (Figure 4 and Table 5).

Table 4.

GO and KEGG pathways analysis of DEGs

Ontology	ID	Description	Count	p-value	q-value
BP	GO:0006910	phagocytosis, recognition	6	5.18621E-07	0.000128832
BP	GO:0006958	complement activation, classical pathway	6	8.72251E-07	0.000128832
BP	GO:0051251	positive regulation of lymphocyte activation	8	9.5367E-07	0.000128832
BP	GO:0002455	humoral immune response mediated by circulating immunoglobulin	6	1.44908E-06	0.000128832
BP	GO:0050853	B cell receptor signaling pathway	6	1.49655E-06	0.000128832
BP	GO:0006956	complement activation	6	1.8088E-06	0.000128832
BP	GO:0006911	phagocytosis, engulfment	6	1.92398E-06	0.000128832
BP	GO:0002696	positive regulation of leukocyte activation	8	2.40207E-06	0.000128832
BP	GO:0099024	plasma membrane invagination	6	2.51903E-06	0.000128832
BP	GO:0072376	protein activation cascade	6	2.5935E-06	0.000128832
BP	GO:0050867	positive regulation of cell activation	8	3.06788E-06	0.000128832
BP	GO:0010324	membrane invagination	6	3.07899E-06	0.000128832
BP	GO:0050871	positive regulation of B cell activation	6	4.05128E-06	0.000156474
BP	GO:0002377	immunoglobulin production	6	1.00483E-05	0.00036038
BP	GO:0016064	immunoglobulin mediated immune response	6	1.17408E-05	0.000368661
BP	GO:0019724	B cell mediated immunity	6	1.2807E-05	0.000368661
BP	GO:0050864	regulation of B cell activation	6	1.30855E-05	0.000368661
BP	GO:0042113	B cell activation	7	1.39209E-05	0.000368661
BP	GO:0008037	cell recognition	6	1.39504E-05	0.000368661
BP	GO:0002449	lymphocyte mediated immunity	7	1.89855E-05	0.000476637
BP	GO:0002460	adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains	7	2.40587E-05	0.00055868
BP	GO:0050851	antigen receptor-mediated signaling pathway	6	2.44788E-05	0.00055868
BP	GO:0002757	immune response-activating signal transduction	7	2.7373E-05	0.00059757
BP	GO:0002764	immune response-regulating signaling pathway	7	3.28184E-05	0.000686596
BP	GO:0002429	immune response-activating cell surface receptor signaling pathway	6	4.27909E-05	0.000859422
BP	GO:0006909	phagocytosis	6	5.06117E-05	0.000972694
BP	GO:0002768	immune response-regulating cell surface receptor signaling pathway	6	5.23052E-05	0.000972694
BP	GO:0002440	production of molecular mediator of immune response	6	5.85939E-05	0.001050725
BP	GO:0006959	humoral immune response	6	6.86141E-05	0.001187982
BP	GO:0042742	defense response to bacterium	6	0.000190505	0.00318845
BP	GO:0007190	activation of adenylate cyclase activity	2	0.000787832	0.012760479
BP	GO:0042311	vasodilation	2	0.002319227	0.036077214
BP	GO:0002683	negative regulation of immune system process	5	0.002371113	0.036077214
BP	GO:0032781	positive regulation of ATPase activity	2	0.002581454	0.038122399
BP	GO:1902476	chloride transmembrane transport	2	0.002717568	0.03898587
BP	GO:0003044	regulation of systemic arterial blood pressure mediated by a chemical signal	2	0.005337652	0.074446204
BP	GO:0070231	T cell apoptotic process	2	0.005919162	0.080325474
BP	GO:0098661	inorganic anion transmembrane transport	2	0.006119164	0.080854322
CC	GO:0042571	immunoglobulin complex, circulating	6	3.29027E-07	1.70125E-05
CC	GO:0019814	immunoglobulin complex	6	3.89442E-07	1.70125E-05
MF	GO:0034987	immunoglobulin receptor binding	6	1.54235E-07	1.28259E-05
MF	GO:0003823	antigen binding	6	1.24389E-06	5.17198E-05
MF	GO:0031690	adrenergic receptor binding	2	0.000591143	0.012836482
MF	GO:0004896	cytokine receptor activity	3	0.000617451	0.012836482
KEGG	mmu04640	Hematopoietic cell lineage	3	0.000320573	0.014172688
KEGG	mmu04630	JAK-STAT signaling pathway	3	0.001695364	0.03747646
KEGG	mmu04923	Regulation of lipolysis in adipocytes	2	0.002894906	0.042661768

BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 3.

GO term enrichment analysis results. (a) The GO term bubble plot of DEGs; (b) The GO term chord plot of DEGs.

Table 5.

The gene lists of selected terms

Ontology	ID	Description	Genes
BP	GO:0006910	phagocytosis, recognition	Ighg3, Ighv5-9, Ighv14-4, Ighv1-42, Ighv4-1, Ighv1-53
BP	GO:0006958	complement activation, classical pathway	Ighg3, Ighv5-9, Ighv14-4, Ighv1-42, Ighv4-1, Ighv1-53
BP	GO:0051251	positive regulation of lymphocyte activation	Il2ra, Il7r, Ighg3, Ighv5-9, Ighv14-4, Ighv1-42, Ighv4-1, Ighv1-53
CC	GO:0042571	immunoglobulin complex, circulating	Ighg3, Ighv5-9, Ighv14-4, Ighv1-42, Ighv4-1, Ighv1-53
CC	GO:0019814	immunoglobulin complex	Ighg3, Ighv5-9, Ighv14-4, Ighv1-42, Ighv4-1, Ighv1-53
MF	GO:0034987	immunoglobulin receptor binding	Ighg3, Ighv5-9, Ighv14-4, Ighv1-42, Ighv4-1, Ighv1-53
MF	GO:0003823	antigen binding	Ighg3, Ighv5-9, Ighv14-4, Ighv1-42, Ighv4-1, Ighv1-53
MF	GO:0031690	adrenergic receptor binding	Adrb3, Gria1
KEGG	mmu04640	Hematopoietic cell lineage	Il11, Il2ra, Il7r
KEGG	mmu04630	JAK-STAT signaling pathway	Il11, Il2ra, Il7r
KEGG	mmu04923	Regulation of lipolysis in adipocytes	Adrb3, Plin1

Figure 4.

KEGG pathway enrichment results. (a) The KEGG pathway bubble plot of DEGs; (b) The KEGG pathway chord plot of DEGs.

The STRING database was used to construct a PPI network of the DEGs, and the results revealed 64 nodes and 156 edges (Figure 5a). Then, Cytohubba plugin [3] was used to screen out the top 10 hub genes associated with PI. These hub genes were Jak1, Jak3, Il2, Stat5a, Il7, Il7r, Il2ra, Il2rb, Il11, and Tslp (Figure 5b and Table 6).

Figure 5.

PI-specific network. (a) The protein–protein interaction network of DEGs constructed with the STRING database; (b) The hub genes with the top 10 degree (Red indicates a higher degree, and yellow indicates a lower degree).

Table 6.

Degree of top 10 hub genes

Rank	Gene ID	Gene name	Degree
1	Jak1	Janus kinase 1	16
2	Jak3	Janus kinase 3	15
3	Il2	interleukin 2	14
4	Stat5a	signal transducer and activator of transcription 5A	13
4	Il7	interleukin 7	13
4	Il7r	interleukin 7 receptor	13
7	Il2ra	interleukin 2 receptor, alpha chain	12
8	Il2rb	interleukin 2 receptor, beta chain	11
8	Il11	interleukin 11	11
8	Tslp	thymic stromal lymphopoietin	11

Target gene-miRNA network and DEG-TF network analysis

To explore the regulation of DEGs at the post-transcriptional stage, we conducted DEG-miRNA network and DEG-TF network analysis. The top three DEGs for miRNAs were Bach2, which could be regulated by 112 miRNAs; Hist1h1d, which could be regulated by 31 miRNAs; and Il2ra, which could be regulated by 31 miRNAs. The miRNA that could modulate the largest number of DEGs (13 genes) was mmu-mir-495-3p (Figure 6a). The top three DEGs for TFs were Hist1h1d, which could be regulated by 25 TFs; Hist1h1a, which could be regulated by 20 TFs; and Fam129c, which could be regulated by 19 TFs (Figure 6b).

Figure 6.

PI-specific network. (a) The network of target gene-miRNA; (b) The network of DEGs-TF. The yellow circle nodes represent the genes, and blue diamond nodes represent the miRNAs/TFs.

To explore the potential therapeutic drugs or molecular compounds of PI, a drug–DEG interaction network was constructed with DGIdb. Finally, as shown in Figure 7, a variety of drugs or molecular compounds were shown to be able to modulate the expression of 12 DEGs (IL2RA, SLC12A3, LGR5, CASR, IL7R, PDE4C, GPR12, PLIN1, ADRB3, GRIA1, IL11, and GPR55). Some drugs or molecular compounds were found to interact with multiple DEGs. For instance, cyclothiazide was found to regulate SLC12A3 and GRIA1.

Figure 7.

The network of drug-gene interaction. The yellow circle nodes represent the genes, and blue diamond nodes represent the drugs.

Confirmation of sequencing results by RT-qPCR

The mRNA levels of the six most differentially expressed genes associated with PI were confirmed by RT-qPCR (Figure 8a). The RNA-seq and RT-qPCR results were compared and their correlations were determined by Pearson correlation coefficient. The results of RT-qPCR were consistent with the RNA-seq results (Figure 8b, r= 0.954), suggesting the reliability of the sequencing results.

Figure 8.

Validation of RNA-seq results by RT-qPCR.

IHC staining of JAK1 and STAT3

The expression of representative genes of the JAK/STAT signaling pathway was determined by IHC staining. Two genes, JAK1 and STAT3, were selected, which also served as hub genes in PI. Both JAK1 and STAT3 were found to be elevated in PI (Figure 9). Overall, the expression of JAK1 and STAT3 as detected by IHC staining was consistent with the RNA-seq based integrative analysis.

Figure 9.

IHC comparison of JAK1 and STAT3 between control and PI groups. All samples were observed under microscopy at ×10 magnification.

Discussion

PI is a prevalent, debilitating musculoskeletal disorder that is often accompanied by refractory pathological changes. In the early stage, loss of bone mass, especially in subchondral bone, and morphological abnormality of femoral trochlea appear. With the progression of PI, cartilage degeneration is aggravated and ultimately results in patellofemoral osteoarthritis (PFOA) [21]. Aberrant gene expression is clearly closely related to many pathological processes in PI [9]. However, the crucial driver genes and molecular mechanisms associated with the condition and its complications remain unclear.

This study performed integrated bioinformatic analysis of changes in the expression of crucial genes to reveal potential pathways and gene interactions involved after the development of PI based on RNA-seq results. In total, we identified 69 DEGs in PI, including 17 upregulated and 52 down-regulated ones.

BP annotation revealed that the DEGs of PI were significantly enriched in the recognition of phagocytosis, complement activation, classical pathway and positive regulation of lymphocyte activation, which are all related to immune responses. Osteoarthritis (OA) is low-grade sterile inflammation, with innate immunity playing an essential role [22,23]. The genes identified here suggested that immune response activation acts as an essential trigger in PI and PFOA.

KEGG enrichment analysis showed that the JAK/STAT signaling pathway and regulation of lipolysis in adipocytes were important pathways in PI. JAK/STAT signaling pathway activation would contribute to bone loss, promote cartilage degeneration, and accelerate joint inflammation and destruction [24-26]. Thus, JAK/STAT signaling pathway may strongly correlate with PI development. Obesity is a population-based risk factor for PI and OA [27,28]. Aberrant lipid metabolism, especially increased adipokines, can cause OA even in non-weight bearing joints [29]. Our results provide new insights regarding the molecular mechanisms for PI.

Ten PI-related hub genes, namely, Jak1, Jak3, Il2, Stat5a, Il7, Il7r, Il2ra, Il2rb, Il11, and Tslp, were identified in this study all of which are involved in the JAK/STAT signaling pathway or inflammation-related pathways. Because PI was characterized by the disruption of bone homeostasis and the onset of PFOA, these genes were consistent with the pathological features of PI. Besides, the JAK/STAT signaling pathway was also enriched in KEGG analysis. Accumulating studies have confirmed that JAK/STAT inhibitors could reduce bone loss in various circumstances, prevent cartilage from degenerating, and even improve cartilage repair [30-32]. Therefore, targeting the JAK/STAT signaling pathway may ameliorate PI and PI-related complications. As mentioned before, inflammation-related genes can be considered excellent PFOA biomarkers. These inflammation-related hub genes would help to better understand the development and progression of PFOA.

In our study, the top three DEGs in the miRNA-gene network were Bach2, Hist1h1d, and Il2ra. Clinical and animal studies have shown that Il2ra is associated with joint destruction and arthritis persistence [33,34]. miR-495-3p, which regulates the greatest number of DEGs, exerts anti-inflammatory effects in different ways in multiple diseases [35-38]. Besides, miR-495-3p also regulates cartilage degeneration in intervertebral disc endplates [39]. However, there is a need for further validation of the function of miR-495-3p in PI. The top three DEGs for the TF-gene network were Hist1h1d, Hist1h1a, and Fam129c. Hist1h1d is a specific gene that can be regulated by most miRNAs or TFs, suggesting that it may serve as an important node in PI and could be a promising target for treatment.

To explore possible effective treatments for PI and its complications, we used the DGIdb database to determine therapeutic drugs or molecular compounds that could restore dysregulated expression of DEGs. Cyclothiazide, an inhibitor of DL-alpha-amino-3-hydroxy-5-methylisoxasole-4-propionate (AMPA) receptor desensitization, has been confirmed to regulate the proliferation and/or differentiation of osteoblasts and chondrocytes [40-42]. Of the drugs mentioned in Figure 8, further studies are urgently needed to evaluate their capacity for clinical application.

Our study has some limitations. Further studies should verify whether the genes identified in our study are involved in PI and explore their potential mechanisms. Besides, given that we only performed transcriptomic analysis of PI, it would be helpful if other omics analysis could be performed. In addition, the sample sizes for RNA-seq were relatively small. Thus, studies with larger-scale sequencing are necessary for further validation. Finally, our RNA-seq was performed in mice. Although the overall homology between mice and humans is 99%, there may be subtle differences between human and murine genetic changes. In vitro studies of human cells or human tissues should be performed to explore their roles in shaping the development of PI.

To the best of our knowledge, this is the first study using RNA-seq technology to identify potential DEGs in PI. We also revealed the PI-related pathways, predicted possible regulatory mechanisms, and explored promising drugs for PI. Thus, our study provided potential targets for future research.

Conclusion

In summary, we first compared the expression of all genes between PI and control mice by RNA-seq and found significant changes of various genes involved in PI development. Subsequently, we performed integrative bioinformatic analysis to identify molecular mechanisms and hub genes in PI. We also revealed the potential pathways, predicted miRNAs or TFs regulating DEGs, and explored promising drugs for PI. Further studies are needed to confirm our findings in order to uncover the exact mechanisms behind PI and develop novel non-surgical treatments for PI and PI-related complications.