Combined bioinformatics analysis reveals gene expression and DNA methylation patterns in osteoarthritis.

Osteoarthritis (OA) is a common type of arthritis, which may cause pain and disability. Alterations in gene expression and DNA methylation have been proven to be associated with the development of OA. The aim of the present study was to identify potential therapeutic targets and associated processes for OA via the combined analysis of gene expression and DNA methylation datasets. The gene expression and DNA methylation profiles were obtained from the Gene Expression Omnibus, and differentially expressed genes (DEGs) and differentially methylated sites (DMSs) were identified in the present study, using R programming software. The enriched functions of DEGs and DMSs were obtained via the Database for Annotation, Visualization and Integrated Discovery. Finally, cross analysis of DEGs and DMSs was performed to identify genes that exhibited differential expression and methylation simultaneously. The protein‑protein interaction (PPI) network of overlaps between DEGs and DMSs was obtained using the Human Protein Reference Database; the topological properties of PPI network overlaps were additionally obtained. Hub genes in the PPI network were further confirmed via reverse transcription‑quantitative polymerase chain reaction (RT‑qPCR). The results of the present study revealed that the majority of DEGs and DMSs were upregulated and hypomethylated in patients with OA, respectively. DEGs and DMSs were primarily involved in inflammatory, immune and gene expression regulation‑associated processes and pathways. Cross analysis revealed 30 genes that exhibited differential expression and methylation in OA simultaneously. Topological analysis of the PPI network revealed that numerous genes, including G protein subunit α1 (GNAI1), runt related transcription factor 2 (RUNX2) and integrin subunit β2 (ITGB2), may be involved in the development of OA. Additionally, RT‑qPCR analysis of GNAI1, RUNX2 and ITGB2 provided further confirmation. Numerous known and novel therapeutic targets were obtained via network analysis. The results of the present study may be beneficial for the diagnosis and treatment of OA.

Introduction

Osteoarthritis (OA) is the most common type of arthritis and a leading cause of pain and disability, which places a great burden on the economy of health and reduces quality of life (1–3). OA involves the degeneration of numerous tissues, including subchondral bone, ligaments, muscle, tendons, and the meniscus and synovium (1). Numerous factors may affect OA progression, including age, gender, obesity, genetics and joint injury (4); however, how these factors affect the development of OA requires further investigation and no effective method has been developed for the relief of pain. In addition, molecular biology studies have identified numerous biomarkers and biological processes that contribute to OA, including the erosion of the extracellular matrix (5), the expression of chemokines (chemokine C-C ligands 9 and 5, and interleukin-8) (6) and the upregulation of inflammatory genes (7). Further investigation into the molecular events associated with cartilage degeneration is required.

Over the past decades, the development of high-throughput technologies has resulted in the large amount of accumulation of omics data for various complex diseases. For OA, gene expression profiling via microarray or high-throughput sequencing has become a promising method for the analysis of the mechanisms underlying its initiation and progression (8). For example, Rasheed et al (9) performed an integrated study of microRNA (miRNA) expression profiles in OA chondrocytes and OA-associated genes, and identified numerous miRNAs associated with the development of OA. Sun et al (10) reported several potential biomarkers for OA via differential expression and network analysis based on gene microarray datasets. Microarray analysis in the study of Loeser et al (11) indicated the link between age-associated differences in gene expression and the development of OA. In addition, epigenetic modifications serve important roles in gene expression regulation, and DNA methylation is one of the most common types of epigenetic modification. Recently, an increasing number of studies have focused on the associations between methylation status and the progression of OA (12,13). In contrast to cancer, in which CpG sites are frequently hypermethylated, the majority of studies investigating OA reported a higher frequency of hypomethylation (14,15). DNA methylation may also affect the allelic imbalance of specific small nucleotide polymorphisms, and thus the development of OA (16). Combined analysis of gene expression and DNA methylation profiles may contribute to the screening of potential biomarkers of OA, early diagnosis and treatment; to the best of our knowledge, an investigation into this is yet to be performed.

