Cross-Disorder Genomics Data Analysis Elucidates a Shared Genetic Basis Between Major Depression and Osteoarthritis Pain

Osteoarthritis (OA) and major depression (MD) are two debilitating disorders that frequently co-occur and affect millions of the elderly each year. Despite the greater symptom severity, poorer clinical outcomes, and increased mortality of the comorbid conditions, we have a limited understanding of their etiologic relationships. In this study, we conducted the first cross-disorder investigations of OA and MD, using genome-wide association data representing over 247K cases and 475K controls. Along with significant positive genome-wide genetic correlations (rg = 0.299 ± 0.026, p = 9.10 × 10–31), Mendelian randomization (MR) analysis identified a bidirectional causal effect between OA and MD (βOA→MD = 0.09, SE = 0.02, z-score p-value < 1.02 × 10–5; βMD→OA = 0.19, SE = 0.026, p < 2.67 × 10–13), indicating genetic variants affecting OA risk are, in part, shared with those influencing MD risk. Cross-disorder meta-analysis of OA and MD identified 56 genomic risk loci (Pmeta ≤ 5 × 10–8), which show heightened expression of the associated genes in the brain and pituitary. Gene-set enrichment analysis highlighted “mechanosensory behavior” genes (GO:0007638; Pgene_set = 2.45 × 10–8) as potential biological mechanisms that simultaneously increase susceptibility to these mental and physical health conditions. Taken together, these findings show that OA and MD share common genetic risk mechanisms, one of which centers on the neural response to the sensation of mechanical stimulus. Further investigation is warranted to elaborate the etiologic mechanisms of the pleiotropic risk genes, as well as to develop early intervention and integrative clinical care of these serious conditions that disproportionally affect the aging population.


INTRODUCTION
Osteoarthritis (OA) is the most frequent form of arthritis, affecting over 32.5 million adults in the United States [Centers for Disease Control and Prevention (CDC), 2020]. The cardinal symptom of OA is debilitating and chronic pain that affects synovial joints, including the knees, hips, and back, often aggravated by stiffness, swelling, and reduced mobility. Patients with OA often struggle with major depression (MD) and anxiety, linked with greater burden of pain severity, functional disability, and increased mortality (Fuller-Thomson et al., 2016;Sharma et al., 2016;Sambamoorthi et al., 2017); a longitudinal study using the United States National Institute of Health OA Initiative data has reported higher odds of developing depressive symptoms in individuals with hip, knee, or multi-site OA [Odds Ratio (OR): 1.43 ∼ 1.72] (Veronese et al., 2016). Similarly, Akintayo et al. (2019) reported that 42% of OA patients suffer from MD (OR: 2.49).
The higher comorbidity of MD and OA has been interpreted in two major ways (Figure 1). Under the psychological, environmental risk model, chronic pain has been identified as a major mediating risk factor leading to the increased risk of depression in OA patients (Mammen and Faulkner, 2013;Veronese et al., 2016). Reduced physical activity and few social interactions are also common in people suffering from OA (Wallis et al., 2013) and especially those affected in the lower parts of the body such as the knee and hip (Stubbs et al., 2015). In contrast, clinical studies have suggested that OA and MD may share some common etiologic mechanisms related to stress, inflammation, and immune responses. MD has been associated with increased inflammatory activation of the immune systems in both the periphery and the central nervous systems (Lee and Giuliani, 2019), while OA patients exhibit higher inflammatory markers (Bonnet and Walsh, 2004) and cortisol levels (Carlesso et al., 2016), indicating altered immune systems. Neurobiologic studies have indicated that the autonomic nervous system has a role in joint homeostasis and OA pathogenesis, suggesting pain sensitization (Clauw and Hassett, 2017;Syx et al., 2018) and neuronal impairment (McDougall, 2019) as potential risk factors connecting OA and MD.
In this study, we examined the etiologic relationships between OA, pain, and MD using cross-disorder genome-wide genetic data analysis. Specifically, we aimed to clarify whether OA and MD share common genetic risk mechanisms that increase susceptibility to both conditions (known as biological pleiotropy), or genetic risk influencing OA and pain mediates increasing risk to MD (known as mediated pleiotropy). To this aim, we applied a series of pleiotropy analyses, using by far the largest genomewide genetic variation data of the two disorders, representing a sample of 77 K OA cases, 171 K MD cases, and 450 K controls. First, we measured genome-wide genetic correlations between OA and MD and inspect their potential causal relationships using Mendelian randomization (MR) analyses. Further, we performed a cross-disorder genome-wide association study (GWAS) to identify specific risk loci, genes, and biologic mechanisms that may influence the two health conditions. Lastly, we investigated five drug-gene interaction databases to clarify clinical implications of the shared etiology. To the best of our knowledge, this is the first investigation of this kind.

