According to Neale Lab genome-wide association meta-analysis (GWAS round 2),¹ the GWAS summary statistics for tea intake (phenotype code: 1488_raw) were based on 349,376 samples of European ancestry from UK Biobank. This GWAS adjusted for age, sex, age 2, sex × age 2, sex × age, and the top twenty principal components. Habitual tea intake was retrieved based on the item: “How many cups of tea do you drink each day?” All instrumental variables were associated with the exposure (e.g., tea consumption) at a genome-wide significance level (P < 5 × 10^–8) with linkage disequilibrium (LD) r² < 0.01 at a 10,000 kb window, which confirmed the independence for the selected genetic variants. In this MR study, 22 independent SNPs with moderate LD were selected as genetic instruments for habitual tea consumption. 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 knee OA, hip OA, RA, and OP. The Summary data was available from a wide-ranging meta-analysis of transethnic GWAS. Over 100,000 participants of European and Asian descent participated in the survey. To reduce population stratification bias, we only analyzed data for European, including knee OA (17,885 controls and 4462 cases), hip OA (50,898 controls and 12,625 cases), and RA (43,923 controls and 14,361 cases) from the large GWAS. The secondary outcomes were specifically diagnosed OP. Data of OP was taken from the FinnGen Data Freeze 2.² On January 14, 2020, FinnGen Data Freeze 2 was released to the public. It includes 96,499 individuals, 16 M variants, and 1,122 endpoints. The dataset contains 4,477 OP cases (93083 controls and 1175 cases). Study protocols related to these data have been released and described in the previous studies (19–21). All outcome summary GWAS data in this study were derived from ‘‘ieu open gwas project.³ “The GWAS IDs for the outcomes are as follows: ebi-a-GCST005813 for Knee OA, ebi-a-GCST005810 for hip OA, ieu-a-832 for RA, finn-b-M13_OSTEOPOROSIS for OP. The relevant ethics committees approved all studies that contributed data to these analyses, and all participants provided written informed consent.

This was a two-sample MR analysis. MR, using genetic variants as instrumental variables, is a genetic epidemiologic method that can avoid potential confounding factors and reverse causation. It is assumed that the instrumental variables are valid if it meet all of the following criteria: (1) strongly associated with tea consumption, (2) not associated with confounding factors, and (3) risk association with knee OA, hip OA, RA, and OP were only detected via tea consumption (22, 23) (Figure 1). We tested for weak instrumental variables using F statistics. The F statistic equals [(n-k-1)/k) × R^2/(1-R^2)], where R² represents the variance in tea consumption explained by the genetic instrument, K represents the number of genetic variations, and N represents the sample size. The R² was calculated as follows: 2 × beta² × EAF × (1-EAF)/2 × beta² × EAF × (1-EAF) + se2 × 2 × N × EAF(1-EAF). EAF represents effect allele frequency (24). When F > 10, it is unlikely that a weak instrument will bias the MR estimates. Nevertheless, there may be pleiotropic genetic variants introduced, thus invalidating the IV.

The Cochran Q-test was applied to check SNPs’ statistical heterogeneity using MR-Egger estimates, with P < 0.05 deemed significantly heterogeneous. Therefore, the IVW method based on random effects was adopted. Several approaches were used to test whether the second and third MR assumptions were violated, including the MR-Egger regression (25), the weighted median (26) and mode (27) methods, an outlier test (MR-PRESSO) (28). MR-Egger is a weighted regression approach that introduces an intercept to accommodate pleiotropy. There was horizontal pleiotropy when the intercept term was away from zero (25). Using this approach, unbiased estimates are achieved in the presence of pleiotropic instruments with the assumption that the magnitude of pleiotropic effects cannot be predicted by the size of the instrumental variables—SNPs related to tea consumption (25).

In addition, funnel plots were used to assess pleiotropy. If the funnel plot is symmetric, there is no potential directional pleiotropy (29). We performed a “leave-one-out” analysis by excluding each SNP in turn. Therefore, individual SNP could be assessed for their influence on causal associations. The PhenoScanner database (version 2.0) was examined for each SNP (30) to determine whether existing risk factors have previously been associated with any significant associations (P < 5 × 10^–8): BMI (31, 32), previous knee injury (32), gout (33), coffee intake (12) for knee OA; Hip dysplasia (34), Steroids use (35), coffee intake (11) for hip OA; hypertension (33), diabetes, smoking (36) for RA; Vitamin D intake, Steroids use, smoking for OP (37). To rule out potential pleiotropic effects, we evaluated the effect of removing these SNPs from the MR estimations.

