Screening and verification of hub genes involved in osteoarthritis using bioinformatics.

Osteoarthritis (OA) is one of the most common causes of disability and its development is associated with numerous factors. A major challenge in the treatment of OA is the lack of early diagnosis. In the present study, a bioinformatics method was employed to filter key genes that may be responsible for the pathogenesis of OA. From the Gene Expression Omnibus database, the datasets GSE55457, GSE12021 and GSE55325 were downloaded, which comprised 59 samples. Of these, 30 samples were from patients diagnosed with osteoarthritis and 29 were normal. Differentially expressed genes (DEGs) were obtained by downloading and analyzing the original data using bioinformatics. The Gene Ontology enrichment and Kyoto Encyclopedia of Genes and Genomes pathways were analyzed using the Database for Annotation, Visualization and Integrated Discovery online database. Protein-protein interaction network analysis was performed using the Search Tool for the Retrieval of Interacting Genes/proteins online database. BSCL2 lipid droplet biogenesis associated, seipin, FOS-like 2, activator protein-1 transcription factor subunit (FOSL2), cyclin-dependent kinase inhibitor 1A (CDKN1A) and kinectin 1 (KTN1) genes were identified as key genes by using Cytoscape software. Functional enrichment revealed that the DEGs were mainly accumulated in the ErbB, MAPK and PI3K-Akt pathways. Reverse transcription-quantitative PCR analysis confirmed a significant reduction in the expression levels of FOSL2, CDKN1A and KTN1 in OA samples. These genes have the potential to become novel diagnostic and therapeutic targets for OA. Copyright: © Xie et al.

Introduction

Osteoarthritis (OA) is a chronic degenerative disease of the joints, characterized by cartilage degeneration, chronic inflammation and decreased normal joint function (1). Its clinical manifestations are mainly local pain and limited activity of the joint (2). OA commonly occurs in females over 55 and males over 65 years of age (3). With economic improvements and the aging of the population, an increasing number of patients are diagnosed with OA, which is now the most common cause of disability worldwide. The number of patients with OA is estimated to be >47 million in the US and is forecast to increase to ~67 million in 2030 (25% of the adult population) (4). The quality of life of patients with OA is significantly lower than that of others and poses great physical, psychological and financial burdens.

Although the pathogenesis of OA has received much clinical attention, genetic factors associated with the development of this disease remain elusive (5). Previously, OA was thought to be a heritable disease, but later studies have not determined any clear hereditary factors in the pathogenesis of OA (6). However, bioinformatics has advanced the study of OA and has identified specific genes involved in disease severity (7). Similar advances have been made in the study of rheumatoid arthritis (8,9). While hub genes involved in OA have been identified as key genes in its pathogenesis, no comprehensive genetic analysis has yet been performed (10). Several studies have performed gene expression profiling of OA samples and screened thousands of differentially expressed genes (DEGs) using high-throughput sequencing technology and other advanced techniques (11,12). However, a comprehensive analysis of all gene data collected is still missing. An integrated bioinformatics approach is able to predict and identify the hub genes involved in ОA.

In the present study, three datasets, GSE12021(13), GSE55457(14) and GSE55235(11), were analyzed. Screening of DEGs between OA patients and normal controls was performed. Gene Ontology (GO) and DEG pathway enrichment, protein-protein interaction (PPI) network and functional module analyses were then performed to explore the underlying molecular mechanisms of the pathogenesis of OA.

Materials and methods

Microarray data

National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO) is a public online repository for high-throughput gene queries and high-throughput gene expression detection for the global research community (15). It was used in the present study to obtain OA-related genes and their expression values were downloaded for further analysis. The GSE12021, GSE55235 and GSE55457 gene expression profiles were downloaded from the GEO database. Subsequently, three datasets from OA patients and normal controls (GSE55235, GSE55457 and GSE77298) were used for analysis; however, the data were downloaded without OA stage identification.

Screening of DEGs by GEO2R

