The European population’s genome-wide association meta-analysis (GWAS) dataset identified SNPs associated with coffee consumption. As for primary genetic instruments, 15 significant SNPs associated with coffee intake were obtained from a GWAS meta-analysis that comprised over 370,000 participants of European ancestry (25). This GWAS adjusted for age, sex, BMI, total energy, the proportion of 24 hours recalls self-reported as capturing “typical intake,” and top twenty principal components (25). Habitual coffee intake was retrieved based on the item: “How many cups of coffee do you drink each day (include decaffeinated coffee)? “ Meanwhile, the effects of SNPs were interpreted as 1% change of coffee consumption per effect allele (25). As for secondary genetic instruments, 10 significant SNPs associated with the cups of coffee consumed per day were identified from another GWAS, including 91,462 participants released by the Coffee and Caffeine Genetics Consortium (CCGC) (26). All instrumental variables were associated with the exposure (e.g., coffee consumption) at a genome-wide significance level (P<5×10^-8) with linkage disequilibrium (LD) r² <0.001 at a 10,000 kb window, which confirmed the independence for the selected genetic variants. In this MR study, 11 and 8 independent SNPs with moderate LD were selected as genetic instruments for habitual coffee consumption after excluding four SNPs (rs117692895, rs4719479, rs12699844, and rs73073176 in chromosome 7) and two SNPs (rs6968554 in chromosome 7; rs247083 in chromosome 15) in primary genetic instruments and secondary genetic instruments, respectively. Detailed information on the relationship between the selected SNPs and exposures is shown in Table 1 .

