Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Abstract

Objective

Synovial inflammation influences the progression of osteoarthritis (OA). Herein, we aimed to identify potential biomarkers and analyze transcriptional regulatory-immune mechanism of synovitis in OA using weighted gene coexpression network analysis (WGCNA).

Design

A data set of OA synovium samples (GSE55235) was analyzed based on WGCNA. The most significant module with OA was identified and function annotation of the module was performed, following which the hub genes of the module were identified using Pearson correlation and a protein-protein interaction network was constructed. A transcriptional regulatory network of hub genes was constructed using the TRRUST database. The immune cell infiltration of OA samples was evaluated using the single-sample Gene Set Enrichment Analysis (ssGSEA) method. The hub genes coexpressed in multiple tissues were then screened out using data sets of synovium, cartilage, chondrocyte, subchondral bone, and synovial fluid samples. Finally, transcriptional factors and coexpressed hub genes were validated via experiments.

Results

The turquoise module of GSE55235 was identified via WGCNA. Functional annotation analysis showed that “mineral absorption” and “FoxO signaling pathway” were mostly enriched in the module. JUN, EGR1, FOSB, and KLF4 acted as central nodes in protein-protein interaction network and transcription factors to connect several target genes. “Activated B cell,” “activated CD4T cell,” “eosinophil,” “neutrophil,” and “type 17 T helper cell” showed high immune infiltration, while FOSB, KLF6, and MYBL2 showed significant negative correlation with type 17 T helper cell.

Conclusions

Our results suggest that the expression level of apolipoprotein D (APOD) was correlated with OA. Furthermore, transcriptional regulatory-immune network was constructed, which may contribute to OA therapy.

Keywords

osteoarthritis weighted gene coexpression network analysis hub genes transcriptional regulatory-immune network multiple tissues

Introduction

Osteoarthritis (OA), the most common musculoskeletal disorder, can lead to disability and reduce the quality of life in the elderly.¹ To date, no available treatment has been effective in preventing OA disease progression, although joint replacement surgery remains an option.² A number of local and systematic risk factors, including sports injury, anatomical abnormity, obesity, genetics, aging, and low-grade systemic inflammation contribute to the initiation and development of OA.³ The pathology of OA involves multiple joint tissues, such as synovium, cartilage, subchondral bone, and synovial fluid. Therefore, OA is considered as a whole joint disease.⁴

Synovium inflammation contributes to the development of OA.⁵ The histological changes in OA synovium include synovial lining hyperplasia, sublining fibrosis, increased angiogenesis, and immune cells infiltration.⁶ Production of proinflammatory cytokines and infiltration of immune cells can aggravate the progression of OA.^5,7 Therefore, the OA synovium may release mediators that could be potential biomarkers and therapeutic targets for OA.^8,9

In the era of big data, transformation of biomedical big data into valuable knowledge has been one of the most important challenges in bioinformatics.¹⁰ Weighted gene coexpression network analysis (WGCNA) can identify biologically significant modules and genes associated with clinical traits.¹¹ WGCNA has been applied in varieties of data, including gene expression profiles, genetic marker data, proteomics data, metabolomic data, and other high-dimensional data.^12,13

In this study, we analyzed synovitis samples from patients with OA using the WGCNA method to identify significant modules, along with function annotation of module genes and hub genes. Then, a protein-protein interaction (PPI) network and a transcriptional regulatory-immune network of hub genes were constructed; hub genes coexpressed in synovium, cartilage, chondrocyte, subchondral bone, and synovial fluid tissues were then screened out. Finally, transcriptional factors and coexpressed hub genes were validated via quantitative reverse transcription PCR (RT-qPCR) and enzyme-linked immunosorbent assay (ELISA).

Materials and Methods

Data Collection from Patients with OA and Control Donors

Data sets were downloaded from the National Center for Biotechnology Information Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/), including GSE55235, GSE55457, GSE82107, GSE89408, GSE1919, GSE114007, GSE57218, GSE110606, and GSE51588. GSE55235, GSE55457, and GSE1919 data sets include OA synovial tissue samples, normal synovial tissue samples, and rheumatoid arthritis (RA) synovial tissue samples. GSE89408 includes OA, normal, RA, and other lesions synovial tissue samples. GSE82107 includes OA and normal synovial tissue samples. OA and normal synovial tissue samples were used for further research.

GSE114007, GSE129147, GSE57218, and GSE110606 data sets include OA and normal cartilage or chondrocyte samples. Another data set of cartilage samples, “Ann Rheum Dis,” was downloaded from the supplementary materials of an original article.¹⁴ GSE51588 includes OA and normal subchondral bone samples. Each sample has lateral and medial tibial plateaus. We divided this group into two data sets and named them as “OA-ND” and “MT-LT.” The OA-ND data set includes OA and normal subchondral bone samples. The MT-LT data set includes OA medial samples and relatively normal lateral samples from OA patients. The proteomics data set of synovial fluid, “Arthritis Rheum” was integrated from proteomic analysis of synovial fluid.¹⁵ The details for all data sets are shown in Supplementary Material S1.

Data Preprocessing and Identification of Differentially Expressed Genes (DEGs)

