Understanding the complex genetic architecture connecting rheumatoid arthritis, osteoporosis and inflammation: discovering causal pathways

Abstract Rheumatoid arthritis (RA) and osteoporosis (OP) are two comorbid complex inflammatory conditions with evidence of shared genetic background and causal relationships. We aimed to clarify the genetic architecture underlying RA and various OP phenotypes while additionally considering an inflammatory component, C-reactive protein (CRP). Genome-wide association study summary statistics were acquired from the GEnetic Factors for OSteoporosis Consortium, Cohorts for Heart and Aging Research Consortium and UK Biobank. Mendelian randomization (MR) was used to detect the presence of causal relationships. Colocalization analysis was performed to determine shared genetic variants between CRP and OP phenotypes. Analysis of pleiotropy between traits owing to shared causal single nucleotide polymorphisms (SNPs) was performed using PL eiotropic A nalysis under CO mposite null hypothesis (PLACO). MR analysis was suggestive of horizontal pleiotropy between RA and OP traits. RA was a significant causal risk factor for CRP (β = 0.027, 95% confidence interval = 0.016–0.038). There was no evidence of CRP→OP causal relationship, but horizontal pleiotropy was apparent. Colocalization established shared genomic regions between CRP and OP, including GCKR and SERPINA1 genes. Pleiotropy arising from shared causal SNPs revealed through the colocalization analysis was all confirmed by PLACO. These genes were found to be involved in the same molecular function ‘protein binding’ (GO:0005515) associated with RA, OP and CRP. We identified three major components explaining the epidemiological relationship among RA, OP and inflammation: (1) Pleiotropy explains a portion of the shared genetic relationship between RA and OP, albeit polygenically; (2) RA contributes to CRP elevation and (3) CRP, which is influenced by RA, demonstrated pleiotropy with OP.


Introduction
Rheumatoid arthritis (RA) is a debilitating autoimmune disease with a prevalence of approximately 1% in the adult population (1), where females are at least twice as likely to be affected than males. The clinical appearance of RA typically manifests symmetrically in the smaller synovial joints of the hand and feet and may also affect large joints (2). Clinical implications of RA also extend to a variety of tissues, including lungs and eyes, and lead to bone and cartilage deterioration including osteoporosis (OP) resulting in a high burden of disease (2). There is a considerable disability and higher mortality (1,3). The complexity of RA arises from its multifactorial nature (4) and it has a heritability of approximately 60% (5).
OP is similarly multifactorial and demonstrates shared risk factors to RA (6). Dual X-ray absorptiometry is the standard technique used to measure bone mineral density (BMD) revealing OP and therefore fracture risk (7), though not the sole indicator of OP. In addition, bone architecture or quality is an important component of OP and fracture prediction (8). The heritability of BMD measures and fracture risk varies by skeletal site measured: 50-85 and 25-58%, respectively (9).
Consistent with previous studies including our own, the epidemiological association between RA and OP phenotypes is well established and has been observed across different samples considering various skeletal sites (10)(11)(12). The association between RA and OP is partially attributed to inflammatory factors altering bone resorption and leading to bone loss resulting from systemic inflammation seen in RA (13,14). Creactive protein (CRP) is an inflammatory blood marker commonly used as a biomarker of RA; higher levels of CRP have been reported in subjects having both RA and OP in comparison to subjects with RA only (15). Yet, regardless of RA and independent of other covariates, an association between CRP and OP has been reported in a large sample (16). Moreover, a longitudinal study suggests that increased CRP levels predict faster BMD decline in women (17). Multiple logistic regression on the ELSA dataset (Wave 6) where RA served as the dependent variable and OP served as an independent covariate in addition to age, BMI and hsCRP. The study sample is composed of 10 183 individuals in total of which 629 were diagnosed with RA and 786 were diagnosed as OP cases; hsCRP levels were measured in 6125 individuals. OR was estimated per SD of the phenotype variation. The independent variables initially included in the analysis were age, sex, BMI, hsCRP, fibrinogen levels, WBC count and OP. Only statistically significant results are reported in the table.
The three conditions-RA, OP and CRP-appear closely related. Thus, the extent to which genetic factors contribute to their correlation and the identification of shared genetic factors may shed light on pathogenic pathways of interest. Preliminary genetic analyses reported pleiotropy between RA and BMD (18) and between CRP and BMD (19). Still, our knowledge is virtually non-existent in regard to the genetic architecture behind the RA/OP/CRP trilateral relationship.
The main aim of this study was to explore the nature and possible causal genetic relationships between these two pathological musculoskeletal conditions and an inf lammatory biomarker. Grounding on extensive clinical data, we explored the following hypotheses: RA, an inf lammatory condition, necessarily contributes to raised CRP and this likely contributes to OP through bone turnover mechanisms and, therefore, the following sequence of pathological events is expected: RA → raised CRP; raised CRP→OP; RA→OP.