In the present study, combined analysis of publicly accessible gene expression and DNA methylation microarray datasets of OA was conducted. Functional enrichment and network analysis was performed for the identification of potential biomarkers. Numerous known and novel targets were obtained and their involvement in OA was further confirmed via reverse transcription-quantitative polymerase chain reaction (RT-qPCR) analysis.

Materials and methods

Microarray datasets

The publicly accessible data were all obtained from the Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo). Gene expression profiles deposited by Klinger et al (17), accession no. GSE43923, containing six samples (three osteophytic cartilage and three corresponding articular cartilage samples from the knee joints of patients with OA) were employed in the present study. The genome-wide expression profiles were quantified using the commercial gene microarray GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA). The DNA methylation profiles (GSE73626) (18) of five hip OA, six knee OA and seven hip healthy cartilage samples were detected via Illumina HumanMethylation450 BeadChip assay (Illumina, Inc., San Diego, CA, USA), which contains >480,000 methylation sites, covering 99% of RefSeq (https://www.ncbi.nlm.nih.gov/refseq/) genes and 96% of University of California, Santa Cruz (http://genome.ucsc.edu/)-defined CpG sites with an average of 17 CpG sites/gene across different genomic regions, including the promoter, 5′ untranslated region (UTR), first exon, gene body, intergenic and 3′UTR.

Microarray data analysis

The present study conducted differential expression analysis for osteophytic and articular cartilage samples from patients with OA. The raw CEL data were imported into R version 3.2.2 (http://www.R-project.org/) and normalized via the affy package (19); subsequently, the limma package (20) was used for the screening of differentially expressed genes (DEGs) with the criteria of fold change >1.5 and false discovery rate (FDR)<0.05. For the methylation dataset, site-level analysis was performed based on the Illumina Methylation Analyzer package (Illumina, Inc.) (21) to obtain the differentially methylated CpG sites (DMSs) between hip/knee OA cartilage and healthy cartilage samples, with thresholds of db value >0.2 and FDR<0.05. DMSs were mapped to the corresponding genes (DMGs) and genomic regions based on the full annotation file of the microarray and following this, cross analysis was performed via the ‘intersect’ function of R version 3.2.2 (http://www.R-project.org/) using DEGs and DMGs to reveal overlapping genes. In addition, differences between distributions of DMSs relative to CpG islands and genes were compared using the χ2 test.

Functional clustering analysis

Investigation into the functions of enriched DEGs and DMGs may improve understanding of their involvement in OA. In the present study, functional clustering analysis of DEGs and DMGs based on the Database for Annotation, Visualization and Integrated Discovery (DAVID; david.abcc.ncifcrf.gov) (22) was conducted. Clusters with an enrichment score >1, and Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.jp/kegg) pathways with P<0.05 were retained in the present study.

Protein-protein interaction network analysis

Genes are likely to function together rather than alone in complex diseases; hub nodes in the network may represent key biomarkers. In the present study, protein-protein interaction (PPI) network analysis was performed to investigate the overlaps between DEGs and DMGs based on the Human Protein Reference Database (HPRD; www.hprd.org) (23). The network was visualized using Cytoscape 3.6.0 software (http://www.cytoscape.org/), and the topological property of every gene [degree (number of direct interactions in the network)] was additionally analyzed for the assessment of their importance.

RT-qPCR

Normal and OA tissues were obtained from the articular cartilage of 26 females and 20 males with a median age of 56.35 years (43.58–69.12 years) between April 2012 and April 2015 in Zibo Central Hospital (Zibo, China). Patients exhibiting temporomandibular joint pain that were not suffering from any form of rheumatic disease or cancer were included in the present study. The study was approved by the ethics committee of Zibo Central Hospital. Written informed consent was obtained from all patients.

Total RNA was extracted from OA and normal tissues (50–100 mg) using an RNeasy Mini kit (Qiagen GmbH, Hilden, Germany) and quantified with a NanoDrop system (Thermo Fisher Scientific, Inc.), and subsequently subjected to RT-qPCR using EasyScript Reverse Transcriptase kit (Promega Corporation, Madison, WI, USA). The temperature protocol used for RT was as follows: 95°C for 10 min, 55°C for 1 min and 68°C for 10 min.

The 7500 Real-Time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.) was used for qPCR. Reactions were conducted in triplicate in each reaction tube using AceQ qPCR SYBR Green Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China). The temperature protocol used for qPCR was as follows: 94°C for 5 min; followed by 40 cycles of 95°C for 10 sec and 60°C for 30 sec; followed by 95°C for 15 sec, 60°C for 60 sec and 95°C for 15 sec. Data were analyzed via the 2−ΔΔCq method using GAPDH as internal control (24). Primer sequences were: G protein subunit α1 (GNAI1) forward, 5′-TTAGGGCTATGGGGAGGTTGA-3′, and reverse, 5′-GGTACTCTCGGGATCTGTTGAAA-3′; runt related transcription factor 2 (RUNX2) forward, 5′-TGGTTACTGTCATGGCGGGTA-3′ and reverse, 5′-TCTCAGATCGTTGAACCTTGCTA-3′; integrin subunit β2 (ITGB2) forward, 5′-TGCGTCCTCTCTCAGGAGTG-3′ and reverse, 5′-GGTCCATGATGTCGTCAGCC-3′; and GAPDH forward, 5′-ACAACTTTGGTATCGTGGAAGG-3′ and reverse, 5′-GCCATCACGCCACAGTTTC-3′.