The GEO2R online database (https://www.ncbi.nlm.nih.gov/gds/) was used to identify DEGs between OA samples and normal controls from three databases (GSE55235, GSE55457 and GSE77298); |log fold change (FC)|>1 and P<0.01 were considered to be DEGs. Venn diagrams were used to screen for common significant differences among DEGs.

GO enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis

For GO analysis, functional terms enriched by the DEGs were determined in three distinct categories representing different biological aspects: Molecular function (MF), biological process (BP) and cellular component (CC) (16). KEGG pathways were determined to analyze gene functions and link the genetic information in the genome with the biological functional information of genes. The Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) was used to perform GO analysis and KEGG analysis.

Construction of PPI network

PPI network construction was performed using the Search Tool for the Retrieval of Interacting Genes/proteins (STRING) database (https://string-db.org/cgi/network.pl). Subsequently, the DEGs were imported into Cytoscape software (http://www.cytoscape.org/) for better visualization. Next, the MCODE plugin was used to screen for key DEGs among all DEGs.

Reverse transcription-quantitative (RT-q)PCR verification

RT-qPCR was used to verify the four key genes. The RNA samples obtained from chondrocytes were extracted and transcribed into complementary DNA by an RNAiso Plus and PrimeScript® RT reagent kit (Takara Bio, Inc.). The DyNAmo SYBR-Green and qPCR systems (Takara Bio, Inc.) were used to analyze gene expression. Fill the ice box with samples. A ribozyme free 200 µl EP tube was placed on ice, after which 1,000 ng RNA samples were added. A total of 1 µl gDNA Eraser and 1 µl 5X gDNA Eraser Buffer was subsequently added to samples. Finally, RNase Free ddH2O was added to samples making a total volume of 10 µl. Amplification 2 min at 42˚C. Then at 4˚C, 4 µl 5X PrimeScript Buffer 2, 4 µl RNase Free ddH2O, 1 µl 1PrimeScrip Enzyme Mix1 and 1 µl RT Prime Mix was added to samples. Amplification 15 min at 37˚C and 5 sec at 85˚C. A total of 20 µl cDNA was obtained by reverse transcription. Subsequently, 5 µl SYBR Premix Ex Taq II with 3.5 µl dd H2O, 0.2 µl ROX, 0.5 µl cDNA, 0.4 µl PCR Forward Primer (10 µM) and 0.4 µl PCR Reverse Primer (10 µM) was added 96-well plates. The plates were run on the fluorescence quantitative PCR instrument, using the following cycling conditions: 95˚C for 30 sec, followed by 40 cycles at 95˚C for 5 sec, 60˚C for 30 sec, 90˚C for 15 sec and 60˚C for 60 sec. The paired primers for the four key genes are listed in Table I. All samples were compared with GAPDH. The 2-ΔΔCq method was used to quantify the relative gene expression levels (17).

Table I

Primer sequences for four hub genes (5'-3').

Gene	Forward primer	Reverse primer
BSCL2	ATGGTCAACGACCCTCCAGTA	GCTGACTGTCGGCATATAGGAA
FOSL2	CAGAAATTCCGGGTAGATATGCC	GGTATGGGTTGGACATGGAGG
CDKN1A	TGTCCGTCAGAACCCATGC	AAAGTCGAAGTTCCATCGCTC
KTN1	AAATGTCTTCGTAGATGAACCCC	TTTGTCAGTTTCGGTCTTCAGTT
GAPDH	GGCACAGTCAAGGCTGAGAATG	ATGGTGGTGAAGACGCCAGTA

BSCL2, BSCL2 lipid droplet biogenesis associated, seipin; FOSL2, FOS-like 2, activator protein-1 transcription factor subunit; CDKN1A, cyclin-dependent kinase inhibitor 1A; KTN1, kinectin 1.

Cell culture and immunofluorescence assays

The collection and culturing of chondrocytes was performed as described in detail in previous studies by our group (18-20). Specimens of the knee joint and femoral head were collected from patients undergoing joint replacement surgery at the Bone and Joint Department of Shenzhen Second People's Hospital (Shenzhen, China) between December 2017 and October 2019. According to the patient's imageological diagnosis, all OA samples obtained were from patients with OA of grades III and IV in Kellgren-Lawrence classification (21). In addition, as the control, the samples of patients with femoral neck fracture on radiographic diagnosis were used. The average age of the patients was 78.44±3.42 years in the control group, including 9 patients (2 males and 7 females), and 66.60±4.53 years in the OA group, including 5 patients (5 females). The following procedures were performed under sterile conditions. The samples were washed with normal saline three times and a surgical blade was used to cut the cartilage tissue blocks to final sizes of approximately 1x1x1 mm. The samples were incubated with collagenase type II (1 mg/ml; cat. no. C6885; Sigma-Aldrich; Merck KGaA) working liquid and oscillated for digestion for 8-12 h at 37˚C. The cell suspension was then divided into two parts. One part was added to complete chondrocyte culture medium (containing 10% FBS; Gibco; Thermo Fisher Scientific, Inc.), 1% penicillin and streptomycin, 1% HEPES, 1% ascorbic acid, 1% proline solution and 1% non-essential amino acids) and cultured in a CO2 incubator. The other part was centrifuged and resuspended and cells were dropped onto a glass slide for subsequent experiments. The present study was approved by the Ethics Committee of Shenzhen Second People's Hospital (Shenzhen, China) and written informed consent was obtained from all subjects.