Epidemiological assessment
The English Longitudinal Study of Ageing (ELSA) dataset comprised 4556 males and 5627 females of European ancestry with a mean age of 66 ± 10 and age range of 28-89 years. Binary multiple logistic regression analysis showed that RA was significantly and independently associated with OP [odds ratio (OR): 1.6627, 95% confidence interval (95% CI) = 1.1417-2.3604, P = 5.94E−03] and high-sensitivity CRP (hsCRP) (OR: 1.211, 95% CI = 1.0271-1.2082, P = 4.95E−03) regardless of age and body mass index (BMI) ( Table 1). Sex, fibrinogen circulating levels and white blood cell count were also included at the entry stage as independent covariates of RA but were not retained in the final equation.

Mendelian randomization analysis
To detect causality, RA served as the exposure variable, using single nucleotide polymorphisms (SNPs) from the meta-analysis as instruments, and a selected OP trait or CRP served as the outcome variable. Next, CRP served as the exposure variable and the selected OP trait as the outcome variable.
Results were considered valid if the heterogeneity assumption was not violated and I 2 Gx > 90% was achieved, thereby suggesting instrument variable (IV) strength and preserving Mendelian randomization (MR) assumptions. In addition to assessing causality and providing the causal estimate, the MR Egger (MRE) test estimates an intercept, which indicates horizontal pleiotropy if significant (P < 0.05).

RA→OP phenotypes
Implementing the inverse-variance-weighted (IVW) approach, RA appeared to causally predict OP and lower BMD measures (P < 0.001 for all; Supplementary Material, Table S1). Upon further investigation, using the MRE approach (Table 2), we clarified that RA appeared to be a causal risk factor for hip BMD (β = −0.036, 95% CI = −0.063 to −0.008, P-value = 1.10E−02), with no evidence of horizontal pleiotropy (P-value = 0.295). However, after Bonferroni correction for multiple comparisons, the association between RA and hip BMD did not achieve a significant MR estimate (α ≤ 0.00625) and, therefore, cannot be considered to be causal at this stage.
The other OP phenotypes did not show causal associations with RA, but each association revealed evidence of horizontal pleiotropy (Table 2). Specifically, in MR analysis of RA to OP fracture, the MRE intercept was highly statistically significant (P-value < 0.001), suggesting the presence of horizontal pleiotropy. Similarly, horizontal pleiotropy was apparent in the OP phenotype (P = 0.001) and each BMD assessment [total BMD P = 0.002, spine BMD P = 0.001, arm BMD P = 0.006 and UK Biobank (UKB) heel BMD P < 0.001].

CRP→OP phenotypes
In assessing IVW, elevated CRP appeared to be causally predictive of each OP trait tested (P ≤ 0.001 for all; Supplementary Material, Table S2). However, using the MRE approach, elevated CRP was not causally predictive of any OP phenotype tested; instead, horizontal pleiotropy appeared to be significantly involved in the association between elevated CRP and either the OP phenotype or OP fracture (P = 0.002 for both; Table 3). Similarly, BMD phenotypes displayed significant pleiotropic effects (total BMD P = 0.002, spine BMD P = 0.021, hip BMD P < 0.001, arm BMD P = 0.006 and UKB heel BMD P < 0.001; Table 3).