Statistical analysis

R version 3.2.2 (http://www.R-project.org/) was used for all of the statistical analysis. The relative mRNA expression levels in the RT-qPCR analysis were presented as the mean ± standard deviation of the three replicates. RT-qPCR data were analyzed using a Student's t-test. P<0.05 was considered to indicate a statistically significant difference.

Results

Gene expression profile analysis

The raw microarray data was normalized and used for the following differential expression analysis. As a result, a total of 466 genes were detected to be DEGs in osteophytic cartilage samples compared with articular samples, which contained 49 downregulated and 417 upregulated genes (Fig. 1A). The two-way supervised clustering indicated notable differences between osteophytic and articular cartilage samples from patients with OA, with blue and red coloring indicating low and high expression levels, respectively (Fig. 1B). The full list of DEGs is provided in Table I.

Figure 1.

Gene expression profile analysis. (A) Distribution of upregulated and downregulated genes. Red plus signs and green triangles represent upregulated and downregulated genes in patients with osteoarthritis, respectively. Black circles are non-DEGs. (B) Two-way supervised clustering of DEGs and samples, with blue and red coloring indicating low and high expression levels, respectively. DEGs, differentially expressed genes.

Table I.

Previously identified biomarkers of osteoarthritis with the PMID of each corresponding reference.

Author, date	Gene	(Refs.)
Mabey T et al, 2014	Angiopoietin 2	(43,44)
Gao W et al, 2013
Valverde-Franco G, 2016	Ephrin B2	(39)
Leijten JCH, 2013	Gremlin 1 DAN family BMP antagonist	(37,45)
Yi J, 2016
Fan Y, 2017	Integrin subunit β2	(46,47)
Hopwood B, 2007
Xiao JL, 2015	Runt related transcription factor 2	(35,48)
Liao L, 2017

DMSs

Comparisons were performed between DMSs in hip/knee cartilage samples and healthy cartilage samples. As presented in Fig. 2A, the β values of OA hip compared with healthy hip tissue, and OA knee compared with healthy knee tissue, of all CpG sites in the microarray were obtained. The β values for OA hip compared with healthy hip tissue, and OA knee compared with healthy knee tissue, of CpG sites satisfied the criterion of P<0.05; the frequencies of hypomethylated sites were increased compared with hypermethylated sites in OA hip and knee samples (Fig. 2B). Additionally, the number of hypomethylated sites with P<0.05 and db>0.2 were increased compared with hypermethylated sites (Fig. 2C). This was consistent with the results of the DEG analysis (the number of upregulated genes was increased compared with downregulated genes), as hypermethylation of the promoter may result in the downregulation of the corresponding gene.

Figure 2.

Scatter plot of β values in the hip and knee cartilage of patients with OA compared with control samples. (A) β values of all the CpG sites in the Illumina HumanMethylation450 BeadChip assay of OA hip cartilage compared with control hip cartilage (left panel), and OA knee cartilage compared with control knee cartilage (right panel). (B) β values of CpG sites satisfying the criterion of P<0.05 of OA hip cartilage compared with control hip cartilage (left panel), and OA knee cartilage compared with control knee cartilage (right panel). (C) β values of CpG sites satisfying the criteria of P<0.05 and δβ>0.2 of OA hip cartilage compared with control hip cartilage (left panel), and OA knee cartilage compared with control knee cartilage (right panel). OA, osteoarthritis.

To improve understanding of the functional significance of DMSs, the functional locations of DMSs were investigated. As presented in Fig. 3A, the majority of the DMSs were reported to be in intergenic, gene body and promoter regions. Additionally, four neighborhood locations were defined in the Illumina HumanMethylation450 BeadChip assay: 31% CpG islands, 23% shores (0–2 kb from canonical CpG islands), 10% shelves (2–4 kb from canonical) and open sea (rest of the sequence). Consistent with the BeadChip assay, the majority of the hypo- and hypermethylated sites in hip and knee cartilage samples were detected in open sea, and following the north and south shores (upstream and downstream shores). To investigate the associations between functional and neighborhood locations with differential methylation status, a χ2 test for data presented in Fig. 3. The results indicated P<5×10−7 in all of the cases, demonstrating that functional and neighborhood locations are associated with differential methylation status. One-way clustering of DMSs of hip and knee cartilage samples is presented in Fig. 4.

Figure 3.

Functional genomic and neighborhood location distributions of differentially methylated sites. (A) Functional genomic distributions of hypomethylation and hypermethylation CpG sites in OA hip and knee cartilage. (B) Neighborhood location distributions of hypomethylation and hypermethylation CpG sites in OA hip and knee cartilage. OA, osteoarthritis; UTR, untranslated region.

Figure 4.

Heatmap of differentially methylated sites in osteoarthritic hip cartilage and knee cartilage samples. (A) Heatmap of DMSs in osteoarthritic hip cartilage. (B) Heatmap of DMSs in knee cartilage samples. DMSs, differentially methylated CpG sites.

Functional clusters

Functional clustering analysis in DAVID resulted in three functional clusters for DEGs and DMGs. DEGs were primarily involved in the GO terms and KEGG pathways that were associated with ‘cell structure’, ‘inflammatory and immune response’, ‘substance synthetic’ and ‘metabolic’. ‘Guanosine 5′-triphosphaate (GTP)ase activity’, ‘gene expression regulation’ and ‘inflammatory and immune response’ were reported to be significantly enriched in DMGs (data not shown).

Network analysis

Network topological properties are important representations of their roles in specific biological processes and diseases. In the present study, 30 overlaps were obtained among DMGs of hip and knee cartilage samples and DEGs; 20 of these overlaps were reported to interact with other genes from the HPRD. PPI networks are presented in Fig. 5. GNAI1 directly interacted with 50 genes, which was markedly higher compared with the degree of the other 19 overlaps in the network, potentially indicating its important roles in the development of OA. Table I includes the five previously identified biomarkers of OA and their corresponding PMID nos.

Figure 5.

Protein-protein interaction network of the overlaps between differentially expressed genes and differentially methylated CpG sites of corresponding genes. Triangles represent overlapped genes and circles are other genes that interacted with these overlaps.

RT-qPCR analysis

A total of three hub genes, GNAI1, RUNX2 and integrin subunit β 2 (ITGB2), were subjected to RT-qPCR analysis for the quantification of their relative abundance in OA and control samples. The results of the RT-qPCR analysis were consistent with the results of the gene microarray analysis; the relative mRNA expression levels of GNAI1, RUNX2 and ITGB2 in OA samples were significantly increased compared with the control samples (Fig. 6).

Figure 6.

RT-qPCR analyses of GNAI1, RUNX2 and ITGB2. (A) Relative mRNA expression levels of GNAI1 in control and OA samples. (B) Relative mRNA expression levels of RUNX2 in control and OA samples. (C) Relative mRNA expression levels of ITGB2 in control and OA samples. The abundance of each of the three genes was quantified by RT-qPCR, and all were consistent with the gene microarray analysis. **P<0.01 vs. normal control. GNAI1, G protein subunit α1; ITGB2, integrin subunit β2; OA, osteoarthritis; RT-qPCR, reverse transcription-quantitative polymerase chain reaction; RUNX2, runt related transcription factor 2.

Discussion

OA is the most common degenerative disease of the synovial membrane, comprising the destruction and loss of articular tissues (25). The initiation and development of OA have been reported to be associated with numerous factors, including age, joint injury, obesity and chronic inflammation, in addition to genetic factors, including epigenetic modification and altered gene or miRNA expression (26). In-depth understanding of gene regulation in OA may contribute to diagnosis and treatment. For this purpose, DNA methylation and gene expression analysis of hip and knee cartilage of OA patients was performed, and potential targets for OA were identified and verified in the present study.

DNA methylation has been revealed to repress gene expression by blocking the sites at the promoter where transcription factors bind; hypermethylation of the promoter is associated with no or low transcription (27). In the present study, the majority of the DEGs were observed to be upregulated and DMSs were hypomethylated, consistent with the roles of DNA methylation. Unlike cancer, which is characterized by the hypermethylation of tumor suppressor genes, genome-wide hypomethylation in OA has been widely observed previously (28–30) and in the present study. It is known that OA is influenced by inflammatory chemokines (31), which was also demonstrated by the functional enrichment analysis of the present study. The hypomethylation and increased expression of a number of inflammation-associated genes have been observed to be associated with OA, including IL-8 (32), nuclear factor-κB (33) and pleckstrin (34), which may serve as therapeutic targets for OA. In the present study, numerous genes that were differentially expressed and methylated simultaneously were identified, including ITGB2, GNAI1 and RUNX2. To the best of our knowledge, the roles of the aforementioned genes in the development of OA have not been investigated; the RT-qPCR analysis conducted in the present study revealed their relative abundance in OA and adjacent tissues, which may indicate their association with the progression of OA. Furthermore, numerous genes were observed to directly interact with other genes in the PPI network, including RUNX2, bleomycin hydrolase and gremlin 1 DAN family BMP antagonist (GREM1). The majority of these genes have been demonstrated to be closely associated with the progression of OA; for example, RUNX2 polymorphisms may affect temporomandibular joint OA in females (35) and may influence the expression of other genes in OA (36). The mRNA expression levels of GREM1 are correlated with OA and may be regulated by OA-associated factors (37); in addition, GREM1 is a key regulator of synoviocyte hyperplasia and invasiveness (38). Valverde-Franco et al (39) reported that ephrin-B2 may be essential for normal bone growth, and its absence may lead to knee and hip OA. In addition, a number of genes with a high degree in the PPI network have not been proven to be associated with the development of OA. For example, GNAI1, also known as Gi, is a protein that can hydrolyze GTP and interact with other proteins, and T cell differentiation may alter its structure (40). Additionally, the altered expression of GNAI1 was observed to be associated with the progression of inflammatory and immune diseases (41,42), and may be considered to be a novel biomarker for OA.

Certain limitations of the present study were noted. Only three genes were analyzed via RT-qPCR and future studies are required to investigate more genes. Furthermore, osteophytes may be considered to contribute to regeneration in OA joints; however, further studies investigating the differences regarding whole genome expression between OA and healthy tissues should be employed for further analysis in the present study.

In conclusion, the present study provided a pipeline for the combined analysis of gene expression and DNA methylation datasets. In addition, numerous known and potential novel markers were proposed in the present study, which may contribute to diagnosis and treatment targets for OA; however, further investigation is required for confirmation of the functions exhibited by these markers.