For the immunofluorescence experiment, methods similar to those of other studies were used (22,23). The steps of immunofluorescence labeling of target protein were as follows: After washing the cells with PBS two times, the cells were fixed with paraformaldehyde for 15 min. The cells were then washed with PBS six times for 5 min each. Triton X100 (0.5%, diluted with PBS) was added, and samples were incubated for 15 min at room temperature. Cells were again washed by PBS six times for 5 min each. The cells were then blocked with 5% bovine serum albumin (Amresco LLC) for 2 h at room temperature and then washed with 0.5% BSA six times for 5 min each. The samples were then incubated with primary antibody (rabbit anti-Collagen II/FITC conjugated antibody; cat. no. bs-10589R-FITC; BIOSS) at 4˚C overnight. As these conjugated fluorescent dye primary antibodies demonstrated a strong specificity, secondary antibodies were not required. Subsequently, the cells were washed with 0.5% BSA six times for 5 min each. Subsequently, the cells were washed by PBS six times for 5 min each. DAPI (Thermo Fisher Scientific, Inc.) was added for 5 min and cells were washed with PBS three times for 5 min each. Images of fluorescently labeled cells were acquired with an LSM800 confocal microscope (Zeiss AG).

Statistical analysis

Prism 8 (GraphPad Software, Inc.) was used to generate figures. Experimental data were statistically analyzed using SPSS 19.0 software (IBM Corp.). Student's t-test was used to assess differences between groups. P<0.05 was considered to indicate a statistically significant difference.

Results

Screening of DEGs

Genes with significant differences in expression (P<0.051 and |log FC|>1) between the OA samples and normal controls were screened in the GEO2R website and 181 DEGs, including 48 upregulated and 133 downregulated genes, were identified (Fig. 1).

Figure 1

Common differentially expressed genes from three datasets. By analyzing the differential genes of GSE12021, GSE55235 and GSE55457 datasets, 181 common differential genes were obtained.

GO functional enrichment analysis

The functions of the 181 DEGs were classified using DAVID. In the GO category MF, these DEGs were significantly enriched in protein-DNA binding, RNA binding and protein homodimerization activity. In addition, in the category CC, the DEGs were significantly enriched in the nucleoplasm, nucleus, cytoplasm and plasma membrane. However, in the category BP, there were fewer genes enriched compared with the other categories (Fig. 2 and Table II). The upregulated DEGs were mainly enriched in the extracellular space and plasma membrane. The downregulated DEGs were significantly enriched in the nucleoplasm, nucleus and cytoplasm, and were involved in protein binding.

Figure 2