Colocalization and gene-set enrichment analyses
While a causal relationship between CRP and OP, CRP->OP phenotypes was not apparent, horizontal pleiotropy was consistently seen. In an attempt to identify potential common genes involved, we performed colocalization analyses, which provided significant strong (PP.H4 ≥ 75%) evidence of shared SNPs in five regions on chromosomes 1, 2, 9, 14 and 22 (Table 4). Moderate evidence of shared SNPs can be seen in Supplementary Material, Table S3, and low evidence of shared SNPs is observed in Supplementary Material, Table S4. The highest posterior probability (PP) of shared common SNPs (PP.H4) was observed between elevated CRP and total body BMD on chromosome 2 amounting to 99.3% and revealing an exonic non-synonymous SNP in gene GCKR. Comparable findings were detected between CRP and arm BMD and spine BMD, however, at lower PPs (PP.H4: 10.5 and 10.4%, respectively; Supplementary Material, Table S4). Subsequently, a 90.0% PP of shared common SNPs was observed between CRP and UKB heel BMD in chromosome 9 mapped to an ncRNA intronic gene, ABO. Similarly, strong evidence of shared SNPs was seen between CRP and several OP-related phenotypes: spine BMD, UKB OP and total body BMD on chromosome 14 (88.0, 84.7 and 84.4%, respectively). All were linked to the exonic non-synonymous SNPs in gene SERPINA1 and the intergenic region of SERPINA2P/SERPINA1 (Table 4). Colocalization was shown between UKB heel BMD and CRP only in an intergenic region of chromosome 1 around LEPR/RN7SL854P and an intergenic region in chromosome 2 around AC093326.3/TMEM18 (PP.H4: 79.3 and 80.9%, respectively, Table 4). In addition, the genomic region of chromosome 22 demonstrated strong evidence of colocalizing SNPs between UKB heel BMD and CRP (76.4%; Table 4) and moderate evidence between OP fracture and CRP (73.0%; Supplementary Material, Table S3) corresponding to two genes, GTPBP1 and SUN2. Similar observations between CRP and hip and spine BMD were seen on chromosome 22 but at much lower PPs (Supplementary Material, Table S4).
Other colocalization regions with a moderate PP.H4 (between 50 and 75%) were found on chromosomes 8 and 1: between spine BMD and CRP in the intronic regions of RP11-115J16. 1   Results of common colocalizing genetic variants and their respective genes seen between OP phenotypes and CRP as determined by high PP of H4 (PP.H4 ≥ 75%) arising from commonly shared SNPs.
Material, Table S3). Colocalization results showing evidence of PP < 50% are given in Supplementary Material, Table S4.
To confirm the pleiotropic relationship between CRP and OP traits seen in Table 4, attributed to SNPs with high PP.H4, the PLeiotropic Analysis under COmposite null hypothesis (PLACO) analysis was carried out. The PLACO analysis revealed that each tested SNP was highly significantly pleiotropic between CRP and the designated OP variable, and the corresponding P-values ranged between 1.06E−09 and 1.17E−35 (Supplementary Material, Table S5).
In assessing the PP of H3 (PP.H3) between CRP and OP fracture and several BMD phenotypes in different study samples, 23 regions were revealed with high PP of at least 80%, though most achieving nearly 100% of colocalization arising from distinct SNPs. All of the results are provided in Supplementary Material, Table S6.
Gene enrichment results, denoted as genes of interest, were only reported for findings consistent with those in Table 4 and were attributed to high PP.H4, or genes with small P-values and PP.H3s approaching 100% (Supplementary Material, Table S7), since high PP.H3s may also be owing to spurious pleiotropy.
Notable for its reoccurrence, high PP.H3s between both UKB heel BMD and total BMD with CRP (PP.H3: 95.6% and 99.6%, respectively) revealed an exonic nonsynonymous SNP in the SERPINA1 gene around the intergenic region of SERPINA2P/SERPINA1 on chromosome 14 (Supplementary Material, Table S7).
We further examined potential colocalizing regions between RA and CRP. Two regions of interest on chromosome 6 were apparent. The first region is between base pairs 30 798 168 and 31 571 218, with a PP of distinct SNPs of 93.1% and linked to the intergenic region of HCP5/MICB (Supplementary Material, Table S8). The second region with PP.H3 of 69.2% was found between the base pairs 31 571 218 and 32 682 664 in the intergenic region of HLA-DRB1/HLA-DQA1.

