Abstract Aims Osteoarthritis (OA) is the most prevalent joint disease. However, the specific and definitive genetic mechanisms of OA are still unclear. Methods Tissue-related transcriptome-wide association studies (TWAS) of hip OA and knee OA were performed utilizing the genome-wide association study (GWAS) data of hip OA and knee OA (including 2,396 hospital-diagnosed hip OA patients versus 9,593 controls, and 4,462 hospital-diagnosed knee OA patients versus 17,885 controls) and gene expression reference to skeletal muscle and blood. The OA-associated genes identified by TWAS were further compared with the differentially expressed genes detected by the messenger RNA (mRNA) expression profiles of hip OA and knee OA. Functional enrichment and annotation analysis of identified genes was performed by the DAVID and FUMAGWAS tools. Results We detected 33 common genes, eight common gene ontology (GO) terms, and one common pathway for hip OA, such as calcium and integrin-binding protein 1 (CIB1) (P[TWAS] = 0.025, FC[mRNA] = -1.575 for skeletal muscle), adrenomedullin (ADM) (P[TWAS] = 0.022, FC[mRNA] = -4.644 for blood), Golgi apparatus (P[TWAS] <0.001, P[mRNA] = 0.012 for blood), and phosphatidylinositol 3' -kinase-protein kinase B (PI3K-Akt) signalling pathway (P[TWAS] = 0.033, P[mRNA] = 0.005 for blood). For knee OA, we detected 24 common genes, eight common GO terms, and two common pathways, such as histocompatibility complex, class II, DR beta 1 (HLA-DRB1) (P[TWAS] = 0.040, FC[mRNA] = 4.062 for skeletal muscle), Follistatin-like 1 (FSTL1) (P[TWAS] = 0.048, FC[mRNA] = 3.000 for blood), cytoplasm (P[TWAS] < 0.001, P[mRNA] = 0.005 for blood), and complement and coagulation cascades (P[TWAS] = 0.017, P[mRNA] = 0.001 for skeletal muscle). Conclusion We identified a group of OA-associated genes and pathways, providing novel clues for understanding the genetic mechanism of OA. Cite this article: Bone Joint Res. 2020;9(3):130–138. Keywords: Hip osteoarthritis, Knee osteoarthritis, Transcriptome-wide association study Article focus * To access the genetic mechanism of osteoarthritis (OA) using the transcriptome-wide association study (TWAS) analysis. * We performed functional enrichment and annotation analysis of the candidate genes associated with OA. Key messages * A total of 33 candidate genes were identified for hip OA, such as calcium and integrin-binding protein 1 (CIB1) and adrenomedullin (ADM). * For knee OA, we detected 24 candidate genes, such as histocompatibility complex, class II, DR beta 1 (HLA-DRB1), and Follistatin-like 1 (FSTL1). Strengths and limitations of this study * TWAS has a boosting power to detect novel disease genes. * One limitation is that TWAS may be too low power to detect the causal loci without cis-expression effects on target disease. Introduction Osteoarthritis (OA) is mainly characterized by the degeneration of articular cartilage in the knee and hip,^[44]1 and knee and hip OA occur in approximately 6% and 3% of Americans aged 30 years or older, respectively.^[45]2 The clinical symptoms of OA include joint pain, swelling, stiffness, and limited motion. The burden of OA is not only on healthcare and society costs of the patients and families, but also includes the patients’ quality of life, mental health, and emotional relationship.^[46]3,[47]4 A review study has shown that the heritability of genetic factors in the development of OA was estimated at approximately 50% in twin studies and familial studies.^[48]5 Through using data across 16.5 million variants from the UK Biobank, a group of OA-associated genes have been identified by genome-wide association studies (GWAS), such as MAP2K6 and ZNF345.^[49]6 However, GWAS usually focus on a few genome-wide significant loci with large genetic effects, while they are likely to miss many biological true positive loci with moderate or weak genetic effects. In addition, a previous study found that a large section of genetic variants identified by GWAS was enriched in non-coding regulatory genomic regions, which leads to difficulty in interpreting the genetic effects.^[50]7 It has been demonstrated that the messenger RNA (mRNA) expression levels were under genetic control.^[51]8 Expression quantitative trait loci (eQTLs) are the genetic variants, which could explain the variations in expression levels of mRNAs.^[52]9 eQTLs provide a novel way to uncover the biological mechanism of identified genetic variants underlying diseases.^[53]10 Integrating the GWAS and eQTL data could boost the power to discover novel disease-associated genes. In the previous research, the transcriptome-wide association study (TWAS) was developed to explore gene-trait relationships, using publicly available GWAS results and eQTL reference datasets.^[54]11 The expression levels of thousands of genes were predicted by TWAS, and subsequently were used to evaluate the associations between gene expression levels and target diseases. Different from GWAS testing millions of single nucleotide polymorphisms (SNPs), TWAS can greatly reduce the burden of multiple comparisons in statistical analysis and enhance the ability of GWAS to detect novel disease genes.^[55]12 As a supplement to traditional GWAS analyses, five novel susceptibility loci associated with cutaneous squamous cell carcinoma (cSCC) were validated by TWAS analysis.^[56]13 To identify novel OA-associated genes, we performed a large-scale integrative analysis of TWAS and mRNA expression profiles for hip OA and knee OA. Using previous large-scale GWAS data and eQTL reference data of skeletal muscle and blood, TWAS was performed to detect novel candidate genes, the predicted expression levels of which were associated with OA. The genes identified by TWAS were further subjected to gene ontology (GO) and pathway enrichment analysis. To further confirm the functional relevance of identified genes, the TWAS results were compared with the mRNA expression profiles of OA to detect common genes, GO terms, and pathways shared by TWAS, and mRNA expression profiles of OA. Methods GWAS summary data of hip and knee OA The GWAS summary data of hip and knee OA were derived from the published studies.^[57]6 Briefly, 2,396 hospital-diagnosed hip OA patients and 9,593 controls, and 4,462 hospital-diagnosed knee OA patients and 17,885 controls were all derived from the UK Biobank. The diagnose standard of hospital-diagnosed OA coding in the UK Biobank comes from the International Classification of Diseases (ICD)-10 code captured from Hospital Episode Statistics (HES) data. Finally, 16,122,076 SNPs for hip OA and 16,309,199 SNPs for knee OA were applied for GWAS analysis. Association testing was conducted by SNPTEST v2.5.2 (University of Oxford, Oxford, UK).^[58]14 A detailed description of the sample characteristics, experimental design, and statistical analysis can be found in the published study.^[59]6 mRNA expression profiles of hip cartilage with OA The mRNA expression profiles data of articular cartilage from hip OA and traumatic femoral neck fracture patients were obtained from a previous study.^[60]15 Briefly, articular cartilage specimens of hip were collected from nine OA patients and ten traumatic femoral neck fracture patients undergoing joint replacement surgery. OA status was confirmed using clinical examination and joint score.^[61]16 The subjects with hip scoring ≤ 1 were considered healthy samples, while those scoring ≥ 5 were classified as case samples. Differentially expressed genes were identified significantly with a fold change (FC) ≥ 1.5 and p < 0.05. The detailed sample characteristics and experimental design can be found in the previous study.^[62]15 mRNA expression profile of knee cartilage with OA The mRNA expression profile data of knee OA were obtained from a previously published study.^[63]17 Briefly, normal human knee cartilage tissues of 18 people without history of joint disease or trauma were procured from tissue banks. OA-affected cartilage specimens were harvested from 20 donors accepting knee arthroplasty surgery. Differentially expressed genes were identified significantly with an adjusted p-value of < 0.05 and log[2]FC ≥ 1. Detailed descriptions of sample characteristics, experiment design, and statistical analysis of this dataset are available in the study by Fisch et al.^[64]17 TWAS of hip OA and knee OA TWAS of hip OA and knee OA were performed by the FUSION software (Harvard T.H. Chan School of Public Health, Boston, Massachusetts, USA) through integrating the OA GWAS summary data and precomputed gene expression weights of skeletal muscle and blood (including peripheral blood and whole blood).^[65]18 A previous study performed TWAS analysis and identified the causal genes to pancreatic cancer (PC).^[66]19 The gene expression weights reference of skeletal muscle and blood were obtained from the FUSION website.^[67]18,[68]20 Peripheral blood and skeletal muscle were also used in previous biological studies of OA.^[69]21,[70]22 The gene expression weights of skeletal muscle were derived from the Genotype-Tissue Expression (GTEx) Project (version 7; National Disease Research Interchange, Philadelphia, Pennsylvania, USA; n = 361).^[71]23 The gene expression weights of peripheral blood and whole blood were collected from 1,245 unrelated control individuals from the Netherlands Twin Registry study and 1,264 subjects from the Young Finns Study.^[72]24–[73]27 For eQTL data, a total of 3,637,328 SNPs (corresponding to 7,408 genes), 1,120,437 SNPs (corresponding to 2,454 genes), and 2,044,474 SNPs (corresponding to 4,700 genes) were used for the TWAS of OA in skeletal muscle, peripheral blood, and whole blood tissues, respectively. Firstly, the gene expression weights were calculated using the prediction models of FUSION. Then, the calculated expression weights were combined with GWAS results to impute association statistics between gene expression levels and target diseases. The SNP-expression weights in the 1 Mb cis loci of the gene for a given gene were computed using Bayesian sparse linear mixed model (BSLMM).^[74]28 The association testing statistics between predicted gene expression and target diseases was calculated as Z[TWAS] = w'Z/(w'Lw)^1/2. ‘Z’ denotes the scores of OA and ‘w’ denotes the weights. ‘L’ denotes the SNP-correlation linkage disequilibrium (LD) matrix. In this analysis, we accounted for LD among SNPs and viewed the imputed gene expression data as a linear model of genotypes with weights. A TWAS p-value was calculated for each gene within skeletal muscle and blood. The genes with p < 0.05 were considered as significant. Then, the hip and knee OA-associated genes identified by TWAS were further compared with the differentially expressed genes detected by the mRNA expression profiles of hip and knee OA, respectively. ‘FC[mRNA]’ and ‘P[mRNA]’ represent the FCs and p-value in mRNA expression profiles. A flowchart illustrating the study design is shown in [75]Figure 1. Fig. 1. [76]Fig. 1 [77]Open in a new tab Detailed flowchart of the study design. eQTL, expression quantitative trait locus; GO, gene ontology; GWAS, genome-wide association study; mRNA, messenger RNA; OA, osteoarthritis; TWAS, transcriptome-wide association study. Functional enrichment and annotation analysis Gene ontology and pathway enrichment analysis of the genes identified by TWAS and mRNA expression profiles were performed by the DAVID tool.^[78]29–[79]31 We compared the analysis results of TWAS and mRNA expression profiles to screen out the common genes, GO terms, and pathways, which were shared by TWAS and mRNA expression profiles. Functional Mapping and Annotation of Genome-wide Association Studies (FUMAGWAS)^[80]32 was used to annotate, prioritize, visualize, and interpret the function of the common genes shared by TWAS and mRNA expression profiles.^[81]33 Results TWAS results of hip and knee OA For hip OA, TWAS identified 182 genes for skeletal muscle and 390 genes for blood with p < 0.05, such as ArfGAP with SH3 domain, ankyrin repeat and PH domain 3 (ASAP3) (P[TWAS] < 0.001 for skeletal muscle), high-density lipoprotein binding protein (HDLBP) (P[TWAS] < 0.001 for skeletal muscle), transcription elongation factor A3 (TCEA3) (P[TWAS] < 0.001 for blood), and serine/threonine kinase 25 (STK25) (P[TWAS] < 0.001 for blood) (Supplementary Table i). For knee OA, TWAS identified 180 genes for skeletal muscle and 410 genes for blood with p < 0.05, such as RP11-347C12.1 (P[TWAS] < 0.001 for skeletal muscle), centrosomal protein 250 (CEP250) (P[TWAS] < 0.001 for skeletal muscle), RWD domain containing 2B (RWDD2B) (P[TWAS] < 0.001 for blood), and ubiquinol-cytochrome c reductase complex assembly factor 1 (UQCC) (P[TWAS] < 0.001 for blood) (Supplementary Table ii). The top ten significant genes of hip and knee OA identified by TWAS are shown in [82]Table I. Table I. List of the top ten significant genes identified by transcriptome-wide association studies for hip and knee osteoarthritis (p < 0.05). OA Tissue Gene CHR Z[TWAS] P[TWAS][83]^* Hip Skeletal muscle ASAP3 1 −4.9331 < 0.001 Hip Blood TCEA3 1 −4.5970 < 0.001 Hip Skeletal muscle HDLBP 2 −4.3684 < 0.001 Hip Blood STK25 2 −4.0764 < 0.001 Hip Skeletal muscle ITIH4-AS1 3 3.8985 < 0.001 Hip Blood ST3GAL3 1 −3.8807 < 0.001 Hip Blood PNO1 2 −3.8675 < 0.001 Hip Skeletal muscle GRINA 8 3.8140 < 0.001 Hip Skeletal muscle RP5-966M1.6 3 3.7270 < 0.001 Hip Blood PMM2 16 −3.7253 < 0.001 Knee Blood RWDD2B 21 −4.1737 < 0.001 Knee Blood UQCC 20 3.9040 < 0.001 Knee Blood N6AMT1 21 −3.8769 < 0.001 Knee Blood MAPK3 16 −3.8468 < 0.001 Knee Skeletal muscle RP11-347C12.1 16 3.8161 < 0.001 Knee Blood SSH1 12 −3.7411 < 0.001 Knee Blood GPX4 19 3.6084 < 0.001 Knee Blood PLIN2 9 −3.5411 < 0.001 Knee Blood CDK7 5 3.4976 < 0.001 Knee Blood MTHFSD 16 −3.4959 < 0.001 [84]Open in a new tab ^* Each PTWAS value was calculated by transcriptome-wide association study analysis. CHR, chromosome; OA, osteoarthritis; TWAS, transcriptome-wide association study. Gene ontology and pathway enrichment analysis For hip OA, GO and pathway enrichment analysis results of the genes identified by TWAS are shown in Supplementary Table iii. DAVID detected 11 GO terms for skeletal muscle and 31 GO terms for blood with p < 0.05, such as UFM1 hydrolase activity (P[TWAS] = 0.016 for skeletal muscle), DNA damage checkpoint (P[TWAS] = 0.019 for skeletal muscle), and membrane (P[TWAS] < 0.001 for blood). Pathway enrichment analysis detected four pathways for blood, such as bile secretion (P[TWAS] = 0.010) and glycosaminoglycan biosynthesis-chondroitin sulfate/dermatan sulfate (P[TWAS] = 0.006). For knee OA, ten GO terms for skeletal muscle and 58 GO terms for blood were detected with p < 0.05, such as protein binding (P[TWAS] < 0.001), poly(A) RNA binding (P[TWAS] = 0.013), and cytosol (P[TWAS] < 0.001) (Supplementary Table iv). Pathway enrichment analysis of the significant genes identified one pathway for skeletal muscle and 14 pathways for blood (P[TWAS] < 0.05), such as complement and coagulation cascades (P[TWAS] = 0.017), influenza A (P[TWAS] = 0.006), and viral carcinogenesis (P[TWAS] = 0.009) (Supplementary Table iv). Comparative analysis of TWAS and mRNA expression profiles We further compared the analysis results of TWAS and mRNA expression profiles. For hip OA, we detected 33 common genes shared by the TWAS and mRNA expression profiles, such as calcium and integrin-binding protein 1 (CIB1) (P[TWAS] = 0.025, FC[mRNA] = -1.575 for skeletal muscle), adrenomedullin (ADM) (P[TWAS] = 0.022, FC[mRNA] = -4.644 for blood), and forkhead box C1 (FOXC1) (P[TWAS] = 0.029, FC[mRNA] = 1.527 for blood) ([85]Table II). In addition, we detected eight common GO terms and one common pathway, such as cell-cell adherens junction (P[TWAS] = 0.037, P[mRNA] = 0.016 for skeletal muscle), Golgi apparatus (P[TWAS] < 0.001, P[mRNA] = 0.012 for blood), and PI3K-Akt signalling pathway (P[TWAS] = 0.033, P[mRNA] = 0.005 for blood) ([86]Table III). The heat map of those common genes of hip OA is shown in [87]Figure 2. Table II. Common genes between the significant genes identified by transcriptome-wide association studies and the differentially expressed genes identified by messenger RNA expression profiles for hip osteoarthritis. Tissue Gene P[TWAS][88]^* P[mRNA]^[89]† FC[mRNA] Skeletal muscle CIB1 0.025 0.006 −1.575 CYP4V2 0.006 0.006 1.782 HRCT1 0.028 0.007 1.636 SLC14A1 0.027 0.009 2.633 STEAP2 0.050 0.006 −1.838 TFPI 0.033 0.006 −2.637 Blood ADM 0.022 0.006 −4.644 ASPH 0.007 0.006 −1.569 EPB41L2 0.038 0.006 2.295 GLT25D2 0.024 0.006 2.816 GNG2 0.003 0.008 1.564 IFITM3 0.001 0.006 −2.309 IL8 0.037 0.006 −2.774 LSM14A 0.031 0.006 1.527 RAB3IP 0.044 0.006 −1.952 SLC14A1 0.015 0.009 2.633 ZHX2 0.041 0.009 −1.603 FOXC1 0.029 0.007 1.527 GALNT12 0.025 0.006 −1.548 GAS1 0.025 0.006 4.096 ID3 0.001 0.006 2.267 ITGB7 0.042 0.006 1.684 KLHL36 0.046 0.006 −1.598 NT5DC1 0.038 0.006 1.609 PAM 0.032 0.006 1.557 RAB21 0.009 0.006 −1.650 RNMT 0.042 0.006 −1.836 RXRA 0.011 0.006 1.889 TCEA3 < 0.001 0.006 1.843 UBE2H 0.007 0.006 −1.693 USP10 0.004 0.006 −1.699 FBL 0.003 0.009 −1.635 RHOBTB3 0.038 0.006 −1.696 [90]Open in a new tab ^* Each P[TWAS] value was calculated by transcriptome-wide association study analysis. ^† Each P[mRNA] value was derived from the published studies. FC, fold change; mRNA, messenger RNA; TWAS, transcriptome-wide association study. Table III. Common gene ontology terms and pathways between the significant genes identified by transcriptome-wide association studies and the differentially expressed genes identified by messenger RNA expression profiles for hip osteoarthritis. Category Tissue Name P[TWAS][91]^* P[mRNA]^[92]† GO Skeletal muscle GO:0005913~cell-cell adherens junction 0.037 0.016 Blood GO:0000139~Golgi membrane < 0.001 0.012 GO:0005794~Golgi apparatus < 0.001 < 0.001 GO:0005515~protein binding 0.008 < 0.001 GO:0005925~focal adhesion 0.018 < 0.001 GO:0070062~extracellular exosome 0.019 < 0.001 GO:0005789~endoplasmic reticulum membrane 0.024 0.002 GO:0005654~nucleoplasm 0.031 0.033 Pathway Blood hsa04151:PI3K-Akt signalling pathway 0.033 0.005 [93]Open in a new tab ^* Each PTWAS value was calculated by gene ontology and pathway analysis using DAVID for the significant genes identified by transcriptome-wide association study analysis. ^† Each PmRNA value was calculated by gene ontology and pathway analysis using DAVID for the significant genes identified by messenger RNA expression profiles. GO, gene ontology; mRNA, messenger RNA; PI3K-Akt, phosphatidylinositol 3’ -kinase-protein kinase B; TWAS, ranscriptome-wide association study. Fig. 2. [94]Fig. 2 [95]Open in a new tab Gene expression heat map of the identified common genes shared by transcriptome-wide association studies (TWAS) and messenger RNA (mRNA) expression data of hip osteoarthritis (OA). For knee OA, 24 common genes were identified, such as major histocompatibility complex, class II, DR beta 1 (HLA-DRB1) (P[TWAS] = 0.040, FC[mRNA] = 4.062 for skeletal muscle), general transcription factor IIE subunit 1 (GTF2E1) (P[TWAS] = 0.043, FC[mRNA] = 2.368 for skeletal muscle), Follistatin-like 1 (FSTL1) (P[TWAS] = 0.048, FC[mRNA] = 3.000 for blood), and beta-1,3-galactosyltransferase 6 (B3GALT6) (P[TWAS] < 0.001, FC[mRNA] = 2.221 for blood) ([96]Table IV). In addition, we detected eight common GO terms and two common pathways, such as protein binding (P[TWAS] < 0.001, P[mRNA] < 0.001 for skeletal muscle), cytoplasm (P[TWAS] < 0.001, P[mRNA] = 0.005 for blood), complement and coagulation cascades (P[TWAS] = 0.017, P[mRNA] = 0.001 for skeletal muscle), and viral myocarditis (P[TWAS] = 0.009, P[mRNA] = 0.050 for blood) ([97]Table V). A heat map showing the expression of common genes of knee OA is shown in [98]Figure 3. Table IV. Common genes identified by comparing the gene ontology and pathway enrichment analysis results of transcriptome-wide association studies and messenger RNA expression profiles for knee osteoarthritis. Tissue Gene P[TWAS][99]^* P[mRNA]^[100]† FC[mRNA] Skeletal muscle GTF2E1 0.043 < 0.001 2.368 HLA-DRB1 0.040 0.045 4.062 RRAGD 0.001 < 0.001 0.403 PDE3A 0.048 0.003 2.826 SERPINA5 0.002 < 0.001 3.078 RGS19 0.002 0.011 3.001 Blood B3GALT6 0.010 < 0.001 2.221 BHLHE40 0.045 < 0.001 0.196 KCNAB1 0.021 < 0.001 0.496 NDUFB2 0.029 0.043 2.270 SPON1 0.007 0.025 2.074 BFSP1 0.005 0.003 2.605 CAMK2N1 0.018 0.035 2.012 CAPZB 0.007 0.003 2.201 DOCK10 0.006 0.002 4.971 FSTL1 0.048 < 0.001 3.000 GM2A 0.005 < 0.001 2.019 HLA-DRB1 0.014 0.048 4.062 ZKSCAN4 0.004 < 0.001 2.318 ZSCAN16 0.008 0.001 2.324 PLIN2 < 0.001 < 0.001 0.499 AFAP1L2 0.019 0.006 0.486 TMEM107 0.037 < 0.001 0.402 RAB31 0.042 < 0.001 3.123 [101]Open in a new tab ^* Each PTWAS value was calculated by transcriptome-wide association study analysis. ^† Each PmRNA value was derived from the published studies. FC, fold change; mRNA, messenger RNA; TWAS, transcriptome-wide association study. Table V. Common gene ontology terms and pathways identified by comparing the gene ontology and pathway enrichment analysis results of transcriptome-wide association studies and messenger RNA expression profiles for knee osteoarthritis. Category Tissue Name P[TWAS][102]^* P[mRNA]^[103]† GO Skeletal muscle GO:0005515~protein binding < 0.001 < 0.001 Blood GO:0005515~protein binding < 0.001 < 0.001 GO:0005737~cytoplasm < 0.001 0.005 GO:0070062~extracellular exosome 0.002 < 0.001 GO:0015629~actin cytoskeleton 0.014 0.026 GO:0006468~protein phosphorylation 0.017 0.033 GO:0060333~interferon-gamma-mediated signalling pathway 0.017 0.047 GO:0046982~protein heterodimerization activity 0.033 0.025 Pathway Skeletal muscle hsa04610:Complement and coagulation cascades 0.017 0.001 Blood hsa05416:Viral myocarditis 0.009 0.050 [104]Open in a new tab ^* Each PTWAS value was calculated by gene ontology and pathway analysis using DAVID for the significant genes identified by transcriptome-wide association study analysis. ^† Each PmRNA value was calculated by gene ontology and pathway analysis using DAVID for the significant genes identified by messenger RNA expression profiles. GO, gene ontology; mRNA, messenger RNA; TWAS, transcriptome-wide association study. Fig. 3. [105]Fig. 3 [106]Open in a new tab Gene expression heat map of the identified common genes shared by transcriptome-wide association studies (TWAS) and messenger RNA (mRNA) expression data of knee osteoarthritis (OA). Discussion To improve the understanding of the genetic aetiology of OA, we performed a TWAS of hip OA and knee OA through integrating the GWAS summary data and precomputed gene expression weights of skeletal muscle and blood. The TWAS results were further compared with mRNA expression profiles of OA cartilage, which detected multiple common genes, GO terms, and pathways shared by TWAS and mRNA expression profiles. One of the important candidate genes identified in this study is ADM. This gene encodes preprohormone, which is cleaved into two biologically active peptides, including ADM and proadrenomedullin N-terminal 20 peptide. Previous studies have shown that ADM is expressed in bone and joint structures including cartilage and synovium.^[107]34,[108]35 ADM is one of the top-ranking differentially expressed genes in OA bone, and plays a role in osteoblast and osteocyte differentiation and function.^[109]36 Downregulation of ADM is capable of inhibiting adipogenesis and osteoblastogenesis.^[110]37 FOXC1 belongs to the forkhead family of transcription factors (TFs), which is characterized by a distinct DNA-binding forkhead domain. A systematic analysis of six Gene Expression Omnibus (GEO) databases for synovial expression profiling identified the top ten TFs covering the most downstream differentially expressed genes and crucial TFs involved in the development of OA, including FOXC1.^[111]38 MicroRNAs have been identified in the development of OA and microRNA (miR)-138 has been reported to be involved with osteogenesis and regulation of chondrocyte phenotype.^[112]39,[113]40 Research has also shown that miR-138-5p promotes cartilage degradation induced by interleukin (IL)-1β in human chondrocytes, possibly by targeting FOXC1.^[114]41 FSTL1/follistatin-related protein (FRP), an extracellular protein, was found in mesenchymal stem cells (MSCs) and cartilage. The FSTL1 mRNA and protein levels in the serum and synovial fluid were significantly higher in OA patients than in controls.^[115]42 Therefore, FSTL1 gene expression may be increased in OA patients.^[116]42 In addition, the findings show that FSTL1 is an important proinflammatory factor in the pathogenesis of OA through activating the canonical NF-κB pathway and enhancing synoviocyte proliferation, which may lead to the development of novel strategies for cartilage repair and be a promising target for the treatment of OA.^[117]43,[118]44 HLA-DRB1 belongs to the major histocompatibility complex (HLA) class II beta chain paralogs. In the previous study, two SNPs (rs7775228 and rs10947262) in a region containing HLA class II/III genes associated with susceptibility to knee OA were identified through a GWAS and a replication, using a total of approximately 4,800 Japanese subjects.^[119]45 However, the rs7775228 and rs10947262 variants were not associated with risk of knee OA in European populations compared with Japanese individuals.^[120]46 The previous study showed that DR2 and DR5 were associated with OA, which hinted at LD between HLA-DRB1 genes and genes involved in the pathogenesis of OA.^[121]47 We also detected eight common GOs associated with hip and knee OA, respectively. For instance, protein binding (GO:0005515) was significant in blood of hip OA (P[TWAS] = 0.008, P[mRNA] < 0.001), skeletal muscle of knee OA (P[TWAS] < 0.001, P[mRNA] < 0.001), and blood of knee OA (P[TWAS] < 0.001, P[mRNA] < 0.001). To our knowledge, cold-inducible RNA-binding protein (CIRP) is a kind of inflammatory cytokine. In addition, research has found that the concentration of synovial fluid CIRP is closely associated with the synovial inflammation of OA and CIRP could be used as a potential marker for synovial inflammation of OA.^[122]48 The binding protein of the growth arrest and DNA damage-inducible protein 45 β (GADD45β) gene is down-regulated in ageing articular cartilage and chondrocyte clusters in OA cartilage.^[123]49 In all, these results further confirm the role of Golgi modifications and apoptosis in OA pathogenesis. Although the heritability of OA is partly explained by GWAS, they still could ignore the genetic variants with expression-trait associations. In this study, TWAS analysis was used to identify novel genes associated with OA in the mRNA expression levels and gene expression profiles were used to verify the results. Compared with previous studies,^[124]50,[125]51 TWAS is prone to spurious prioritization based on the expression data from OA-related tissues. In addition, TWAS was more accurate in prioritizing candidate causal genes than simple baselines. TWAS has been successful in identifying many genes and has the required boosting power to detect novel disease genes.^[126]12 There are three limitations that need to be noted in this study. Firstly, TWAS was developed to identify the genes, the regulated expression of which is associated with target diseases. TWAS may be too low power to detect the causal loci without cis-expression effects on target disease. Secondly, the objective of this study was to scan novel candidate genes related with OA. Further functional studies are needed to confirm our findings and clarify the potential biological mechanisms underlying OA, as detailed in a previous experimental study by Zhang et al.^[127]52 Thirdly, fracture may have had an impact on the hip cartilage expression profiles. However, in previous studies the individuals undergoing hip arthroplasty following femoral neck fracture were selected as the normal controls to perform the expression profiles.^[128]53,[129]54 Therefore, the hip cartilage expression profile datasets, in which the patients with femoral neck fracture were used as the control cartilage specimens,^[130]15 could be used in this study analysis. Given these limitations, our results should be interpreted with caution and further studies are needed to confirm our findings. In conclusion, this study combined TWAS and gene expression profiling datasets to identify the candidate gene associated with OA. We have identified 33 and 24 common genes in hip OA and knee OA, respectively, such as FOXC1, ADM, FSTL1, and HLA-DRB1. In view of these limitations, the results should be explained with caution. Therefore, we need further studies to verify our findings and reveal the potential effect of identified genes in the development of OA. Footnotes Author Contributions: X. Qi: Wrote the manuscript, Designed the study, Analyzed the data. F. Yu: Wrote the manuscript. Y. Wen: Wrote the manuscript. P. Li: Wrote the manuscript. B. Cheng: Wrote the manuscript. M. Ma: Wrote the manuscript. S. Cheng: Wrote the manuscript. L. Zhang: Wrote the manuscript. C. Liang: Wrote the manuscript. L. Liu: Wrote the manuscript. F. Zhang: Wrote the manuscript, Designed the study. Ethical review statement: This study did not require ethical approval. Follow us [131]@BoneJointRes Supplementary Material Supplementary tables showing: transcriptome-wide association study results of hip osteoarthritis; transcriptome-wide association study results of knee osteoarthritis; gene ontology and pathway enrichment analysis results of the significant genes identified by transcriptome-wide association studies for hip osteoarthritis; gene ontology and pathway enrichment analysis results of the significant genes identified by transcriptome-wide association studies for knee osteoarthritis. Funding statement The authors report institutional grants related to this study (paid to Xi'an Jiaotong University) from the National Natural Science Foundation of China, the Key projects of international cooperation among governments in scientific and technological innovation, the Fundamental Research Funds for the Central Universities, and the Natural Science Basic Research Plan in Shaanxi Province of China. No benefits in any form have been received or will be received from a commercial party related directly or indirectly to the subject of this article. References