GWAS Datasets
We obtained the latest GWAS summary statistics for OA (Tachmazidou et al., 2019), which investigated a total of A B FIGURE 1 | Two distinct but not mutually exclusive risk models for OA, pain, and MD. (A) Psychological, environmental risk model attributes various negative health outcomes of OA as a major contributing factor to increasing risk to MD. (B) Shared biological risk model proposes that OA and MD are affected by some common biological mechanisms related to stress response, chronic inflammation, and articular/sensory neuronal dysfunction. 455,221 individuals (77,052 OA cases and 378,169 controls). We obtained GWAS summary statistics for MD from the Psychiatric Genomics Consortium (PGC) (Howard et al., 2019). This study was based on a total of 170,756 MD cases and 329,443 controls. Table 1 summarizes the statistics of the two GWAS datasets. For further details, we refer to the original publications (Howard et al., 2019;Tachmazidou et al., 2019).

Quality Control and Standardization of GWAS Datasets
We applied unified quality control (QC) of the two GWAS summary statistics. Original datasets included 8.91 M SNPs for OA and 8.48 M SNPs for MD. First, chromosome coordinates were standardized to the reference genome GRCh37/hg19. We removed SNPs that were rare (minor allele frequency < 1%), non-biallelic, lack of rs IDs, not present in the 1,000 Genomes Project Phase 3 European samples, and had incompatible alleles. We also removed palindromic SNPs, with allele frequency differences less than 15%. After QC, 7.85 M and 7.75 M SNPs retained for OA and MD, respectively. We restricted analysis to autosomal chromosomes.

Estimation of Polygenicity, Discoverability, and Genetic Correlation
To compare the genetic architecture of MD and OA, we estimated the polygenicity (defined as the fraction of independent causal SNPs, or π) and discoverability (the variance of causal SNP effect sizes, or σ 2 β ) using a univariate Gaussian mixture model implemented in MiXeR (Holland et al., 2020). Genomewide genetic correlation was estimated using LDSC (linkage disequilibrium score regression) with pre-computed European LD scores (Bulik-Sullivan et al., 2015). For measuring a liability-based SNP heritability, we used a population prevalence compiled from the literature: 13.46% for OA (Tachmazidou et al., 2019) and 15% for MD (Lee et al., 2013). We used the GCTA-GREML power calculator (Visscher et al., 2014) to calculate the statistical power of estimating genetic correlation of at least 0.2 between OA and MD with default setting (α = 0.05). We ran the sign-concordance test between OA and MD using the binomial test in R, assessing how often genome-wide significant SNPs for one phenotype show the same directional effect with the other phenotype. We expect sign concordance of 50% between any random phenotypes, and the estimated binomial probability represents how unlikely it is to observe the assessed concordance between two phenotypes by chance.