Gene ontology
Since the colocalization analysis revealed a number of exonic SNPs that were significant in CRP and OP phenotypes, we examined whether they could be assigned to a set of genetic factors connected by a common pathway or share a gene ontology (GO) term denoting a common genetic topology. GO enrichment analysis was carried out for RA, CRP and all OP phenotypes. Next, the GO terms were mapped and annotated per phenotype and were further assessed using Fisher's exact test to generate top scores for GO terms present per selected class. Conversely, we identified whether at least one GO term included SERPINA1 and/or GCKR and observed that only GO:0005515 term describes both genes.
Fisher's test conducted for the GO molecular functions class using the classic algorithm revealed that GO:0005515 is significantly over-represented in SNPs associated with RA (P = 3.30E−10), CRP (P = 5.90E−07) and several OP traits, including UKB OP (P = 1.40E−15), total BMD (P = 2.40E−14), spine BMD (P = 1.40E−19), hip BMD (P = 1.20E−10) and UKB heel BMD (P = 3.10E−06), annotated over 13 000 genes (Supplementary Material, Table S9). Arm BMD and OP fracture variables of the GEnetic Factors for OSteoporosis (GEFOS) Consortium did not reveal significant associations with the GO term examined. Supplementary Material, Table S9 summarizes the ranking of the node, the respective annotations and Fisher's tests. GO:0005515 term was among the top three nodes for each variable that was significantly associated under the molecular functional class. The other two classes did not display statistically significant findings. GO:0005515 denotes protein binding and, therefore, demonstrates that RA, OP and CRP share a common genetic pathway that is involved in the molecular functional level of protein binding. The GO universe that explained this particular GO term is associated with both the SERPINA1 and GCKR genes (along with another ∼13 000 genes), which showed high colocalization between CRP and OP phenotypes and revealed exonic SNPs in both genes.
In addition to the two major exonic SNPs that were highly colocalized between CRP and OP phenotypes, we identified genes demonstrating colocalization between CRP and OP arising from common and distinct SNPs (Table 4 and Supplementary Material, Table S7) and were also a component of GO:0005515. We observed several genes that contribute to GO:0005515 from chromosomes 1, 2, 8, 11, 16, 18 and 22, including LEPR1, TMEM18, PPP1CB, SPDYA, TRMT61B, MSRA, TRPS1, C11orf49, FTO, PTPN2, LDLRAD4, RNMT, GTPBP1 and SUN2. We also examined genes in colocalization between RA and CRP and found two additional genes that are components of GO:0005515, including HLA-DRB1 and HLA-DQA1.