GO analysis results of hub genes. GO analysis classified the differentially expressed genes into three groups: Molecular function, biological process and cellular component. GO, Gene Ontology.

Table II

Significant enrichment of DEGs in osteoarthritis.

A, Downregulated DEGs
Category	Term	Description	Gene count	P-value
BP	GO:0071850	Mitotic cell cycle arrest	2	8.41x10-2
BP	GO:0045821	Positive regulation of glycolytic process	2	9.02x10-2
BP	GO:0010839	Negative regulation of keratinocyte proliferation	2	9.02x10-2
BP	GO:0070498	Interleukin-1-mediated signaling pathway	2	9.02x10-2
BP	GO:0097009	Energy homeostasis	2	9.64x10-2
CC	GO:0005654	Nucleoplasm	39	2.45x10-6
CC	GO:0005634	Nucleus	56	5.71x10-5
CC	GO:0005737	Cytoplasm	48	4.78x10-3
CC	GO:0005667	Transcription factor complex	5	3.59x10-2
CC	GO:0000790	Nuclear chromatin	5	3.59x10-2
MF	GO:0005515	Protein binding	80	3.11x10-4
MF	GO:0042803	Protein homodimerization activity	15	4.81x10-4
MF	GO:0044822	poly(A) RNA binding	19	6.78x10-4
MF	GO:0000166	Nucleotide binding	9	2.78x10-3
MF	GO:0003677	DNA binding	22	4.73x10-3
B, Upregulated DEGs
Category	Term	Description	Gene count	P-value
BP	GO:0007155	Cell adhesion	6	4.97x10-3
BP	GO:0016525	Negative regulation of angiogenesis	3	1.00x10-2
BP	GO:0030336	Negative regulation of cell migration	3	2.25x10-3
BP	GO:0071257	Cellular response to electrical stimulus	2	2.89x10-2
BP	GO:0007399	Nervous system development	4	3.27x10-2
CC	GO:0005615	Extracellular space	11	6.34x10-4
CC	GO:0005578	Proteinaceous extracellular matrix	5	3.02x10-3
CC	GO:0005886	Plasma membrane	18	5.37x10-3
CC	GO:0005788	Endoplasmic reticulum lumen	4	9.15x10-3
CC	GO:0030133	Transport vesicle	3	1.93x10-2
MF	GO:0001786	Phosphatidylserine binding	3	2.97x10-3
MF	GO:0008201	Heparin binding	4	5.12x10-3

BP, biological process; CC, cellular component; MF, molecular function; GO, Gene Ontology; DEG, differentially expressed gene.

KEGG signaling pathway analysis

DEGs were mainly enriched in the MAPK and PI3K-Akt signaling pathways (Fig. 3) according to the KEGG analysis. The signaling pathways of downregulated DEGs were mainly enriched in osteoclast differentiation, MAPK signaling, TNF signaling, Epstein-Barr virus infection and insulin signaling. Upregulated genes were mainly enriched in influenza A signaling (Table III).

Figure 3

Kyoto Encyclopedia of Genes and Genomes analysis results of differentially expressed genes. The color gradient represents the P-value; the size of the dots represents the gene number.

Table III

Kyoto Encyclopedia of Genes and Genomes analysis result of DEGs.

A, Downregulated DEGs
Term	Description	Gene count	P-value
hsa04380	Osteoclast differentiation	7	9.17x10-4
hsa04010	MAPK signaling pathway	9	1.46x10-3
hsa04668	TNF signaling pathway	6	2.30x10-3
hsa05169	Epstein-Barr virus infection	6	4.05x10-3
hsa04910	Insulin signaling pathway	6	6.81x10-3
B, Upregulated DEGs
Term	Description	Gene count	P-value
hsa05164	Influenza A	3	8.21x10-3

DEG, differentially expressed gene; hsa, Homo sapiens.

PPI network builds and MCODE analysis

