2.1 Data sources for overall OA, knee OA, and HIP OA
We obtained pooled-level GWAS data for overall OA from a meta-analysis of European ancestry participants. OA patients were defined based on UK Biobank data, a large biomedical database and biobank that collects genotype, phenotype, lifestyle data, medical records, and biological samples from 500,000 volunteers in the UK. Individuals were considered to have OA if they had either (1) a physician diagnosis of OA, (2) joint replacement surgery, or (3) joint pain or stiffness and were aged ≥45 years. The overall OA dataset included 327,918 European individuals, comprising 30,727 OA patients and 297,191 controls. [12]
The knee OA and hip OA data were obtained from two cohorts, UK Biobank and arcOGEN, with patient inclusion criteria being hospital-documented or self-reported osteoarthritis of the knee or hip and no hospital-documented or self-reported osteoarthritis of the other joint. The knee OA dataset included 402,994 European individuals, consisting of 24,955 knee OA patients and 378,039 controls. The hip OA dataset included 393,743 European individuals, of which 15,704 were hip OA patients and 378,039 were controls. [13]
2.2 Data source for brain cortex surficial area and cortex thickness
Magnetic resonance imaging (MRI) data from 51,665 individuals were used in this study, obtained from 60 cohorts consisting of 33,992 individuals of European origin and 17,673 individuals of non-European origin (enrolled cohorts were listed in Supplementary Table S4). Only individuals of European origin were analyzed for this study. The MRI data underwent quality control and preprocessing to exclude data with significant abnormalities or missing values, as well as individuals with duplicates or relatives. For the GWAS data on cortical structures, we examined surface area and mean thickness of the whole cerebral cortex, as well as surface area and thickness of 34 cortical regions with known functional specialization. These 34 brain regions were distinguished using a gradient and boundary algorithm that detects points of change in multimodal data on the cerebral cortical surface. We only analyzed the 34 specific brain regions that were corrected for global surface area or thickness as covariates, in order to exclude the confounding effect of global on regional factors and to more accurately identify genetic variants affecting specific brain regions.[14]
Overall, knee OA, and hip OA data were used to determine thickness (TH) and surface area (SA) for the entire cerebral cortex. We analyzed TH and SA of 34 cortical regions with and without global brain-weighted estimates, yielding 138 results. Data that included globally weighted estimates indicated that SA and TH in specific regions were related to SA and TH in the whole brain, whereas data without globally weighted estimates indicated that SA and TH measurements in specific regions were not related to whole brain SA and TH.
2.3. Selection of genetic instruments
To explore the causal relationship between osteoarthritis and cerebral cortical structures, we selected three genetic data sets for osteoarthritis: overall OA, hip OA, and knee OA. To ensure the validity of the instrumental variables, we checked the three basic model assumptions of MR [15]. Firstly, the genetic variants must be strongly associated with the target exposure. Secondly, the genetic variants should not be associated with any confounders of the exposure-outcome relationship. Lastly, the genetic variants can only affect the outcome through the target exposure and not through other pathways. These assumptions ensure that the genetic variants can be used as valid instrumental variables to simulate a randomized controlled trial and thus infer a causal relationship between exposure and outcome.
Initially, we selected independent SNPs that were strongly associated with brain cortical structure (p-values < 5e-08), but the number of such SNPs was too low for our analysis. Thus, we relaxed the threshold to 5e-06 to screen our instrumental variables. We then performed a clustering procedure with R2<0.001 and window size=10000kb on European ancestral individuals from the 1000 Genomes Project to exclude SNPs with strong linkage disequilibrium (LD). Subsequently, we removed SNPs that had a weaker association with brain cortical structure (p-value > 5e-06). We applied MR Pleiotropy Residual Sum and Outlier (MR PRESSO) to correct for the presence of horizontal pleiotropy. This method uses genome-wide summary statistics to perform linear regression analysis on multiple instrumental variables, then calculates the residual sum of squares and compares it with the data from random simulations to detect the presence of horizontal pleiotropy [16]. After finding horizontal pleiotropy, we improved the accuracy of causal inference by correcting for the effect of outliers. The specific analysis process can be seen in the figure 1.
2.4. Mendelian randomization analyses
Three MR analysis methods were utilized in our study to address heterogeneity and pleiotropy issues: random-effect inverse-variance weighted (IVW), MR Egger, and weighted median. We excluded SNPs associated with brain cortical structures and outliers detected by MR-PRESSO before using IVW as the primary outcome. IVW is the most established and widely used method because it provides the most efficient and accurate causal estimates if there is no horizontal pleiotropy of the instruments [17]. It regresses the exposure effect on the outcome effect without an intercept term, using the inverse of the outcome variance as the weight. The instruments must satisfy the three core assumptions of MR, i.e., no horizontal pleiotropy, to obtain correct causal estimates. Weighted median is another method that only requires that most of the instruments (more than 50%) satisfy the three core assumptions of MR to obtain correct causal estimates. MR Egger is a method that can detect and correct for horizontal pleiotropy and the intercept term can be used to assess its magnitude. We further assessed horizontal pleiotropy for significant estimates using MR Egger intercept test and funnel plot analysis. Cochran's Q test was used to detect heterogeneity, and funnel plots were used to visualize possible directional pleiotropy. To exclude potential confounding factors, such as body mass index, obesity, waist-to-hip ratio, smoking, etc., we searched PhenoScanner and removed any SNPs that showed genome-wide association with any of these factors [18].
2.5 Statistic
The statistical analyses were conducted using TwoSampleMR (version 1.0.3) and MRPRESSO (version 6.1) in R (version 0.4.25). The global level tests were considered significant at a two-sided P value of 0.05. For region-level analyses, we corrected for multiple testing using the Bonferroni method for 408 MR estimates, resulting in a significance threshold of P < 1.22 × 10-4. We also considered results with a P value < 0.05 as nominally significant.
2.6 Ethical statement
This study did not involve human or animal experiments, but only used existing genetic data. We used genetic data from UK Biobank and ENIGMA Consortium, which were collected and used with the informed consent of the participants and in accordance with the data sharing and confidentiality principles. UK Biobank has obtained Research Tissue Bank (RTB) approval from the North West Multi-centre Research Ethics Committee (MREC), so we did not need to apply for separate ethical review. We cited the sources of these data and followed the usage rules of UK Biobank and ENIGMA Consortium.