Discussion
The main aim of this study was to clarify the relationship between RA and OP and to additionally clarify whether increased CRP serves a role in the relationship and to what extent it is grounded on clinical data. The latter suggests that RA, an inf lammatory condition, necessarily contributes to raised CRP and this likely leads to OP through bone turnover mechanisms (17). These aims were accomplished by implementing MR, which ascertains whether a causal relationship is present and may explain the relationship between RA and OP. In addition, we determined which genes are involved in the pleiotropic relationship that appeared between CRP and OP.
MRE analyses between RA and varying OP phenotypes were suggestive of pleiotropy, which was previously shown by us (20) using polygenic risk score analysis in addition to colocalization and gene enrichment analyses and it is consistent with previous reports (18,21). Several loci have been discovered in relation to both phenotypes, particularly in the major histocompatibility complex (18,21). However, these pleiotropic regions appear to be polygenic with no clear contributing SNP. Here, we confirmed that the suggested common causal SNPs were pleiotropic through assessing PLACO.
Risk factors common to both RA and OP are well known (22) and include heightened RA disease activity, increased RA disease duration and the quality of life of RA subjects (12,23) and showed that RA subjects were at a greater risk for overall bone fracture and fragility fracture at the vertebral, hip, forearm and proximal humerus (23). However, recent studies highlighted the inflammatory component involved in the relationship between RA and OP, which is exacerbated with RA disease activity (14,22,24).
Elevated CRP is a response to systemic inflammation elicited by RA and is correlated with other RA disease markers, including bone and joint erosion, early morning stiffness and swollen joints and, in particular, bone deterioration (25). As such, MR results where RA causally predicts CRP are expected.
Importantly, however, our data consistently indicated that CRP-OP relationship is mostly attributable to horizontal pleiotropy. A shared genetic profile between CRP and OP traits was suggestive using PLACO and was subsequently confirmed using MRE, further suggesting that pleiotropy may explain the genetic relationship between CRP and OP.
The suspected pathophysiology, including CRP elevation in RA subjects, was suggested to arise from prompting RANKL as well as osteoclastogenesis, thereby eliciting bone resorption and subsequent BMD deterioration (25,26). The present results further clarify the strong and complex genetic relationship between RA and OP, which likely occurs via elevated CRP.
We implemented a colocalization analysis and functional annotation to investigate potential pleiotropic genes between CRP and OP focusing on non-synonymous SNPs, which are potentially deleterious for protein structure and function, in addition to demonstrating a high degree of colocalization. Two major nonsynonymous SNPs of genes, GCKR and SERPINA1, emerged consistently between CRP and various OP phenotypes. A previous study similarly demonstrated that hsCRP and OP phenotypes do not share a causal association; however, the study also identified GCKR and ABO as pleiotropic genes (19). Though this appears consistent with our findings, ABO was considered to be associated with interleukin-6, a different but closely linked inflammatory marker (19). GCKR codes for a glucokinase regulator, which is involved in metabolic modification and is associated with several phenotypes, including bone and inflammatory phenotypes (27). SERPINA1 is influenced by GCKR and is similarly involved in metabolic alteration and suggestively contributes to the pathway involving metabolic aberration through inflammatory mechanisms arising from alpha-1 antitrypsin deficiency (28). Taken together, these genes, which were prominent and pleiotropic across CRP and OP, suggest a potential pathway involving metabolic syndrome. Perhaps, these genes may be influenced by BMI and other latent variable(s), which may mediate the relationship between CRP and OP phenotypes. PPs were not strong in examining all OP phenotypes; however, differences among OP phenotypes are most probably attributable to smaller sample sizes, resulting in an under-powered assessment.
The HLA region was suggestive of pleiotropy between RA and CRP, specifically genes HCP5/MICB and region HLA-DRB1/HLA-DQA1. RA with high inf lammatory markers in the presence of HLA-DRB1 alleles was shown to have a greater risk of cardiovascular complications (29) and increased disease severity of RA (30).
To summarize the complex genetic architecture connecting RA and OP, we propose a simplified model that explains these relations. We speculate that while a portion of the genetic inf luence in RA and OP manifestation may be phenotype-specific, pleiotropic genes contribute polygenically to the development of RA and OP. In addition, elevated CRP results from systemic inf lammation, in this case from RA, and is pleiotropically associated with low BMD and OP fracture. The role of CRP further demonstrates the complex multifactorial nature of RA and RA-induced OP.

Limitations
We limited our genome-wide association study (GWAS) datasets to samples of European ancestry, thus our findings may not be applicable to the global population. The GEFOS OP fracture dataset includes data from the UKB, which may lead to bias and overinf lated results. While the use of GWAS datasets from differing data sources may contribute to reducing power, it allows for limiting bias in a two-sample MR study. The effect of behavior and lifestyle among other potential confounders cannot be evaluated through MR studies. We could not determine SNP independence prior to conducting MR analyses since the RA phenotype is binary in nature, so we were unable to extract SNPs that may serve as confounders. Despite this limitation, MR findings were similar across different BMD GWAS summary statistics. Different BMD samples were used, as well as two different BMD measurements with different sample sizes, which lends strength to the findings, while recognizing that the anatomical sites show differing susceptibility to OP. The SNPs identified in GWAS are not necessarily the true causal SNPs. Thus, fine mapping would be advantageous to determine the extent to which the SNPs alter or contribute to the phenotype.