Mendelian Randomization
Mendelian randomization examines whether relationships between an exposure and an outcome are causal (Minca et al., 2020). We performed bi-directional, two-sample MR analysis between OA and MD using GSMR (Generalized Summary-data-based Mendelian Randomization) (Zhu et al., 2018). This multi-step method uses genome-wide significant SNPs for individual traits (p ≤ 5 × 10 −8 ) as instrumental variables in each unidirectional test. GSMR takes into account linkage disequilibrium (LD) between genetic variants and sample variance in genetic effects on exposure and outcome for each SNP. We estimate the causal effect (b xy ) of the exposure (b zx ) on the outcome (b zy ) as:b xy =b zy /b zx We used default settings for the method with the 1,000 genomesbased European samples for estimating LD. To avoid potential bias and test-statistic inflation due to pleiotropy, we used the HEIDI-outlier (heterogeneity in dependent instruments) method to detect putative pleiotropic SNPs.

Cross-Disorder Meta-Analysis
We conducted meta-analysis of OA and MD using METAL (Willer et al., 2010). This method adjusts potential sample overlaps by first estimating the covariance of the Z-score distribution and assigning proper weights while summing Z-scores. Optimal weights were previously represented by Lin and Sullivan's work (Lin and Sullivan, 2009): where w k is the weight for the k th study, e is a K × 1 vector of 1's and is the estimated covariance matrix of (Z 1 , . . . Z k ). To estimate the covariance matrix, Z-scores were truncated to remove non-null effects using the absolute Z-score threshold of 1. We then meta-analyzed the two datasets using: wherer kl is the estimated correlation between Z-scores of the k th and l th studies under the null. A total of 7,433,563 SNPs were meta-analyzed. We used LDSC and genomic inflation factors to assess inflation due to potential confounding factors of the meta-analysis results (Clayton et al., 2005). Number of cases and controls, total sample size, population and sample prevalence, and PubMed ID for two GWAS for major depression (MD) and osteoarthritis (OA).

Functional Annotation of Risk Loci Using FUMA
With summary-level association statistic results from crossdisorder meta-analysis as input, we performed functional annotation of risk loci using FUMA (Watanabe et al., 2017). We used default settings, while excluding the MHC region (chr6:25-35 Mb) for annotations and gene-set enrichment analysis. Gene mapping of SNPs was based on functional consequences on genes from positional, eQTL, and chromatin interaction mapping data in the brain. For annotations, we included protein coding and non-coding RNAs.

Gene-Based Association Test and Gene-Set Enrichment Analysis
Gene-set enrichment analysis was performed using MAGMA v1.6 (20 kb gene window) (de Leeuw et al., 2015) and the MSigDB v7.0 Gene Ontology (GO) data (Liberzon et al., 2011). We also ran gene-level analysis using cross-disorder meta-analysis summary statistics. Significant genes were designated using a p-value threshold P meta ≤ 2.64 × 10 −6 (=0.05/18,939 genes) for adjusting multiple testing of 18,939 genes at alpha = 0.05. Among the identified genes with significant association, we prioritized the ones that show suggestive association with both OA and MD (P OA ≤ 1 × 10 −4 , P MD ≤ 1 × 10 −4 ). Tissue-specific gene expression analysis was conducted with GTEX data (v8) (GTEx Consortium, 2020).

Regional Brain and Cell-Type-Specific Expression Analysis
We examined whether shared OA/MD risk genes show celltype specific gene expression using six brain cell type RNA sequencing data: astrocyte, endothelial, microglia, neuron, oligodendrocyte, and OPC (Darmanis et al., 2015). To designate significant RNA expression for each cell type, we used two criteria suggested in the original publication: mean expression threshold of at least five and a Z-score of 2. We also investigated spatiotemporal gene expression using the Human Brain Transcriptome (Kang et al., 2011).