To further analyze the relationships among all DEGs, PPI networks were constructed, including 161 nodes and 375 edges. The network included 181 protein interactions with combined scores of >0.4 (Fig. 4A). By further analyzing the PPI networks, four modules were detected using the MCODE plugin in Cytoscape. According to the MCODE score, key genes in four modules were screened, including BSCL2 lipid droplet biogenesis associated, seipin (BSCL2), FOS-like 2, activator protein-1 transcription factor subunit (FOSL2), cyclin-dependent kinase inhibitor 1A (CDKN1A) and kinectin 1 (KTN1) (Fig. 4B-D).

Figure 4

PPI networks constructed by Search Tool for the Retrieval of Interacting Genes/proteins. (A) PPI network constructed with the upregulated and downregulated differentially expressed genes. Red nodes represent upregulated genes and blue nodes represent downregulated genes. The potential interaction networks of (B) cluster 1 containing 11 molecules and (C) cluster 2 containing 14 molecules and (D) cluster 3 containing 8 molecules identified from Molecular Complex Detection (MCODE). Three significant modules were obtained from the PPI network. PPI, protein-protein interaction.

Validation of key genes

To verify the results, the expression levels of four key genes were determined in human articular chondrocytes (Fig. 5). First, normal chondrocytes (Ctrl group) and osteoarthritis chondrocytes (OA group) were obtained from patients undergoing joint replacement and subsequently, cells were cultured and stained with toluidine blue (Fig. 5A and B). To further identify the two types of cells, immunofluorescence was used to observe differences in the expression of type II collagen, which was decreased significantly in OA chondrocytes (Fig. 5C). Subsequently, the expression levels of the four key genes were determined by RT-qPCR. The results indicated that the expression level of BSCL2 in OA samples was increased but with no significant difference (P>0.05, n=3; Fig. 5F), while the expression levels of FOSL2, CDKN1A and KTN1 were significantly decreased. The differences between OA and normal chondrocytes in the expression levels of the four key genes were consistent with the analytical results of PCR (P<0.05, n=3; Fig. 5D, E and G).

Figure 5

RT-qPCR results of four key genes. (A and B) Growth of normal chondrocytes (A) under a normal microscope and (B) staining with toluidine blue (passage number of the cells is P1; scale bar, 20 µm). (C) Difference in expression of collagen II (green) between normal chondrocytes and OA chondrocytes under confocal microscopy (scale bar, 20 µm). (D-G) Validation of the differential expression of four hub genes between the OA group and the normal group by RT-qPCR. (D) FOSL2, (E) CDKN1A, (F) BSCL2 and (G) KTN1. All samples were normalized to the expression of GAPDH and the relative expression levels of each gene were analyzed using the 2-ΔΔCq method. *P<0.05. OA, osteoarthritis; Ctrl, control; BSCL2, BSCL2 lipid droplet biogenesis associated, seipin; FOSL2, FOS-like 2, activator protein-1 transcription factor subunit; CDKN1A, cyclin-dependent kinase inhibitor 1A; KTN1, kinectin 1.

Discussion