Concluding remarks
The complex relationship between RA and OP inclusive of inf lammatory mechanisms should be clarified using molecular techniques. These three phenotypes share genetic factors that play a role in protein binding. Further investigation may reveal the intricacies of the shared molecular pathway potentially involved in the development of the debilitating musculoskeletal conditions and may lead to insight regarding pathways and potential therapeutics.

Study design
We conducted analyses using several large datasets and implemented modern genetic epidemiology techniques involving several consecutive analytical steps: (1) We established the relationship between the three variables of interest (RA, OP and CRP) independent of possible confounders. (2) We estimated causality between the study phenotypes, assuming that RA is the causal factor for OP and elevated CRP; CRP is causal factor for OP. (3) We examined the shared genomic regions for colocalization to reveal pleiotropy while conducting gene-set enrichment analysis to identify the potential common genes. (4) We confirmed pleiotropic effects for SNPs shared by CRP and OP traits, which showed high PP of being causal to both traits using a statistical pleiotropy test, PLACO. (5) We performed GO analysis to describe common genetic pathways.

Dataset for phenomenological analysis
A cross-sectional study was performed using data from the ELSA Wave 6 (31). The dataset consisted of 10 183 participants of European descent (4744 males and 5857 females), including 7731 individuals with measurements of hsCRP among other serological markers, collected during 2012 and 2013. Of these, 661 participants were diagnosed with RA and 819 with OP.