Search of Drug-Gene Interaction Databases
We investigated drug-gene interactions related to OA and MD using five databases: DGIdb v3.0 (Cotto et al., 2017), DrugBank v5.0 (Wishart et al., 2017), PharmGKB (Whirl-Carrillo et al., 2012), STITCH v5.0 (Szklarczyk et al., 2015), and TTD 2020 (Wang et al., 2019). First, we compiled a list of OA-and MD-associated genes identified from disease-specific gene-based association analyses, shared OA/MD risk genes obtained from cross-disorder meta-analysis, and member genes of the gene-sets for which we identified significant cross-phenotype association in gene-set enrichment analyses described in the previous sections. Additionally, we included a list of OA and MD-associated genes compiled from the original GWAS studies (Howard et al., 2019;Tachmazidou et al., 2019).

Osteoarthritis and Major Depression Share Significant Genome-Wide Genetic Correlations
We have assembled GWAS datasets of OA (N = 77,052 cases and 378,169 controls) and MD (N = 170,756 cases and 329,443 controls), which represent by far the largest publicly available studies of the two disorders ( Table 1). Statistical power analysis indicated >99% power for detecting SNP-based heritability and genetic correlation in our datasets (Supplementary Table 1). Significant SNP-based heritability was estimated for both disorders on the liability scale: OA (h 2 = 9.4%, SE = 0.004) and MD (h 2 = 8.5%, SE = 0.003). The disorders were also highly polygenic, involving thousands of causal risk variants; approximately 8.7 K variants were estimated to causally influence OA and 13.7 K influence MD. Sign tests of genome-wide significant SNPs indicated that 67 out of 90 genome-wide significant SNPs for MD share the same directional effect with OA (binomial sign test P MD → OA = 3.80 × 10 −6 ), while 33 out of 51 OA index SNPs share direction with MD (P OA → MD = 4.89 × 10 −2 ). We also identified a substantial and statistically significant level of positive genome-wide genetic correlation between the two disorders (r g = 0.30 ± 0.026, p = 9.10 × 10 −31 ).

Mendelian Randomization Predicts Shared Genetic Etiology Between OA and MD
To infer potential causality in the relationship between OA and MD, we performed bi-directional MR analyses using a GSMR method. As input, we used 39 and 64 linkage disequilibrium (LD)-independent SNPs that show genome-wide significant association with OA and MD, respectively. We ran GSMR under a default-setting while removing SNPs displaying horizontal pleiotropy (HEIDI outlier p-value < 0.01). Figure 2 illustrates the bidirectional MR analysis results. When using MD as the risk factor and OA as the outcome, we found an effect size β MD→OA of 0.192 (SE = 0.026, p < 2.67 × 10 −13 ). The reverse-directional MR test estimated a significant but more modest effect of OA on MD (β OA→MD = 0.09 ± 0.02, p < 1.02 × 10 −5 ).

Cross Disorder Meta-Analysis Highlights a Role of Mechanosensory Behaviors
To elucidate shared genetic risk mechanisms, we conducted a cross-disorder meta-analysis of OA and MD. Figure 3 summarizes the genome-wide meta-analysis results of 7.43 M SNPs. LDSC analysis indicates that the observed inflation in the QQ-plot is due to underlying polygenic signals of the two disorders rather than confounding factors (LDSC intercept = 0.9816 ± 0.01). Genomic inflation factor λ 1,000 was 1.001. After LD-based pruning, we identified 56 genomic risk loci represented by 63 lead SNPs (SNP-based association P meta ≤ 5 × 10 −8 ; Supplementary Table 2). The identified loci were annotated using various functional genomics data FIGURE 2 | Bidirectional Mendelian randomization analysis results. The X-axis represents the estimated effect of the SNPs associated with the exposure (b zx ), while the Y-axis represents the estimated effect of the same SNPs on outcome (b zy ). Each cross marks independent instruments (i.e., SNPs) used to test for causality. The ratio,b zy /b zx is an estimate of the mediation effect of the exposure on the outcome assuming that under a causal model, if an exposure has a causal effect on an outcome, all instruments causally associated with the exposure will have an effect on the outcome expected to be identical.
We examined cell type-specific gene expression of the 43 OA/MD risk genes. More than a half of the OA/MD risk genes showed highly enriched expression in neurons (24 out of 43 genes), followed by astrocytes (nine out of 43) (Supplementary Table 7). Many genes showed significant expression in all six brain regions we assessed (cerebellar cortex, mediodorsal nucleus of the thalamus, striatum, amygdala, hippocampus, and neocortex) (Supplementary Table 8). More than three quarters of the brain-expressed genes were also significantly expressed both pre-and postnatally, suggesting their significant role throughout the lifespan.