The primary outcome in this study was the clinically diagnosed OA. The summary-level data for OA were obtained from the large GWAS conducted by the UK Biobank, enrolling 50,508 participants of European ancestry (10,083 with cases and 40,425 controls). Moreover, the secondary outcomes were specific diagnosed knee OA (4462 cases and 17,885 controls) and hip OA (12,625 cases and 50,898 controls), and self-reported OA (12,658 cases and 50,898 controls). Study protocols related to these data have been released and described in the previous studies (27–29). All summary data can be obtained from the UK Medical Research Council Integrative Epidemiology Unit Open GWAS Project database (httep://gwas.mrcieu.au.uk). The relevant ethics committees approved all studies that contributed data to these analyses, and all participants provided written informed consent.

The conventional MR method was applied in this study (https://mrcieu.github.io/TwoSampleMR/articles/perform_mr.html#mr-methods-1) (30). The random-effects inverse-variance weighted (IVW) model was conducted to examine the causal association, and this approach was considered as the main analysis because of the potential observed heterogeneity (31, 32). The IVW method combines individual MR effects across SNPs to derive an overall weighted effects of the potential causal association. Furthermore, the forest plots were showed to visualize the MR-derived odds ratio (OR) of overall, knee, hip, and self-reported OA risk for 1% or one cup per day increase in genetically predicted coffee intake. In the sensitivity analyses, “leave-one-out” method was performed to estimate that the causal association was reliant on any single SNP. In addition, Steiger-MR was used to test whether the SNPs explained significantly more variance in exposure than outcome (the opposite may indicate reverse causation).

The IVW method assumes that all genetic variants should satisfy three assumptions for the instrumental variables. 1) strongly associated with coffee intake, 2) not associated with confounders of the association between coffee intake and OA, and 3) the association with OA risk was only found via coffee consumption (33, 34). Thus, the F statistics were used to test for weak instrumental variables. F statistics=((n-k-1)/k)(R^2/(1-R^2)), where R² is the variance in coffee consumption explained by the genetic instrument, k is the number of genetic variants, and n is the sample size. The strengths of the primary (11 SNPs) and Secondary (8 SNPs) genetic instruments used for analyses were 70.59 and 28.46, respectively. F > 10 was proven to use strong genetic instruments in MR study. The Cochran’s Q test was used to quantify the heterogeneity in effect sizes between the genetic instruments (35), which may indicate horizontal pleiotropy that could violate the third MR assumption. Potential violation of the second and third MR assumptions was tested using several approaches such as the MR-Egger regression (36), the weighted median (37) and mode (38) methods, and the MR pleiotropy residual sum and outlier test (MR-PRESSO) (39). The association between the selected SNPs and exposures was validated in the PhenoScanner database (http://www.phenoscanner.medschl.cam.ac.uk/) ( Supplementary Table 1 ). SNP associated with traits other than coffee consumption at the GWAS significance level was documented.

The MR-Egger approach is an adaption for Egger regression that allows for directional pleiotropy by introducing an intercept in the weighted regression model. Horizontal pleiotropy was indicated when values were away from zero for the intercept term (36). Based on this approach, unbiased estimates are performed in the presence of pleiotropic instruments assuming that the magnitude of pleiotropic effects is independent of the size of the instrumental variables—SNPs associated with coffee intake (36).

Weighted median method orders the estimates in MR using each instrument weighted for the inverse of their variance, and the median result is selected and shows the single MR estimate with confidence intervals based on bootstrapping technique (37). The weighted median requires and assumes that at least half of the instruments are valid (40).

The mode-based causal estimate consistently estimates the true causal effect when the most instruments with consistent MR estimates are valid (38).

The “twosampleMR” package (version 0.5.5) and R software version 3.6.1 were used for all statistical analyses. Bonferroni correction (P = 0.05/2 exposures/4 outcomes = 0.0625) was implemented for multiple comparisons. Two-sided P values were computed, with P < 0.00625 regarded as statistically significant.

The primary and secondary genetic instruments using IVW with random effect analyses provide strong evidence for the causal association between genetically predicted coffee consumption and overall OA (OR:1.009, 95%CI:1.003-1.016; OR:1.270, 95%CI:1.099-1.469) ( Figures 1A, B ). The estimates for each SNP on overall OA were shown in Supplementary Figures 1 , 2 . The scatter plot for two analyses was attached as Supplementary Figures 3 , 4 .

As for primary genetics instruments, Cochran’s Q statistic suggested a slight sign of heterogeneity: Q value (df)=18.34 (9), P=0.03 for MR Egger method; Q value (df)= 19.56 (10), P=0.03 for IVW method. However, the significant effect was support for Weight Median methods that can be given priority when the model with heterogeneity but no pleiotropy (OR:1.008, 95%CI:1.002-1.014). In addition, there was no indication of pleiotropy when the intercept was derived from the MR-Egger regression (Egger intercept:0.01, P value:0.46). Furthermore, no outliers were detected with potential pleiotropy using the MR-PRESSO method. Moreover, the P-value for the MR-PRESSO Global test was 0.073, which also suggested no significant sign of heterogeneity. The leave-one-out analysis results suggested no influential SNPs from the causal link between coffee intake and OA in the replication analyses ( Supplementary Figure 5 ). Steiger filtering identified 2 instruments explaining more variation in total OA than in caffe consumption. Removing those slightly attenuated the identified effect estimate, which was still suggestive of a significant effect of genetic liability to higher cups of coffee intake per day on risk of total OA in the primary analyses ( Supplementary Table 2 ).

In the secondary analyses, Cochran’s Q statistic suggested no sign of heterogeneity: Q value (df)=4.50 (6), P=0.61 for MR-Egger method; Q value (df)= 7.39 (7), P=0.38 for IVW method. The associations were consistent with using the weighted median method (OR = 1.228, 95% CI: 1.025–1.472; P = 0.026). In addition, there was no indication of pleiotropy when the intercept was derived from the MR-Egger regression (Egger intercept:0.02, P value:0.14). Furthermore, no outliers were detected with potential pleiotropy using the MR-PRESSO method. The leave-one-out analysis results suggested no influential SNPs from the causal link between coffee intake and OA in the replication analyses ( Supplementary Figure 6 ). Steiger filtering suggested that the direction of the effect was correct for all the coffee consumption instruments in the secondary analyses.

As for primary genetic instruments, IVW with random effect suggested marginally significant evidence for the causal association between coffee consumption and OA (OR:1.01,95% CI:0.999-1.025). However, strong heterogeneity was found in both MR-Egger (Q value:32.81, P value<0.0001) and IVW methods (Q value: 32.96, P value<0.0001), respectively. The hypothesis of the pleiotropy test was satisfied in the fitting model. By MR PRESSO method, one SNP (rs34060476) was detected as an outlier, and the Global test was 0.012. Based on leave-one-out methods, we found 4 SNPs (rs4410790, rs597045, rs34060476, and rs1956218) that were potential sources of heterogeneity ( Supplementary Figure 11 ). After excluding these SNPs (rs4410790, rs597045, rs34060476 and rs1956218), results from the Weight Median and IVW methods suggested a strong causal association (OR:1.019,95%CI:1.008-1.030; OR:1.023, 95%CI:1.009-1.038) ( Figure 2A ). The forest and scatter plots were performed in Supplementary Figures 7 , 9 , respectively.

In the secondary analyses, the causal relationship was supported for the IVW method (OR:1.410, 95% CI:1.085-1.832), and The effect was slightly attenuated with the Weight Median method, yielding a null causal estimation between coffee consumption and OA (OR:1.251, 95% CI: 0.923-1.696) ( Figure 2B ). The forest and scatter plots were performed in Supplementary Figures 8 , 10 , respectively. There was no heterogeneity and outlier in the model (P >0.05). The leave-one-out analyses revealed that the result was essentially consistent after removing SNP in turn ( Supplementary Figure 12 ).

Steiger filtering suggested that 3 and 1 SNPs associated with knee OA explained more variation in coffee consumption and these were removed in primary and secondary dataset, respectively. The exclusion of these SNPs, despite attenuating the two analyses effect estimate, were suggestive of an effect of genetic liability to higher coffee consumption on knee OA ( Supplementary Table 2 ).

The causal association between coffee consumption and hip OA was not found based on IVW analyses, either primary or secondary genetic instruments (OR:1.012, 95% CI: 0.999-1.024; OR: 1.199, 95% CI: 0.896-1.606) ( Figures 3A, B ). The hypothesis of both heterogeneity and pleiotropy was testified (P >0.05). The forest and scatter plots were performed in Supplementary Figures 13 – 16 , and the plot for leave-one-out analyses was shown in Supplementary Figures 17 , 18 . Steiger filtering which explain more variation in the outcome across two analyses, the effect estimate attenuated further, providing limited evidence of a direct effect of genetic liability to higher coffee consumption on the risk of hip OA ( Supplementary Table 2 ).

A significant causal relationship was found between genetically predicted coffee consumption and self-reported OA in primary and secondary genetic instruments ( Figures 4A, B ). In primary analyses, the OR for OA was 1.007 (95% CI:1.003-1.011) using the IVW method and 1.007 (95% CI:1.002-1.012) using the Weight Median method; In secondary analyses, the OR for OA was 1.249 (95%CI:1.101-1.417) using the IVW method and 1.223 (95% CI:1.044-1.433) using weighted median. There was no significant heterogeneity and pleiotropy (P >0.05). None outlier was confirmed in MR PRESSO analyses. Detailed forest and scatter plots were shown in Supplementary Figures 19 – 22 . The plot for leave-one-out analyses was shown in Supplementary Figures 23 , 24 .

Steiger filtering suggested that the effect direction was correct for all the coffee consumption SNPs, and no SNP was removal ( Supplementary Table 2 ).