All analyses were conducted using the “TwoSampleMR,” “MRPRESSO” R packages in RStudio version 3.6.3. The Bonferroni adjustment was used to correct for multiple comparisons (P-value: 0.05/4 outcomes = 0.0125). We computed two-side P-values, with P < 0.0125 regarded as statistically significant.

All 22 SNPs together as the instrument together have an F-statistic more than 10 meeting the first assumption for MR study (Table 1). The PhenoScanner database showed one SNPs (rs10741694) associated with RA confounding factors (hypertension, diabetes, smoking), four other SNPs (rs11022751; rs1481012; rs2472297; rs4410790) related to the confounding factors (BMI; previous knee injury; coffee intake) of the knee OA, two other SNPs (rs2472297; rs4410790) related to the confounding factors (coffee intake) of the hip OA. In order to satisfy assumption 2, these SNPs will be eliminated when the corresponding analysis is performed. We used a range of sensitivity analyses to satisfy assumption 3.

We found four SNPs (rs11022751; rs1481012; rs2472297; rs4410790) related to the confounding factors (BMI; previous knee injury; coffee intake) of the knee OA through the PhenoScanner database (Supplementary Table 1). Weeding out the four SNPs, There was no evidence for the causal association between genetically predicted tea intake and knee OA. The OR was 1.11 (95% CI: 0.79–1.55), using IVW with random effect in MR (Figure 2). The weighted Median methods (OR: 1.18, 95% CI: 0.81,1.74) did not support a causal relationship between tea consumption and the risk of knee OA. The estimates for each SNP and scatter plot for two analyses on knee OA were shown in Supplementary Figures 1, 2.

In MR-Egger regression, the intercept term indicated that horizontal pleiotropy was unlikely to affect the result (Egger intercept: 0.01, P-value: 0.77). Additionally, As the funnel plot was visually symmetrical, this indicated that pleiotropy had not been established. SNPs did show slight heterogeneity according to Cochran’s Q-test [Q-value (df) = 25.08 (13), P = 0.03 for MR Egger method; Q-value (df) = 25.25 (15), P = 0.05 for IVW method]. “Leave-one-out” analysis demonstrated that overall risk estimates were not significantly affected by any single SNP (Supplementary Figure 3).

After excluding 2 SNPs (rs2472297; rs4410790) associated with coffee consumption, Based on IVW analyses, no causal relationship was found between tea consumption and hip OA (OR: 1.20, 95% CI: 0.84–1.71) (Figure 3). In MR-Egger regression, the intercept term indicated no horizontal pleiotropy affecting the result (Egger intercept: −0.07, P-value: 0.08). As for genetics instruments, Cochran’s Q statistic suggested no heterogeneity: Q-value (df) = 16.31 (16), P = 0.43 for MR Egger method; Q-value (df) = 19.82 (17), P = 0.28 for IVW method. The forest and scatter plots can be seen in Supplementary Figures 4, 5 and the leave-one-out plot in Supplementary Figure 6.

The PhenoScanner database showed 2 SNPs (rs10741694; rs397074) associated with RA confounding factors. Excluding the two SNPs. Based on genetic instruments using IVW with random effect analyses, there was no causal link between genetically predicted tea consumption and RA (OR: 1.24 95% CI: 0.81–1.91) (Figure 4). The forest plots can be seen in Supplementary Figure 7. The scatter plot for the two analyses was attached as Supplementary Figure 8 and the leave-one-out plot in Supplementary Figure 9.

As for genetics instruments, Cochran’s Q statistic suggested slight heterogeneity: Q-value (df) = 20.72 (12), P = 0.05 for MR Egger method. Additionally, the intercept of the MR-Egger regression showed no signs of pleiotropy (Egger intercept: 0.05, P-value: 0.09).

The risk of OP did not increase with higher levels of tea consumption due to genetic predisposition. The causal ORs for each cup increased in tea consumption were 1.11 (95% CI: 0.89,1.39) in the IVW estimates. At the same time, MR-Egger and Weighted median methods also confirm no causal relationship between tea consumption and OP, and the P-value were 0.24 and 0.11, respectively (Figure 5). No heterogeneity found in Cochran’s Q statistic [Q-value (df) = 17.47 (15), P = 0.29 for MR Egger method; Q-value (df) = 18.38 (16), P = 0.30 for IVW method]. Detailed forest and scatter plots were shown in Supplementary Figures 10, 11. The leave-one-out analysis shows that none of the single SNP substantially affects the overall risk estimation (Supplementary Figure 12).