DISCUSSION
While estimated prevalence varies depending on the sample size, population, and measurement tools, epidemiology studies have consistently reported the increased comorbidity between OA and MD (Gandhi et al., 2015;Hooten, 2016;Tsuji et al., 2019). The health burden of comorbid OA and MD is surprisingly harmful, spanning from greater symptom severity (de Heer et al., 2014) to poorer treatment outcomes (Roughan et al., 2020) to increased suicide attempts (Fuller-Thomson et al., 2016). Using the GWAS datasets of over 697 K individuals, our study shows for the first time that the increased comorbidity of OA and MD is, at least in part, due to a shared genetic etiology.
Cross-disorder meta-analysis of OA and MD revealed important insights into potential genetic mechanisms shared between the two conditions. First, we identified 63 LDindependent genome-wide significant SNPs associated with OA and MD. Implicated risk genes showed heightened expression in the brain and pituitary but not in other tissue types, suggesting that the shared genetic mechanisms between the two conditions are involved in the central nervous system and the endocrine system. Specific to the brain, we found significant enrichment of the shared risk genes in multiple brain regions, including the frontal cortex, the anterior cingulate cortex (ACC), and amygdala. Neuroimaging studies have reported a number of brain regions that play an important role in depression and pain perception, many of which converge to the ACC and amygdala (Sellmeijer et al., 2018;Serafini et al., 2020). In fact, Cottam et al. (2016) found that OA pain intensity correlates with blood flow in brain regions including the ACC and amygdala. Similarly, Brown et al. (2020) observed altered amygdala connection density in depressed patients. Our findings support the implication of the ACC and amygdala, along with other brain regions, in the shared genetic basis connecting OA, pain, and MD. Top 20 shared OA/MD risk genes out of 42 statistically significant shared risk genes (gene-based P meta ≤ 2.64 × 10 −6 , P OA ≤ 1 × 10 −4 , P MD ≤ 1 × 10 −4 ). Table  includes gene symbol, locus, full name, and p-values for individual analyses (OA, MD) and meta-analysis.
FIGURE 4 | Venn diagram summarizing the overlap of gene-based analysis results from OA, MD, and cross-disorder GWAS. The numbers in the Venn diagram represent the number of shared and/or distinct genes identified as statistically significant in three gene-based association analyses after Bonferroni correction (p < 2.64e-06). The total number of OA-, MD-, and cross-disorder associated genes are listed in the parentheses. There are eight genes for which all three gene-based analysis showed genome-wide significant association (listed in yellow box). In the right blue box, we list 43 genes, which showed genome-wide significant association in meta-analysis-based gene association analysis, and suggestive association at p < 1e-04 in OA-and MD-specific association.
Secondly, our analysis revealed for the first time a statistically significant association of "mechanosensory behavior" genes with OA and MD. This gene set comprises 18 genes, and our GWAS data had coverage of 15 genes. Mechanosensory behavior genes are involved in mechanical stimuli, such as emotional processes and behaviors related to pain and mechanosensation. Interestingly, autism spectrum disorder (ASD) mouse models with dysfunctional peripheral mechanosensory neurons have demonstrated touch hypersensitivity and altered tactile discrimination, contributing to difficulties in social interaction and anxiety symptoms (Orefice et al., 2016). Mechanosensation may also play a role in exacerbating pain through chronic stress, increasing synaptic efficiency in the amygdala, which combines nociceptive and affective information (Li et al., 2017). Additionally, certain mechanosensory nerves have been found to trigger the release of oxytocin in response to low intensity cutaneous stimulation and may serve to reward physical contact with others, reduce physiological arousal, and inhibit pain (Fingleton et al., 2015;Clauw and Hassett, 2017;Walker et al., 2017).
Third, although our cross-disorder meta-analysis did not find specific enrichment of genes involved in neuroinflammation, investigations of drug-gene interaction data suggest a potential biologic link relating OA and MD through the Hypothalamic-Pituitary-Adrenal (HPA) axis. Notably, CRHR1, the OA risk gene identified in the latest GWAS by Tachmazidou et al. (2019), plays a key role in the HPA axis. This gene is related to longterm depression and neuroactive ligand-receptor interaction in the brain (Schierloh et al., 2007), while interacting with 19 drugs indicated for MD and other mood disorders, anxiety disorders, and alcohol and drug dependence (Liu et al., 2007;Geng et al., 2014). The role of the HPA axis and implication of neuroinflammation in OA and psychosocial conditions, including anxiety and depression, has been well established (Loggia et al., 2015;Albrecht et al., 2018Albrecht et al., , 2019a. Carlesso et al. (2016) reported higher cortisol production in OA, which may be associated with deleterious effects on bone density and muscle mass (Sato et al., 2018). Similarly, elevated cortisol levels are associated with persistent depressive symptoms, with strong links to somatic symptoms like poor sleep and low energy (Iob et al., 2020). In fact, symptoms of cortisol dysfunction include fatigue, depression, and pain (Hannibal and Bishop, 2014), and patients with MD also frequently experience unexplained painful physical symptoms (UPPS) (Jaracz et al., 2016).
Lastly, we identified several OA-specific risk genes that interact with various psychiatric medications, suggesting the importance of integrative clinical care. The OA risk gene CYP1A2 (P OA = 2.47 × 10 −10 , P MD = 0.782) encodes a cytochrome P450 superfamily enzyme that catalyzes reactions involved in the metabolism of drugs, chemicals, and various substrates (Kapelyukh et al., 2019). CYP1A2 interacts with multiple psychiatric medications (e.g., clozapine, escitalopram, and paroxetine) indicated for MD, schizophrenia, anxiety disorders, and post-traumatic stress disorder (Kohlrausch et al., 2013;Kuo et al., 2013). Another OA risk gene ITIH3 (P OA = 2.59 × 10 −7 , P MD = 5.27 × 10 −6 ) interacts with antipsychotics clozapine (Brandl et al., 2016) and has been linked to depression (Amare et al., 2020), schizophrenia (Ikeda et al., 2018), and ASD (Xie et al., 2020), among other disorders. This gene has also shown genome-wide significant and independently replicated SNP-based associations with anxiety (Nagel et al., 2018) and risk-taking behavior (Karlsson Linnér et al., 2019). Our drug-gene interaction search also revealed a link between the OA-risk gene SLC39A8 (P OA = 2.14 × 10 −7 , P MD = 0.080) and risk of developing alcoholism . SLC39A8 encodes a ubiquitously expressed metal cation transporter (Nebert and Liu, 2019). This gene exhibits extensive pleiotropy and is associated with a range of diseases including schizophrenia, Parkinson's disease, and inflammatory bowel disease (Nebert and Liu, 2019). DPEP1 is another OA risk gene (P OA = 1.32 × 10 −8 , P MD = 0.643), encoding an enzyme that hydrolyzes a range of substrates and participates in liver and lung neutrophil recruitment, interacts with glutamic acid, which may ameliorate symptoms of schizophrenia and help control alcoholism (Hammarström, 1983;Murphy and Gijón, 2007;Vance and Vance, 2008;Buczynski et al., 2009;Claesson, 2009;Choudhury et al., 2019). In addition, DPEP1 has also shown associations with mathematical ability (Lee et al., 2018) and white matter microstructure (Zhao et al., 2019). Further investigation is warranted to decipher the functional role of OA-risk genes in the brain, as well as their interactions with psychiatric medications.
There are several limitations of the present study that we recognize. First, our bi-directional MR analysis indicates that there is no single directional causal effect between OA and MD. It is important to note that like all statistical methods, MR has key assumptions, such as no association between instrumental variables and any potential confounder (e.g., assortative mating). Within-family MR designs will be helpful to clarify issues related to potential uncontrolled confounders (Brumpton et al., 2020). Secondly, shared genetic effects can occur due to genuine pleiotropy, but there is also the possibility that LD among genetic loci or population stratifications (Haycock et al., 2016) may induce artificial or spurious pleiotropy. Robust genome-wide genetic correlations and inspection of our cross-disorder analysis data, however, overrule that these phenomena are the major factors driving the shared genetic relationships between OA and MD. Third, it is possible that comorbid cases of MD and OA may have contributed to pleiotropy we detected. This is a challenging issue for any large-scale genomic studies where comprehensive data on disorder comorbidity are not available. On the other hand, when multiple phenotypes overlap in their genetic risk architecture, we would expect increased comorbidity between patients with corresponding disorders. Thus, to the extent that comorbidity is related to genetic overlap, it is not confound but a reflection of the pleiotropy we are detecting. It will be important for future studies to investigate the impact of comorbidity in pleiotropy analysis. Fourth, stringent QC and standardization is necessary for cross-disorder analysis of two separately conducted GWAS summary statistics. We thus applied strict QC procedures that excluded about 5% of SNPs. This was necessary but may have resulted in the loss of power for subsequent gene and pathwaybased analysis. OA and MD datasets also have different sample sizes, which may affect statistical power to detect pleiotropic risk loci. Fifth, both OA and MD have sex and age-dependent diseaseprevalence and symptom severity. While we could not evaluate the relevance of sex and age in the shared genetic basis of the two disorders due to a lack of data, it is an important direction of our future research. Lastly, we analyzed samples of European ancestry, so our results may not generalize to other populations. In fact, the genetic architecture of OA and MD in other populations remains largely unexplored. Concerted efforts are warranted for collecting and interrogating non-European cohort samples, with ultimate aims of reducing the health disparity.
In summary, we have performed a cross-disorder genomewide data analysis to elaborate the genetic relationships between OA and MD, two chronic and disabling conditions that affect millions of people worldwide. Our findings demonstrate that OA and MD share common genetic risk mechanisms, one of which centers on the neural response to the sensation of mechanical stimulus. Further investigation is warranted to elaborate the etiologic mechanisms of the mechanosensory behavior genes, as well as to develop effective strategies for early intervention and integrative clinical care of these serious health conditions.

AUTHOR CONTRIBUTIONS
PL conceived the project idea. PL, J-YJ, JS, ML, and MF developed the study. SB, J-YJ, NN, and MS ran the data analysis. All authors contributed to data interpretation, manuscript drafting and editing, and agreed to the submission of this manuscript.

FUNDING
For this research, PL was supported by National Institute of Mental Health (NIMH) R00 MH101367 and R01 MH119243. ML was supported by National Institute of Neurological Disorders and Stroke (NINDS) R01 NS094306.

ACKNOWLEDGMENTS
We would like to thank the participants who donated DNAs for various genetics research of common and complex human diseases. We would like to appreciate the clinical and scientific teams that processed, analyzed, and publicly shared the summary statistics of individual disease GWAS datasets, including the Psychiatric Genomics Consortium, United Kingdom Biobank, and the Understanding the Society Study. We would also like to thank the investigators for their dedications in genetics research and the open data science policy. Statistical analyses were carried out on the partner's research computing cluster servers and highperformance computing clusters hosted by the Broad Institute of MIT and Harvard. We would also further like to thank Christel Bastida for her constructive comments and the reviewers of the journal for insightful suggestions and editorial support.