The GEO data set files were downloaded and processed by the GEOquery package using R language (version 3.6.3, http://r-project.org/).¹⁶ The probes corresponding to multiple genes were deleted. Considering the existence of one gene corresponding to multiple probes, the average expression of multiple probes was taken as the expression value of the gene. Differentially expressed genes (DEGs) were analyzed using “limma” or “edgeR” package for all data sets, according to the gene expression data.^17,18 The cutoff criteria of DEGs was |fold change (FC)|>2 and P < 0.05. Data filtering for differentially proteomic analysis of synovial fluid was performed using |FC|≥1.2 and P < 0.05. The results of GSE110606, “Ann Rheum Dis,” and “Arthritis Rheum” were directly downloaded from Supplementary Materials.^14,15,19 The details for all data sets preprocessing are shown in Supplementary Material S2.

Construction of Coexpression Modules Associated with OA Based on WGCNA

The R package WGCNA was used to construct a coexpression network. First, samples of synovium data set GSE55235 were clustered to detect and eliminate outliers using the flashClust function.²⁰ The independence and average connectivity degrees of the coexpression networks were mainly determined by a power value.²¹ Then, suitable soft thresholding power value was selected using the pickSoftThreshold function of the WGCNA algorithm if the independence degree was >0.85.²² Then, a hierarchical clustering tree was generated with the power value with a module minimum size cutoff of 30.²³ Genes with similar expression data were clustered and formed a branch of the tree, representing a gene module.

Identification of Clinically Significant Modules and Functional Annotation

In this process, gene significance (GS), module eigengene (ME), and module membership were calculated. GS was defined as the absolute value of the correlation between a gene and a trait. ME was the first principal component of the given module. Module membership was defined as the correlation between ME and the gene expression profile. The significance between ME and a clinical trait was measured using Pearson correlation. We selected the module of highest correlation coefficient with OA trait.¹¹ Functional annotation analysis of module genes, including gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes pathways (KEGG), was conducted using the online tool DAVID6.8 (https://david.ncifcrf.gov/). GO terms contain biological process, cell component, and molecular function. P < 0.05 was considered as the threshold, and the top seven significant results were visualized using R.

Identification of Hub Genes and Construction of PPI Network

Hub genes are defined as genes that have high connectivity inside coexpression modules. The hub genes of the significant module were determined through an absolute value of the gene module membership >0.7. Furthermore, all hub genes were mapped into the Search Tool for the Retrieval of Interacting Genes (STRING) database (v11.0; https://string-db.org/). A combined score of >0.4 was considered significant. The PPI network was visualized using Cytoscape software (version 3.7.1), and the top 20 hub genes were calculated using the cytoHubba plugin through the maximal clique centrality method.²⁴

Construction of Transcriptional Regulatory Network of Hub Genes and Calculation of Immune Cell Infiltration Score in OA

First, we downloaded 8,472 regulatory relationships between human transcription factors and target genes from TRRUST (v2) database to construct a background transcriptional regulatory network, and all hub genes of synovium OA data set GSE 55235 were mapped into this network.²⁵ We extracted the one-step neighborhoods of the hub nodes to build a transcriptional regulatory network of hub genes. A list of transcription factors was obtained from AnimalTFDB, including 1,557 known human transcription factors.²⁶ Cytoscape software was used to network visualization in which hub nodes, transcription factors, and target genes were marked in different ways in the network. Then, we got a metagene that contained 782 genes presenting 28 immune cell types (Suppl. Material S3).²⁷ The immune cell infiltration score of 28 immune cell types of GSE55235 synovium samples was determined using single-sample Gene Set Enrichment Analysis (ssGSEA) software implemented in the R GSVA package.^27,28

Construction of Transcription Factor and Immune Cell Network

The transcription factors in the hub genes transcriptional regulatory network were selected. We calculated the Pearson correlation coefficient between the expression of transcription factors and the infiltration score of immune cell types. The significant correlation between transcription factors and cell types was retained to construct a network (P < 0.01, |cor| > 0.6), which was displayed using Cytoscape visualization tool. The transcription factors and cells were displayed with different node shapes, and the correlation between them was expressed using different-colored lines.

Screening of Coexpressed Hub Genes in Multiple Tissues

DEGs and differentially expressed proteins of synovium, cartilage, chondrocyte, subchondral bone, and synovial fluid tissues were selected for hub genes validation. First, we obtained intersections of synovium GEO data sets, including GSE55457, GSE82107, GSE89408, and GSE1919. Second, intersections of cartilage and chondrocyte samples were presented, including GSE114007, GSE57218, “Ann Rheum Dis,” GSE129147, and GSE110606. Moreover, we integrated intersections of GSE51588 (OA-ND), GSE51588(MT-LT), and “Arthritis Rheum.” Data analysis and visualization were produced using an online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

Validation of Transcription Factors and Coexpressed Hub Genes with RT-qPCR and ELISA Experiments

To further validate the expression of transcription factors and coexpressed hub genes, articular cartilage, synovium, subchondral bone, and synovial fluid of OA samples were obtained at the time of total joint arthroplasty (age 53-70 years, Kellgren-Lawrence 4 grade). Control samples of synovial fluid were collected from patients receiving arthroscopic surgery, including meniscectomy and free body removal (age 42-54 years, Kellgren-Lawrence 0-1 grade). This study was approved by the institutional review board, and informed consent was obtained from each donor.

Isolation of fibroblast-like synoviocyte (FLS) and chondrocyte was conducted as previously reported.^29,30 Passage 8 of FLS and passage 2 of chondrocyte were used for further experiments.³¹ FLS and chondrocytes were seeded in 6-well plates with or without stimulation with 10-ng/mL interleukin-1β (IL-1β) for 48 hours. Medium supernatant, subchondral bone samples, and synovial fluid samples were preserved in liquid nitrogen until being utilized. RNA isolation, reverse transcription cDNA synthesis, and quantitative PCR (RT-qPCR) were undertaken following protocols (ESscience Biotech, Shanghai, China; Bio-Rad, Hercules, CA). The primer sequences of transcription factors and coexpressed hub genes were listed in Table 1. The data were analyzed using 2-ΔΔCq method and GAPDH served as internal controls. ELISA was performed to confirm the protein expression of coexpressed hub genes in medium supernatant and synovial fluid samples (Elabscience Biotech, Wuhan, China; Tecan, Switzerland). All experiments repeated at least three times.

Table 1.

Primers Used for RT-qPCR Amplification.

Gene	Primers	Sequence 5-3
JUN	Forward	TGAGTGACCGCGACTTTTCA
JUN	Reverse	TTTCTCTAAGAGCGCACGCA
MYC	Forward	CCACAGCAAACCTCCTCACA
MYC	Reverse	TCCAACTTGACCCTCTTGGC
FOSB	Forward	GCGCCGGGAACGAAATAAAC
FOSB	Reverse	CAACTGATCTGTCTCCGCCT
EGR1	Forward	GGTCAGTGGCCTAGTGAGC
EGR1	Reverse	GTGCCGCTGAGTAAATGGGA
ATF3	Forward	GCCATTGGAGAGCTGTCTT
ATF3	Reverse	GGGCCATCTGGAACATAAG
APOD	Forward	CCAC CCCA GTTAACCT CACA
APOD	Reverse	CCACTGTTTCTGGAGGGAGA
GAPDH	Forward	CACTCAGACCCCCACCACAC
GAPDH	Reverse	GATACATGACAAGGTGCGG CT

Statistical Analysis

The results are presented as mean ± the standard error of the mean (SEM). Statistical analysis was performed using GraphPad Prism 7.0 software (GraphPad Software Inc., San Diego, CA). The differences between groups were determined using the Student t test. P < 0.05 was considered statistically significant.

Results

Results of Data Processing

An overall flow chart of the study design is presented in Figure 1 . The synovium data set GSE55235 of OA was examined using DEG analysis, WGCNA analysis, function annotation, identification of hub genes, PPI interaction, and immune infiltration. The hub genes were screened with data sets of multiple tissues, while transcription factors and coexpressed hub genes were validated with RT-qPCR and ELISA experiments. A total of 919 genes were differentially expressed in the data set of GSE55235, including 454 upregulated genes and 465 downregulated genes (Suppl. Material S4). Findings related to the DEGs are presented as a heatmap and volcano plot ( Fig. 2 ). For GSE110606 chondrocyte samples, DEGs were directly downloaded from the supplementary materials of an original article.¹⁹ For the GSE51588 data set, DEGs were selected from OA-ND and MT-LT subchondral bone samples. For proteomics analysis, a final list of 63 proteins that were differentially expressed in late OA patients were identified, excluding the proteins that were inconsistent across the two experiment sets.¹⁵ For all data sets, the DEGs and the proteomics analysis details are shown in Supplementary Material S2.

Figure 1.

Flow diagram of the study design. The synovium data set GSE55235 of osteoarthritis was analyzed using differentially expressed genes analysis, weighted gene coexpression network analysis, function annotation, identification of hub genes, protein-protein interaction network analysis, and step-by-step immune infiltration. Then, the hub genes were screened with data sets of multiple tissues. Finally, transcriptional factors and coexpressed hub genes in multiple tissues were validated with experiments. WGCNA = weighted gene coexpression network analysis; PPI = protein-protein interaction.

Figure 2.

Screening of differentially expressed genes (DEGs) from the OA synovium data set GSE55235. (A) Heatmap of DEGs between OA and normal tissues. (B) Volcano plot of DEGs between OA and normal tissues. |Fold Change| > 2 and P < 0.05 were considered as cutoff criteria. OA = osteoarthritis.

Construction of Modules Associated with OA Based on WGCNA

For sample clustering analysis, no outlier samples were found using the flashClust function in the synovium data set GSE55235 ( Fig. 3A ). Gene coexpression modules were established for the expression of DEGs using the WGCNA method. An experienced soft threshold value of 6, which was the lowest power value for an independence degree of up to 0.85, was selected to construct a scale-free network ( Fig. 3B and C). The hierarchical cluster dendrogram included seven modules, excluding the gray module ( Fig. 3D ). Different colors represented different modules; the grey module was a set of genes that could not be aggregated into other modules. The range of genes, from 246 in the turquoise module to 60 in the black module, and the detailed number of genes in each module are shown in Supplementary Material S5.

Figure 3.

Construction of modules based on weighted gene coexpression network analysis. (A) Sample clustering of GSE55235 to detect outliers. (B) and (C) Soft power value screening of differentially expression genes for OA related coexpression modules. (D) Establishment of the cluster dendrogram in OA coexpression modules. OA = osteoarthritis.

Identification of the Most Significant Modules Associated with OA

A heatmap was drawn based on the interactive relationships of seven modules with the OA traits. The turquoise module was significantly correlated to the disease status (Cor = 0.95, P = 2e−10) ( Fig. 4A ). Module and OA status clustering showed that the turquoise module was the most correlated with OA trait ( Fig. 4B ). Analysis of the GS modules also showed that the turquoise module was correlated with the OA trait ( Fig. 4C ). Next, a correlation was identified between GS and module membership within modules, as shown in Figure 4D (Cor = 0.91 and P < 0.001). A highly significant correlation was observed between GS and module membership in this module.

Figure 4.

The combination of heatmap (A) and module clustering (B) indicates that the turquoise module was highly related to OA. (C) Analysis of GS across modules showed that the turquoise module was significant. (D) A scatterplot of GS for OA versus module membership in the turquoise module. OA = osteoarthritis; GS, gene significance.

Functional Annotation Results of Module Genes

To further explore the functions of turquoise module, the online tool DAVID was used to investigate the GO terms and KEGG pathways. The biological processes of the turquoise module were mainly enriched in cellular response to zinc ion, negative regulation of growth, cellular response to interleukin-1, nucleus, transcription factor activity, and RNA polymerase II core promoter proximal region sequence-specific binding ( Fig. 5A-C ). The KEGG pathway analysis revealed that the turquoise module was mostly enriched in mineral absorption, FoxO signaling pathway, MAPK signaling pathway, tumor necrosis factor (TNF) signaling pathway, osteoclast differentiation, and cytokine-cytokine receptor interaction ( Fig. 5D ). Metallothioneins were involved in the mineral absorption pathway (Suppl. Material S6).

Figure 5.

Function annotation of the turquoise module genes showing (A) biological process enrichment results, (B) cell component enrichment results, (C) molecular function enrichment results, and (D) Kyoto Encyclopedia of Genes and Genomes enrichment results. The size of the bubble indicates the gene number, while the color indicates the enrichment significance (−log10[P value]). cAMP = cyclic advenosine monophosphato; TNF = tumor necrosis factor.

Identification of Hub Genes and Construction of the PPI Network

Genes with high connectivity inside coexpression modules were considered as hub genes. A total of 143 hub genes in the turquoise module were identified with a module membership of >0.7. All hub genes are shown in Table 2 . To further investigate the PPIs among these hub genes, a PPI network was constructed using the STRING database. A total of 251 interactions were included with a combined score of >0.4 ( Fig. 6A ). The top 20 genes based on the maximal clique centrality method included JUN, ATF3, DUSP1, EGR, CDKN1A, FOSB, and MYC ( Fig. 6B ).

Table 2.

Hub Genes of the Turquoise Module.

Genes	KME	Genes	KME	Genes	KME	Genes	KME
ACKR3	0.798811	DEPP1	0.755117	LTC4S	−0.89995	SAV1	0.881266
ADAMTS1	0.904704	DNAJB1	0.823873	MAFF	0.913598	SCML1	0.814431
ADCY7	−0.71436	DOCK2	−0.81717	MAP3K8	0.834823	SEL1L3	−0.75612
Apolipoprotein D	0.741071	DPYS	−0.74154	MCL1	0.897231	SLC16A7	0.898654
APOLD1	0.902364	DUOX2	0.709594	MSX1	0.780386	SLC38A1	0.745017
ARL4A	0.762658	DUSP1	0.764781	MT1E	0.783728	SLC39A14	0.798964
ATF3	0.871799	EDNRB	0.780818	MT1F	0.777147	SLC5A2	−0.70182
ATXN1	−0.73369	EGR1	0.889559	MT1G	0.837243	SNRPA1	0.769283
BAG3	0.889848	EMP1	0.789792	MT1X	0.727683	SOCS3	0.763069
BCL6	0.9686	ENC1	−0.78235	MYC	0.910212	SPEN	0.861098
BIN3-IT1	0.792184	ERO1B	0.710164	NAGK	−0.92199	SRSF7	0.85942
BNIP3	0.774912	FCER1A	−0.79403	NKX3−2	−0.76917	SYTL2	−0.7913
BTG1	0.946997	FNDC4	0.783447	NPAS2	0.82257	TAPBPL	−0.74701
BTG2	0.924282	FOSB	0.905201	NR4A2	0.866944	TGIF1	0.754771
BTG3	0.951607	FOSL2	0.884781	OGFRL1	0.803917	THAP9−AS1	0.874231
BTN3A2	−0.82328	GADD45B	0.842	PELI1	0.766379	TIPARP	0.881802
C1RL	0.824866	GLT8D2	−0.87364	PHC2	0.74955	TMCO3	−0.82199
CBR1	−0.79062	GTPBP4	0.7844	PIGA	0.905694	NEDD9	0.799223
CCDC59	0.906281	HAS1	0.777987	PLA2G2E	0.750245	TNS3	−0.7699
CCL4	0.752835	HNRNPA1	0.856242	PLEKHA2	−0.80238	TPPP3	−0.85197
CCL8	0.755761	ICAM3	−0.73487	PLGRKT	−0.76878	TRAC	−0.78495
CCNL1	0.922631	IER2	0.791182	PLK2	0.720737	TRBC1	−0.74455
CCR5	−0.74034	IER5	0.85267	PMAIP1	0.871238	TRIB1	0.830818
CD2	−0.79831	IFRD1	0.855785	PNMA8A	−0.74178	TSKU	0.801417
CDC42EP4	0.821751	IGHG1	−0.82646	PNRC1	0.90628	TTLL1	−0.84309
CDH19	0.714417	IL11RA	−0.84117	POMT1	−0.76739	UAP1	0.85628
CDKN1A	0.908705	IL1B	0.875524	RAB23	−0.73901	VASH1	−0.75789
CEBPB	0.920827	JUN	0.751925	RAI1	−0.82005	WEE1	0.85174
CEBPD	0.906308	JUND	0.895574	RBPMS	0.803561	WWC1	−0.75226
CEP55	−0.71515	KCNMA1	−0.72659	REV3L	0.717088	ZBED8	−0.72634
CHD1	0.822113	KCNN4	−0.72303	RGS13	−0.70188	ZBTB43	0.745008
COQ10B	0.728188	KDM3A	0.904149	RHOB	0.936837	ZNF395	0.852841
CRY1	0.748101	KLF10	0.718001	RHOBTB1	−0.89084	ZNF671	−0.74945
CX3CR1	−0.82106	KLF2	0.886407	RLF	0.85195	ZNF74	−0.70301
CXCL3	0.813085	KLF4	0.959921	RND3	0.946211
DDX3Y	0.757571	LINC00312	0.78322	RORA	0.752516

Figure 6.

Protein-protein interaction network construction of hub genes in the turquoise module. (A) Overview of the network with a combined score of >0.4. The larger the size of the node, the greater the number of connected nodes. (B) The top 20 genes are calculated based on the maximal clique centrality method. The redder colors of the node indicate greater connectedness.

Transcriptional Regulation and Immune Characterization of Hub Genes

A total of 8,472 relationships between transcription factors and target genes were downloaded from TRRUST database, which contained 796 transcription factors and 2,492 target genes, while 90 hub genes were mapped to the TRRUST transcriptional regulatory networks. In addition, one-step neighbourhoods of hub genes were selected to construct a transcriptional regulatory network of hub genes. Finally, this subnetwork contained 3,110 edges and 649 nodes (Suppl. Material S7-S10). In the network, hub genes such as JUN, EGR1, JUND, KLF4, and MYC acted as central nodes and transcription factors to connect a large number of target genes. The biggest degree node, SP1, regulated by EGR1, also regulated the highest number of target genes in the network as a transcription factor (Fig. 7A). The immune cell infiltration scores of each sample were calculated using the ssGSEA method and was visualized using R function heatmap. Enrichment results showed that “activated B cell,” “activated CD4T cell,” “eosinophil,” “neutrophil,” and “type 17 T helper cell” showed high immune infiltration in most OA samples (Fig. 7B). The network of transcription factors and cell types was calculated using the Pearson correlation among transcription factors and cell types in the hub genes regulatory network. The final network contained 28 cell type nodes, 159 transcription factor nodes, and 417 edges. The relationship between transcription factors and cell types was divided into positive correlation and negative correlation. For example, CEBPB and FOSL1 transcription factors were positively correlated with immature B cell infiltration. In contrast, increase in the expression of SRF and KLF4 resulted in decrease in the infiltration of immature B cell ( Fig. 7C ). The top 10 significantly correlated transcription factors and cell types are shown in Figure 7D. FOSB, KLF6, and MYBL2 were significantly positively correlated with effector memory CD8T cell (P < 0.01). Type 17 T helper cell and FOSB, KLF6, and MYBL2 had significant negative correlation (P < 0.01, |cor| > 0.8) (Fig. 7D and E).

Figure 7.

Transcriptional regulation and immune characterization of hub genes. (A) Transcriptional regulatory network of hub genes. Yellow round nodes represent hub genes, blue round nodes represent target genes, and the square nodes with a red border represent transcription factors. (B) A heatmap of immune infiltration scores in OA samples. The color of the heatmap from blue to red indicates the increase in immune infiltration scores. (C) Transcription factor and immune cell type network. Yellow round nodes represent cell types, green triangle nodes represent transcription factors, yellow lines represent positive correlation, and grey lines represent negative correlation. (D) A correlation between the transcription factor and immune cell type. The color and proportion of pie chart indicate the Pearson correlation coefficient. (E) The correlation plot of “type 17 T helper cell” and expression of FOSB and KLF6. Each yellow point represents an OA sample. OA = osteoarthritis.

Screening of Coexpressed Hub Genes in Multiple Tissues

A total of 12 data sets were used for validation of hub genes, including synovium sample data sets of GSE55457, GSE82107, GSE89408, and GSE1919; cartilage/chondrocyte sample data sets of “Ann Rheum Dis,” GSE114007, GSE129147, GSE57218, and GSE110606; subchondral bone sample data sets of GSE51588 (OA-ND) and GSE51588 (MT-LT); and synovial fluid sample data set of “Arthritis Rheum.” A Venn diagram showed that APOD was the only gene in synovium, synovial fluid, and subchondral bone tissues (Fig. 8A and C). In addition, APOD and SERPINE2 were coexpressed in cartilage/chondrocyte samples ( Fig. 8B ). In conclusion, APOD was the only hub gene coexpressed and downregulated in multiple tissues found in the data sets analyzed. The different expression value of APOD in different data sets is shown in Table 3.

Figure 8.

Screening coexpressed genes with multiple tissue data sets. (A) Venn diagram reveals that APOD was the only gene in synovium interactions. (B) Two genes were obtained in cartilage and chondrocyte interactions. (C) APOD was the only gene in subchondral bone and synovial fluid interactions. APOD = apolipoprotein D.

Table 3.

Different Expression of Apolipoprotein D in All Data Sets.

Data Set	Log Fold Change	P Value
GSE55235	−4.4189	4.23E−19
GSE55457	−1.4999	0.008406
GSE82107	−1.3686	0.007558
GSE89408	−1.6691	1.43E−06
GSE1919	−4.3451	7.67E−21
GSE114007	−3.2731	2.17E−05
Ann Rheum Dis	−1.5733	2.36214E−13
GSE129147	−2.1458	0.001328
GSE57218	−3.3741	1.16E−07
GSE110606	−1.83	0.0033
GSE51588(OA–ND)	−1.5523	0.003996
GSE5158(MT–LT)	−1.4024	0.000627
Arthritis Rheum	−0.51	<0.05

Validation of Transcription Factors and Coexpressed Hub Genes

Transcription factors such as JUN, MYC, FOSB, EGR1, and ATF3 were validated using RT-qPCR. The expression of APOD was identified using RT-qPCR and ELISA in different tissues. We isolated RNA from relatively normal lateral tibia and sclerosis area of the medial tibia (Fig. 9E). Results showed that the mRNA expression of JUN, FOSB, and EGR1 was significantly decreased in IL-1β stimulated FLS (Fig. 9A). The mRNA expression of APOD was downregulated in FLS, chondrocyte, and subchondral bone samples with OA (Fig. 9B-D). The protein expression of APOD in medium supernatant of IL-1β stimulated FLS and chondrocyte and OA synovial fluid was also downregulated ( Fig. 9F-H ).

Figure 9.

Validation of transcription factors and APOD expression. (A) Relative mRNA expression levels of transcription factors in FLS. N = 3 for each sample group. (B-C) APOD was downregulated in FLS and chondrocyte stimulated by IL-1β. N = 3 for each sample group. (D) APOD was downregulated in lesioned medial tibia than relatively normal lateral tibia. N =10 and n = 9 for medial and lateral sample groups, respectively. (E) We isolated total ribonucleic acid from relatively normal lateral tibia and the sclerosis area of the medial tibia. (F-G) Relative protein expression of APOD in medium supernatant of FLS and chondrocyte stimulated by IL-1β. N = 3 for each sample group. (H) Relative protein expression of APOD in synovial fluid. N = 5 for both osteoarthritis samples and control samples. Data are shown as mean ± standard error. ns = not significant; IL-1β = interleukin-1β; APOD = apolipoprotein D; FLS = fibroblast-like synoviocyte; *P < .05, **P < .01, ***P < .001.

Discussion

Abnormal synovial changes can occur at different stages of OA.³² Inflammation of the synovial membrane can produce proinflammatory and catabolic mediators such as cytokines, prostaglandin E2, nitric oxide, and neuropeptides, which can accelerate the breakdown of cartilage.⁹ Indeed, an increasing number of studies have indicated that immune infiltration plays important roles in OA progression.³³ The synovial membranes of OA patients often exhibit immune infiltration, such as in macrophages, T cells, mast cells, B cells, plasma cells, natural killer cells, dendritic cells, and granulocytes.³⁴ Based on the WGCNA method for synovium tissue transcription analysis, the most significant module associated with the OA trait was identified and APOD was identified a potential biomarker for OA diagnosis. Furthermore, transcriptional factors including JUN, MYC, EGR1, JUND, KLF4, FOSB, and CEBPB were identified, and a transcriptional regulatory-immune network was constructed.

In the OA synovium data set, the turquoise module was the most significant and the correlation value was 0.95. The KEGG pathways of the module were mostly enriched in mineral absorption, FoxO signaling pathway, and MAPK signaling pathway. Metallothioneins (MT1M, MT1F, MT1G, MT1H, MT1X, and MT1E) are involved in the mineral absorption pathway. Evidence has shown metallothioneins play critical roles in immunity, infections, protection against oxidative stress, and regulating immune cells,^35,36 and may have potential in the treatment of OA. FoxO can promote cell cycle arrest, DNA damage repair, autophagy, and scavenging of free radicals, and protect cells against various stressors such as oxidative stress.³⁷ Expression and activation of FoxO play important roles in postnatal cartilage development and protect against OA-associated cartilage damage.³⁸ In this study, the genes involved in the FoxO signaling pathway were all downregulated. Modulation of FoxO pathway in synovium may be advantageous for OA amelioration.

Hub genes such as JUN, EGR1, JUND, KLF4, MYC, and CEBPB not only played important roles in the PPI network but also acted as central nodes and transcription factors in the transcriptional regulatory network. For example, JUN and JUND, members of the activator protein-1 family, a group of transcription factors,³⁹ are involved in cell proliferation, differentiation, apoptosis, and survival via participation in the transcriptional network, and immune responses regulation.^40,41 As transcription factors, EGR1, KLF4, MYC, and CEBPB also perform diverse transcriptional functions in many cell physiological functions and cellular processes.^42
-45 Overexpressed EGR1 in chondrocytes could accelerate chondrocyte hypertrophy, prevent COl2A1 expression, and promote the release of inflammatory factors, such as MMP9 and MMP13.⁴⁶ Notably, endoplasmic reticulum stress can kinetically control CEBPB isoforms in fibroblasts. Furthermore, CEBPB pathways can be involved in antitumor immunity mediated by the myeloid-derived suppressor cells that are inhibited in triple-negative breast cancer.⁴⁷ However, the function and mechanism of these transcription factors in OA require further research.

Previous studies have demonstrated strong correlations between the infiltration of immune cells and OA.^48,49 The inflammatory membranes of OA synovium are often infiltrated by immune cells.³⁴ Th17 is an immune cell which is implicated in the pathogenesis of most common autoimmune diseases, including RA.⁵⁰ Recent findings have shown that FOSB is an essential regulator that determines the production of transforming growth factor (TGF)-β.⁵¹ TGF-β is a regulator of CD4+ T cell differentiation and determines the balance of Th17 and Treg cells for autoimmunity.⁵² KLF6 is a potential regulator of Th17 cells. Iliopoulos et al.⁵³ suggested that KLF6 is a novel protein that was deregulated in OA. Our results demonstrated that the expression of hub genes was regulated by several transcription factors that are significantly correlated with immune cell infiltration. Further research is required to provide guidance for future biological studies and clinical immunotherapy in OA.

Several genes, such as IL-1β and MMP-13, are coexpressed in cartilage, synovium, and subchondral bone tissues and play important roles in the development of OA.^5,54,55 Herein, we used 12 public data sets and experiment validation of mRNA and protein levels to demonstrate that APOD is a hub gene coexpressed and downregulated in the diseased condition as compared to normal tissue. Tew et al.⁵⁶ demonstrated that APOD expression was downregulated not only in patients with OA but also in the highly loaded and more physically damaged cartilage. Furthermore, APOD was downregulated in OA cartilage and synovium, as previously demonstrated.^57
-61 In animal experiments, APOD was also correlated with OA in synovial fluid samples and chondrocytes.^62,63 APOD is a secreted glycosylated protein and an atypical lipoprotein, which can bind to several small molecules, such as arachidonic acid, steroids, and cholesterol.⁶⁴ APOD is neuro- and cardio-protective and has important functions such as oxidation resistance, anti-inflammation, and stress resistance.⁶⁵ Furthermore, in aged mice, APOD deficiency can cause high bone turnover, low bone mass, and impair osteoblastic function.⁶⁶ Although in the small samples included in this study, the lower APOD levels correlated with OA, it remains to be confirmed in a large sample size. Also, further investigations into the potential role and underlying mechanism of APOD in OA development remain warranted.

Conclusion

The most significant module associated with OA trait was identified based on the WGCNA method for synovium tissue transcription analysis, and APOD was validated as a hub gene coexpressed and downregulated in multiple tissues. Furthermore, transcriptional regulatory-immune network was constructed, which may help to improve immune therapy for OA.

Supplemental Material

sj-csv-1-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-csv-1-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-csv-2-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-csv-2-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-cys-1-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-cys-1-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-cys1-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-cys1-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-docx-1-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-docx-1-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-docx-2-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-docx-2-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-docx-3-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-docx-3-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-pdf-1-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-pdf-1-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-xlsx-1-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-xlsx-1-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Supplemental Material

sj-xlsx-2-car-10.1177_19476035211053824 – Supplemental material for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis

Supplemental material, sj-xlsx-2-car-10.1177_19476035211053824 for Apolipoprotein D as a Potential Biomarker and Construction of a Transcriptional Regulatory-Immune Network Associated with Osteoarthritis by Weighted Gene Coexpression Network Analysis by Yong Qin, Jia Li, Yonggang Zhou, Chengliang Yin, Yi Li, Ming Chen, Yinqiao Du, Tiejian Li and Jinglong Yan in CARTILAGE

Footnotes

Author Contributions

Jinglong Yan, Yong Qin, and Yonggang Zhou conceived and designed the research study; Chengliang Yin and Yong Qin performed data analysis; Jia Li, Yi Li, Ming Chen, Yinqiao Du, and Tiejian Li collected patient samples and performed the experimental validation; Yong Qin drafted the original article; Jinglong Yan and Yonggang Zhou revised the article. All authors read and approved the final article to be published.

Acknowledgments and Funding

We thanked Professor Desi Shang for providing metagene that represented 28 immune cell types. The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was supported by the National Natural Science Foundation of China (grant number:82072472) and Medical Big Data and AI R&D Project of Chinese PLA General Hospital (2019MBD-001).

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Ethical Approval

This study was approved by the Ethics Committee of the Second Affiliated Hospital of Harbin Medical University (ethical approval number: ky2020-078). Written informed consent was obtained from each tissue donor included in this study.

ORCID iD

Jinglong Yan

Supplementary material for this article is available on the Cartilage website at .

References

Shi

Cheng

Zhang

Zhao

Shi

, et al. A small molecule promotes cartilage extracellular matrix generation and inhibits osteoarthritis development. Nat Commun. 2019;10:1914. doi:10.1038/s41467-019-09839-x.

Postler

Goronzy

Günther

Lange

Redeker

Schmitt

, et al. Which disease-related factors influence patients’ and physicians’ willingness to consider joint replacement in hip and knee OA? Results of a questionnaire survey linked to claims data. BMC Musculoskelet Disord. 2020;21:352. doi:10.1186/s12891-020-03368-1.

Martel-Pelletier

Barr

Cicuttini

Conaghan

Cooper

Goldring

, et al. Osteoarthritis. Nat Rev Dis Primers. 2016;2:16072. doi:10.1038/nrdp.2016.72.

Hügle

Geurts

. What drives osteoarthritis? Synovial versus subchondral bone pathology. Rheumatology (Oxford). 2017;56:1461-71. doi:10.1093/rheumatology/kew389.

Han

Fang

Tan

Jiang

Gong

Wang

, et al. The emerging role of fibroblast-like synoviocytes-mediated synovitis in osteoarthritis: an update. J Cell Mol Med. 2020;24:9518-32. doi:10.1111/jcmm.15669.

Mathiessen

Conaghan

. Synovitis in osteoarthritis: current understanding with therapeutic implications. Arthritis Res Ther. 2017;19:18. doi:10.1186/s13075-017-1229-9.

Lopes

Filiberti

Husain

Humphrey

. Immune contributions to osteoarthritis. Curr Osteoporos Rep. 2017;15:593-600. doi:10.1007/s11914-017-0411-y.

Klimak

Nims

Pferdehirt

Collins

Harasymowicz

Oswald

, et al. Immunoengineering the next generation of arthritis therapies. Acta Biomater. 2021;133:74-86. doi:10.1016/j.actbio.2021.03.062.

Sellam

Berenbaum

. The role of synovitis in pathophysiology and clinical symptoms of osteoarthritis. Nat Rev Rheumatol. 2010;6:625-35. doi:10.1038/nrrheum.2010.159.

10.

Min

Lee

Yoon

. Deep learning in bioinformatics. Brief Bioinform. 2017;18:851-69. doi:10.1093/bib/bbw068.

11.

Nangraj

Selvaraj

Kaliamurthi

Kaushik

Cho

Wei

. Integrated PPI- and WGCNA-retrieval of hub gene signatures shared between Barrett’s esophagus and esophageal adenocarcinoma. Front Pharmacol. 2020;11:881. doi:10.3389/fphar.2020.00881.

12.

Langfelder

Horvath

. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. doi:10.1186/1471-2105-9-559.

13.

Pei

Chen

Zhang

. WGCNA application to proteomic and metabolomic data analysis. Methods Enzymol. 2017;585:135-58. doi:10.1016/bs.mie.2016.09.016.

14.

Coutinho de Almeida

Ramos

Mahfouz

den Hollander

Lakenberg

Houtman

, et al. RNA sequencing data integration reveals an miRNA interactome of osteoarthritis cartilage. Ann Rheum Dis. 2019;78:270-7. doi:10.1136/annrheumdis-2018-213882.

15.

Ritter

Subbaiah

Bebek

Crish

Scanzello

Krastins

, et al. Proteomic analysis of synovial fluid from the osteoarthritic knee: comparison with transcriptome analyses of joint tissues. Arthritis Rheum. 2013;65:981-92. doi:10.1002/art.37823.

16.

Davis

Meltzer

. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics. 2007;23:1846-7. doi:10.1093/bioinformatics/btm254.

17.

Ritchie

Phipson

Law

Shi

, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47. doi:10.1093/nar/gkv007.

18.

Robinson

McCarthy

Smyth

. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139-40. doi:10.1093/bioinformatics/btp616.

19.

Chen

Chang

Hsu

Chen

Kuo

. Systematic analysis of transcriptomic profile of chondrocytes in osteoarthritic knee using next-generation sequencing and bioinformatics. J Clin Med. 2018;7:535. doi:10.3390/jcm7120535.

20.

Huang

Qiu

Chen

Tian

Han

, et al. Weighted gene co-expression network analysis identifies key genes from extracellular vesicles as potential prognostic biomarkers for congenital pulmonary stenosis. Mol Med Rep. 2020;22:2528-36. doi:10.3892/mmr.2020.11332.

21.

Guo

Xiao

Fang

Zhang

, et al. Identification of key pathways and genes in different types of chronic kidney disease based on WGCNA. Mol Med Rep. 2019;20:2245-57. doi:10.3892/mmr.2019.10443.

22.

Liu

Zhao

Chen

. Identification of key gene modules in human osteosarcoma by co-expression analysis weighted gene co-expression network analysis (WGCNA). J Cell Biochem. 2017;118:3953-9. doi:10.1002/jcb.26050.

23.

Langfelder

Zhang

Horvath

. Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics. 2008;24:719-20. doi:10.1093/bioinformatics/btm563.

24.

Shannon

Markiel

Ozier

Baliga

Wang

Ramage

, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498-504. doi:10.1101/gr.1239303.

25.

Han

Cho

Lee

Yun

Kim

Bae

, et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Res. 2018;46:D380-6. doi:10.1093/nar/gkx1013.

26.

Miao

Jia

Zhang

Guo

. AnimalTFDB 3.0: a comprehensive resource for annotation and prediction of animal transcription factors. Nucleic Acids Res. 2019;47:D33-8. doi:10.1093/nar/gky822.

27.

Chen

Qian

Chen

Zhang

, et al. Characterizing heterogeneity of non-small cell lung tumour microenvironment to identify signature prognostic genes. J Cell Mol Med. 2020;24:14608-18. doi:10.1111/jcmm.16092.

28.

Zhang

Zhao

Dai

Cheng

Gong

Feng

, et al. Immune landscape of colorectal cancer tumor microenvironment from different primary tumor location. Front Immunol. 2018;9:1578. doi:10.3389/fimmu.2018.01578.

29.

Zeng

Wang

Huang

Qiu

Xiao

Shi

, et al. Halofuginone inhibits TNF-α-induced the migration and proliferation of fibroblast-like synoviocytes from rheumatoid arthritis patients. Int Immunopharmacol. 2017;43:187-94. doi:10.1016/j.intimp.2016.12.016.

30.

Wei

Zhou

Wei

Zhang

Fleming

Terek

, et al. Activation of Indian hedgehog promotes chondrocyte hypertrophy and upregulation of MMP-13 in human osteoarthritic cartilage. Osteoarthr Cartil. 2012;20:755-63. doi:10.1016/j.joca.2012.03.010.

31.

Pérez-García

Carrión

Gutiérrez-Cañas

González-Álvaro

Gomariz

Juarranz

. VIP and CRF reduce ADAMTS expression and function in osteoarthritis synovial fibroblasts. J Cell Mol Med. 2016;20:678-87. doi:10.1111/jcmm.12777.

32.

Wang

Hunter

Jin

Ding

. The importance of synovial inflammation in osteoarthritis: current evidence from imaging assessments and clinical trials. Osteoarthr Cartil. 2018;26:165-74. doi:10.1016/j.joca.2017.11.015.

33.

Yuan

Xie

Wang

Shen

Feng

Liu

, et al. Screening of osteoarthritis diagnostic markers based on immune-related genes and immune infiltration. Sci Rep. 2021;11:7032. doi:10.1038/s41598-021-86319-7.

34.

Luo

Zhu

Lei

. T cells in osteoarthritis: alterations and beyond. Front Immunol. 2017;8:356. doi:10.3389/fimmu.2017.00356.

35.

Subramanian Vignesh

Deepe

Jr . Metallothioneins: emerging modulators in immunity and infection. Int J Mol Sci. 2017;18:2197. doi:10.3390/ijms18102197.

36.

Sun

Ding

Liu

Chen

, et al. Metallothionein-1 suppresses rheumatoid arthritis pathogenesis by shifting the Th17/Treg balance. Eur J Immunol. 2018;48:1550-62. doi:10.1002/eji.201747151.

37.

Almeida

Porter

. Sirtuins and FoxOs in osteoporosis and osteoarthritis. Bone. 2019;121:284-92. doi:10.1016/j.bone.2019.01.018.

38.

Akasaki

Hasegawa

Saito

Asahara

Iwamoto

Lotz

. Dysregulated FOXO transcription factors in articular cartilage in aging and osteoarthritis. Osteoarthr Cartil. 2014;22:162-70. doi:10.1016/j.joca.2013.11.004.

39.

Papoudou-Bai

Hatzimichael

Barbouti

Kanavaros

. Expression patterns of the activator protein-1 (AP-1) family members in lymphoid neoplasms. Clin Exp Med. 2017;17:291-304. doi:10.1007/s10238-016-0436-z.

40.

Papavassiliou

Musti

. The multifaceted output of c-Jun biological activity: focus at the junction of CD8 T cell activation and exhaustion. Cells. 2020;9:2470. doi:10.3390/cells9112470.

41.

Atsaves

Leventaki

Rassidakis

Claret

. AP-1 transcription factors as regulators of immune responses in cancer. Cancers (Basel). 2019;11:1037. doi:10.3390/cancers11071037.

42.

Wang

Guo

Chen

Zhao

. The role of the transcription factor EGR1 in cancer. Front Oncol. 2021;11:642547. doi:10.3389/fonc.2021.642547.

43.

Ghaleb

Yang

. Krüppel-like factor 4 (KLF4): what we currently know. Gene. 2017;611:27-37. doi:10.1016/j.gene.2017.02.025.

44.

Fallah

Brundage

Allegakoen

Shajahan-Haq

. MYC-driven pathways in breast cancer subtypes. Biomolecules. 2017;7:53. doi:10.3390/biom7030053.

45.

Yang

Qiu

. Caveolin-1 promotes pituitary adenoma cells migration and invasion by regulating the interaction between EGR1 and KLF5. Exp Cell Res. 2018;367:7-14. doi:10.1016/j.yexcr.2018.01.008.

46.

Sun

Huang

Pan

Xie

, et al. EGR1 promotes the cartilage degeneration and hypertrophy by activating the Krüppel-like factor 5 and β-catenin signaling. Biochim Biophys Acta Mol Basis Dis. 2019;1865:2490-503. doi:10.1016/j.bbadis.2019.06.010.

47.

Tanikawa

Kryczek

Xia

, et al. Aerobic glycolysis controls myeloid-derived suppressor cells and tumor immunity via a specific CEBPB isoform in triple-negative breast cancer. Cell Metab. 2018;28:87-103.e6. doi:10.1016/j.cmet.2018.04.022.

48.

Sun

Zhang

Song

Yin

Wang

, et al. IgM(+)CD27(+) B cells possessed regulatory function and represented the main source of B cell-derived IL-10 in the synovial fluid of osteoarthritis patients. Hum Immunol. 2019;80:263-9. doi:10.1016/j.humimm.2019.02.007.

49.

Zhu

Zhang

Jiang

Liu

Huang

Wei

, et al. Alterations in peripheral T cell and B cell subsets in patients with osteoarthritis. Clin Rheumatol. 2020;39:523-32. doi:10.1007/s10067-019-04768-y.

50.

Yang

Sundrud

Skepner

Yamagata

. Targeting Th17 cells in autoimmune diseases. Trends Pharmacol Sci. 2014;35:493-500. doi:10.1016/j.tips.2014.07.006.

51.

Barrett

Millena

Khan

. TGF-β effects on prostate cancer cell migration and invasion require FosB. Prostate. 2017;77:72-81. doi:10.1002/pros.23250.

52.

Lee

. The balance of Th17 versus Treg cells in autoimmunity. Int J Mol Sci. 2018;19:730. doi:10.3390/ijms19030730.

53.

Iliopoulos

Malizos

Oikonomou

Tsezou

. Integrative microRNA and proteomic approaches identify novel osteoarthritis genes and their collaborative metabolic and inflammatory networks. PLoS ONE. 2008;3:e3740. doi:10.1371/journal.pone.0003740.

54.

Chien

Tsai

Liu

Huang

Lin

Yang

, et al. Noggin inhibits IL-1β and BMP-2 expression, and attenuates cartilage degeneration and subchondral bone destruction in experimental osteoarthritis. Cells. 2020;9:927. doi:10.3390/cells9040927.

55.

Lin

Shao

Zeng

Zhao

Fang

Wang

, et al. Blocking PI3K/AKT signaling inhibits bone sclerosis in subchondral bone and attenuates post-traumatic osteoarthritis. J Cell Physiol. 2018;233:6135-47. doi:10.1002/jcp.26460.

56.

Tew

Clegg

Brew

Redmond

Hardingham

. SOX9 transduction of a human chondrocytic cell line identifies novel genes regulated in primary human chondrocytes and in osteoarthritis. Arthritis Res Ther. 2007;9:R107. doi:10.1186/ar2311.

57.

Arii

Kumon

Sugahara

Nakatani

Ikeda

Suehiro

, et al. Edaravone inhibits collagen-induced arthritis possibly through suppression of nuclear factor-kappa B. Mol Immunol. 2008;45:463-9. doi:10.1016/j.molimm.2007.05.020.

58.

Collins-Racie

Yang

Arai

Majumdar

Nagpal

, et al. Global analysis of nuclear receptor expression and dysregulation in human osteoarthritic articular cartilage: reduced LXR signaling contributes to catabolic metabolism typical of osteoarthritis. Osteoarthr Cartil. 2009;17:832-42. doi:10.1016/j.joca.2008.12.011.

59.

Tew

McDermott

Fentem

Peffers

Clegg

. Transcriptome-wide analysis of messenger RNA decay in normal and osteoarthritic human articular chondrocytes. Arthritis Rheumatol. 2014;66:3052-3061. doi:10.1002/art.38849.

60.

Neidlin

Dimitrakopoulou

Alexopoulos

. Multi-tissue network analysis for drug prioritization in knee osteoarthritis. Sci Rep. 2019;9:15176. doi:10.1038/s41598-019-51627-6.

61.

Zhang

Chen

Cui

Wang

Mao

, et al. Lnc-SAMD14-4 can regulate expression of the COL1A1 and COL1A2 in human chondrocytes. Peerj. 2019;7:e7491. doi:10.7717/peerj.7491.

62.

Garner

Kuroki

Stoker

Cook

. Expression of proteins in serum, synovial fluid, synovial membrane, and articular cartilage samples obtained from dogs with stifle joint osteoarthritis secondary to cranial cruciate ligament disease and dogs without stifle joint arthritis. Am J Vet Res. 2013;74:386-94. doi:10.2460/ajvr.74.3.386.

63.

Ashwell

O’Nan

Gonda

Mente

. Gene expression profiling of chondrocytes from a porcine impact injury model. Osteoarthr Cartil. 2008;16:936-46. doi:10.1016/j.joca.2007.12.012.

64.

Rassart

Desmarais

Najyb

Bergeron

Mounier

. Apolipoprotein D. Gene. 2020;756:144874. doi:10.1016/j.gene.2020.144874.

65.

El-Darzi

Mast

Petrov

Dao

Astafev

Saadane

, et al. Studies of ApoD(-/-) and ApoD(-/-)ApoE(-/-) mice uncover the APOD significance for retinal metabolism, function, and status of chorioretinal blood vessels. Cell Mol Life Sci. 2021;78:963-83. doi:10.1007/s00018-020-03546-3.

66.

Martineau

Najyb

Signor

Rassart

Moreau

. Apolipoprotein D deficiency is associated to high bone turnover, low bone mass and impaired osteoblastic function in aged female mice. Metabolism. 2016;65:1247-58. doi:10.1016/j.metabol.2016.05.007.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.05 MB

0.37 MB

0.08 MB

0.86 MB

0.02 MB

0.01 MB

0.10 MB

0.02 MB

0.03 MB