OA is considered to be the most common cause of disability (24); it seriously affects the quality of life, creating physical, psychological and economic burdens for patients. The diagnosis and treatment of OA require an interdisciplinary approach. Bioinformatics is a tool widely used to predict potential target genes for numerous diseases. In the present study, three sets of data from patients with OA and normal controls were analyzed by using GEO2R. Between the OA and control groups, 181 DEGs were identified, including 48 upregulated and 133 downregulated genes. The upregulated DEGs were mainly enriched in the extracellular space and plasma membrane. Downregulated DEGs were significantly enriched in functional terms associated with the nucleoplasm, nucleus, cytoplasm and protein binding, and this explains that the main difference between patients with OA and individuals without OA is in the cartilage cells. Cellular changes are an important part of the process of OA. OA is a chronic degenerative disease characterized by cartilage defects and chondrocytes are the only cell type in articular cartilage (25). Cartilage cell apoptosis is a key part of OA (26). The present study suggested that the downregulated gene CDKN1A, which regulates the occurrence and development of apoptosis, may participate in the development of OA. Further analysis indicated that these DEGs were mainly enriched in the ErbB, MAPK and PI3K-Akt signaling pathways. These pathways are associated with chondrocyte apoptosis in OA. MAPK modulates joint inflammation and joint damage (27). The major pathological changes in OA are cartilage damage and poor healing (1), so the MAPK pathway is able to promote the aggravation of OA or rheumatoid arthritis to exert an important role in the pathogenesis of arthritis (26,28). Erb activates PI3K to generate phosphatidylinositol-3,4,5-triphosphate that effectively activates the Akt pathway. Apoptosis-associated proteins are competitively inhibited by the Akt pathway. Thus, recruiting this gene may promote cell survival and inhibit cell apoptosis (29). There is a definite correlation between the degree of cartilage damage and chondrocyte apoptosis (26). Therefore, inhibiting chondrocyte apoptosis may be effective in regulating cartilage degeneration during OA. Animal experiments have indicated that the use of MAPK inhibitors effectively improves inflammation and joint degeneration in mice with OA (30), which demonstrates the potential of inhibiting MAPK in OA treatment. The PI3K/Akt pathway is associated with TNF-α-induced activation of OA fibroblast-like synoviocytes, which may be involved in OA pathogenesis (31). It promotes cartilage degeneration, subchondral bone function damage and inflammation during OA (32). The use of MAPK inhibitors reduces the pathological changes of OA (31). CDKN1A is able to activate the AKT pathway and induce MAPK8 to participate in the inactivation of MAPK, further reducing apoptosis (33,34). It has been reported that CDKN1A is significantly downregulated in the synovium of arthritis patients and is associated with inhibition of chondrocyte proliferation (35), suggesting that upregulation of CDKN1A may have a positive role in the early treatment of OA.

By contrast, the expression of BSCL2 in OA samples was significantly increased compared with that in normal samples. BSCL2 has a key role in lipogenesis, lipid metabolism and lipid droplet synthesis (36). The loss of BSCL2 may lead to serious disorders in metabolic dysfunction and a significant reduction in fat (37), whereas upregulation of BSCL2 may lead to weight gain. Obesity is a major cause of OA (38), and thus, upregulation of BSCL2 may induce its development. Furthermore, various transforming (FOS) proteins affect the physiology of chondrocytes, osteoblasts and osteoclasts (39). FOSL2 is the key regulator of leptin expression in fat cells (40) and its deficiency is able to promote obesity. In the present study, its low expression in OA samples was noted, which is consistent with obesity being a leading cause of OA (38). Recent studies have indicated that FOSL2 is suppressed in the early hypertrophy state of chondrocytes, suggesting that this gene is strongly associated with the early initiation of OA (41). Therefore, overexpressed BSCL2 and decreased expression of FOSL2 in OA samples promotes obesity in patients with OA, and their regulation may have a positive role in the prevention of OA. Furthermore, KTN1 is a receptor on the endoplasmic reticulum (42) that has an important part in adjusting protein biosynthesis in cells (43). A proteomics analysis of human mesenchymal stem cells undergoing inhibited chondrogenesis indicated reduced levels of KTN1(44), suggesting an association of KTN1 with chondrogenesis enhancement. Thus, enhancing KTN1 may induce chondrocyte proliferation.

Of note, the present study had a limitation: When NCBI GEO was used to obtain the microarray data, the information on the OA stage was unclear or the diagnosis of the OA classification by different doctors was biased. Therefore, it is difficult to ensure that all patients whose data were downloaded had the same OA classification. In order to confirm certain key genes that regulate the occurrence of osteoarthritis identified in bioinformatics screenings, more data and further experimental verification are still required.

A limitation of the present study was that cells were prepared using collagenase II. Further experiments are therefore required to determine whether this may have affected results.

In conclusion, bioinformatics and experimental data suggested that BSCL2, FOSL2, CDKN1A and KTN1 are key DEGs in OA compared with normal samples. Therapeutic targeting of these genes may positively contribute to the treatment and prevention of OA and consequent disability. One limitation of the present study is the limited amount of collected data. Further analyses and experiments are required prior to performing some potential targets.