Datasets for genetic analysis
In total, 11 GWAS datasets were used. Summary statistics for GWAS conducted by the Neale Laboratory ( We collected additional GWAS summary statistics to conduct genetic analyses in a two-sample setting whenever applicable. We obtained another RA GWAS conducted by Okada et al. (33), which consisted of metaanalyzed data and was restricted to subjects of European ancestry (33) and comprised of 14 361 RA cases and 43 923 controls from 18 data sources with over 8 million imputed SNPs. Each RA subject was diagnosed in accordance with the revised 1987 criteria of the American Rheumatism Association or by a professional rheumatologist.
Finally, CRP GWAS summary statistics (quantitative variable) were obtained from the Cohorts for Heart and Aging Research (CHARGE) consortium (https://www. chargeconsortium.com/) Inf lammation Working Group (37) and it encompassed 204 402 participants gathered from 78 studies and 10 million SNPs.

Association between RA, OP and CRP
Statistical assessment of the ELSA phenotypes was conducted using R software environment (http://www.Rproject.org/) to determine the extent of the association of RA (as the dependent variable) with OP and CRP circulating levels (the independent variables) in addition to covariates age, sex and BMI. The glm function from the R stats package was used to conduct logistic regression analysis, implementing the pairwise deletion approach in case of missing data.

Mendelian randomization
To test for the potential causal effect between two phenotypes, we performed MR analysis using GWAS summary statistics and employed the MendelianRandomization package in R (38). Figure 1 describes the MR methodology.
Several MR approaches are available with the intent to satisfy different assumptions, including IVW and MRE (39). Except for MRE, we found most other MR approaches to correspond closely with IVW. Thus, we initially focused on MRE and IVW approaches.
The MRE approach is advantageous because it is a robust approach that assumes and detects pleiotropy (40). While it is more conservative, it provides a more pronounced bias in a two-sample setting with overlap or one-sample setting, as it is prone to type I error. A component of the MRE approach determines instrument validity, the I 2 Gx sensitivity statistic, which was used to measure the instrumental variable dilution or bias (41), and it was suggested to be at least 90% in a two-sample analysis (40).
IVs, or genetic variants of the causal trait, were chosen following the compilation of the GWAS summary statistics of the exposure and outcome phenotypes. The P-value threshold of the effect of the exposure variable was restricted to P < 5E−08 regardless of the P-value attributed to the outcome phenotype. Subsequently, linkage disequilibrium (LD) clumping was performed to remove SNPs in LD by utilizing the ld_clump function available in the MRCIEU/ieugwasr R package (42) with a parameter of clumping R 2 = 0.01. This procedure along with the harmonization of the GWAS summary statistics of the two traits being compared were performed (43). In assessing RA/CRP as a risk factor for BMD, IVs were restricted to associated beta estimates of opposing signs between the two traits; conversely, in assessing whether RA is a risk factor of OP binary phenotypes or CRP or if CRP is a risk factor of OP binary phenotypes, IVs were restricted to associated beta estimates of the same signs for both traits.
We chose the meta-analyzed RA GWAS and CHARGE CRP GWAS instead of their UKB GWAS counterparts since we preferred to perform two-sample MR tests while avoiding sample overlap.

Colocalization analysis and gene enrichment analysis
Pleiotropy may explain the genetic relationship between OP phenotypes and CRP, so we conducted colocalization analysis to further decipher the regions that may contain shared SNPs (44). We examined SNPs that appeared to colocalize either through mutual common SNPs or distinct SNPs within the same gene. The colco.abf function from the coloc package in R was used (45). This colocalization method is based on Bayesian modeling and includes testing of the PPs of five mutually exclusive hypotheses (44). We scanned for SNPs achieving GWAS significance threshold for CRP and defined the corresponding genomic region for colocalization analysis by LD Detect blocks (46). LD detect regions were selected as boundaries of prominent SNPs, which include generated independent LD blocs (46). These boundaries are the regions where a higher frequency of recombination takes place. Subsequently, the PPs of each of the five hypotheses were evaluated in relation to OP phenotypes.
We focused on Hypothesis 4 (H4) and Hypothesis 3 (H3), which assume the existence of causal SNP(s) attributed to both traits in a given genomic region. However, H4 estimates the PP of shared common SNPs affecting both phenotypes, whereas H3 estimates the PP of distinct SNPs within the same region, although common SNPs may also be present (44). The PP was considered as having strong evidence if it was at least 75%, while a PP of at least 50% was cautiously suggestive of evidence of different causal SNPs mapped to the designated genomic region (45).
The coloc.abf function additionally generates an output with SNP.PP.H4, which includes the PP of each SNP being the causally conditional shared SNP, corresponding to H4, between the two examined traits. Following a PP of H4 > 50%, we focused on the SNP with the highest PP to explain the suggested colocalized relationship, while prioritizing and/or also reporting exonic SNPs.
Next, we used the functional mapping and annotation GWAS platform (47) for gene enrichment analysis to define the genes in common between the two phenotypes, while focusing on SNPs with the smallest P-value for both traits.

Statistical presence of pleiotropy
PLACO was implemented to determine evidence of pleiotropy of the selected shared SNP between two traits (48). SNPs were identified through the colocalization analysis and were subsequently assessed using PLACO to confirm the pleiotropic relationship between the traits, which was suggestively attributed to the SNP. A significant P-value from PLACO is suggestive of a common genetic framework between the traits. PLACO utilizes GWAS summary statistics to generate the pleiotropic test, which is subsequently tested considering the Z scores (beta/effect) of each trait in relation to the SNP in question (48). PLACO is limited to common causal SNPs and was, therefore, only applied to SNPs showing high PP of H4 from the colocalization analysis.

Gene ontology
GO enrichment analysis was carried out to decipher at least one potential model that would describe genetic topology of one GO term or node and can be calculated under three classes: biological process, molecular function and cellular component (49). Each GO term explains a cluster of genes that contribute to the process, function or component examined. Thus, we used the topGO R package in Bioconductor (50) to identify highly associated GO terms under the three classes with RA, OP and CRP variables. Algorithmic parameters were set to 'classic', and Fisher's exact test statistic was assessed.

Data availability
Data were acquired from databases and consortia ,such as GEFOS, UKB, Neale Laboratory, CHARGE and ELSA. Access to these data may be provided from the aforementioned sources.

Supplementary Material
Supplementary Material is available at HMG online.