Hypothyroidism and rheumatoid arthritis: a two-sample Mendelian randomization study

Background: Meta-analysis of genome-wide association studies (GWAS) data showed that the relationship between hypothyroidism and rheumatoid arthritis (RA) risk remains under debate. This study is conducted to test the causal relationship of hypothyroidism and RA.

Methods: A two-sample Mendelian randomization (TSMR) analysis was employed to estimate the causality of hypothyroidism and rheumatoid arthritis in European ancestry and Asian ancestry. Integrating the effects generated by TSMR, functional annotations and noncoding variant prediction framework were applied to analyze and interpret the functional instrument variants (IVs).

Results: The results of the inverse variance weighted method showed a strong significant causal relationship between hypothyroidism and risk of RA in European ancestry [odds ratio (OR) = 1.96; 95% confidence interval (CI) 1.49, 2.58; p < 0.001]. The outcomes of MR-Egger, weighted median, weighted mode, and simple mode also showed that hypothyroidism was significantly associated with increased risk of RA in European ancestry. The MR-PRESSO method also showed significant results [Outlier-corrected Causal Estimate = 0.70; standard error (SE) = 0.06; p < 0.001]. An independent dataset and an Asian ancestry dataset were applied to estimate and obtain the coincident results. Furthermore, we integrated the effect of variants in TSMR analysis, functional annotations, and prediction methods to pinpoint the single-nucleotide polymorphism (SNP) rs4409785 as one of the causal variants, which suggested that this variant could impact the binding of CTCF-cohesin and play a vital role in immune cells.

Conclusion: In this study, we prove that hypothyroidism is significantly causally associated with increased RA risk, which has not been shown in previous studies. Furthermore, we pinpoint the potential causal variants in RA.

Background

Rheumatoid arthritis (RA) is one of the most common autoimmune and chronic inflammatory diseases. Almost 0.5% to 1% of the world population are affected by RA, which leads to bone damage, disability, and premature mortality (1, 2). Both genetic and environmental triggers may contribute to the onset and progress of RA (3). Thyroid-related diseases including hypothyroidism and hyperthyroidism can be divided into subclinical stage and overt, and autoimmune thyroid disease can also occur in patients with other systemic autoimmune disorders such as RA (4–7). The association study of Graves’ disease patients with other autoimmune diseases showed that RA is one of the most frequently observed diseases (8). Graves’ disease is one of the causes of secondary hypophysitis and an autoimmune pattern such as lacking T regulator cells in hypophysitis is also seen with RA and other autoimmune conditions (9). A Mendelian randomization (MR) study found the causal relationship between Graves’ disease and RA (10).

Previous studies showed that the rate of thyroid dysfunction diseases, especially hypothyroidism, was significantly increased in RA patients, compared with the prevalence in controls (11). McCoy et al. reported that the different presence of hypothyroid disease between cohorts of patients with or without RA was not significant at the time of RA diagnosis in the earlier stage (12). Recently, the meta-analysis also showed that the risk of developing thyroid dysfunction particularly hypothyroidism was increased in RA patients (13). However, in many cases, the earlier diagnosed RA or hypothyroidism does not mean that it begins earlier; the relationship between hypothyroidism and RA remains elusive.

MR uses the genetic markers as instruments of the exposure to assess the causal relationship between exposure and outcome (14). To decipher the causal relationship between hypothyroidism and RA, we adopted a two-sample Mendelian randomization (TSMR) that employed “MR-Egger”, “Weighted median”, “Inverse variance weighted”, “Simple mode”, and “Weighted mode” methods (Figure 1). Genome-wide association studies (GWAS) have identified a number of variants associated with diseases and traits that fall into noncoding regions (15). Though TSMR could infer the causal relationship between two diseases or traits, the causal variants in the instrument variants (IVs) are still hard to pinpoint. Furthermore, we also pinpointed the potential causal variants of IVs by integrating genomic and epigenetic annotations. Therefore, our study integrated MR and functional genomics data to investigate causal relationships between hypothyroidism and RA, and to pinpoint the potential causal variants to provide some clinical helpful insights for the diagnosis and therapy of RA patients.

FIGURE 1

Figure 1 Schematic of an MR analysis. The causal estimation of hypothyroidism and RA. The first key assumption is that the genetic variant selected as the IVs should be significantly associated with exposure. The second key assumption is that the IVs should not be associated with any confounding factors. The third key assumption is that the IVs should affect the outcome merely through exposure but not any alternative pathways. SNP, single-nucleotide polymorphism.

Methods

Five independent summary datasets were used to conduct MR analysis separately. The first summary data of RA were collected from the FinnGen R5 release (16), which included 6,236 cases and 147,221 controls of European ancestry in this dataset, and the IEU open GWAS project (ID: “finn-b-M13_RHEUMA”); the controls of the “M13_RHEUMA” were excluded from the listed items in M13_ARTHROPATHIES. The second summary data of RA were collected from IEU open GWAS project (ID: “ebi-a-GCST000679”), which included 5,539 cases and 20,169 controls of European ancestry. The RA cases in the second summary data of RA met the 1987 American College of Rheumatology (ACR) criteria for diagnosis of RA and were autoantibody-positive individuals; controls were matched to BRASS RA cases using principal components analysis from GWAS data from three separate studies (17). The hypothyroidism summary data were collected from the FinnGen R5 release (16) where hypothyroidism was defined as “hypothyroidism, drug reimbursement” (ID: “finn-b-HYPOTHY_REIMB”), meaning the phenotype of drug reimbursement of hypothyroidism patients, which included 7,183 cases, and 59,893 controls of European ancestry. The second hypothyroidism dataset defined as self-reported hypothyroidism/myxoedema was collected from IEU open GWAS project (ID: “ukb-a-77”), which included 16,373 cases and 320,783 controls in UK Biobank with filed Data-Field 20002 (https://biobank.ctsu.ox.ac.uk/crystal/field.cgi?id=20002). In addition, to analyze the treatment drug levothyroxine sodium relationship, we also used the hypothyroidism dataset with treatment/medication code: levothyroxine sodium summary data were collected from IEU open GWAS project (ID: “ukb-a-190”), which included 13,717 cases and 323,442 controls in UK Biobank with filed Data-Field 20003 (https://biobank.ctsu.ox.ac.uk/crystal/field.cgi?id=20003). The ancestry of the above datasets is European; findings were verified on the Asian ancestry using another independent dataset (18). GWAS summary statistics on RA with Asian ancestry were downloaded from the GWAS Catalog with 5,348 cases and 173,268 controls (http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST90018001-GCST90019000/GCST90018690/harmonised/34594039-GCST90018690-EFO_0000685.h.tsv.gz). GWAS summary statistics on hypothyroidism with Asian ancestry was downloaded from the GWAS Catalog with 1,114 cases and 172,656 controls (http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST90018001-GCST90019000/GCST90018642/harmonised/34594039-GCST90018642-EFO_0004705.h.tsv.gz).

To further verify the conclusion obtained from the above datasets, we adopted three recently larger datasets. We used the summary data of RA (finngen_R8_M13_RHEUMA, https://storage.googleapis.com/finngen-public-data-r8/summary_stats/finngen_R8_M13_RHEUMA.gz), which included 11,178 cases and 221,323 controls, and hypothyroidism summary data (finngen_R8_HYPOTHY_REIMB, https://storage.googleapis.com/finngen-public-data-r8/summary_stats/finngen_R8_HYPOTHY_REIMB.gz), which included 10,943 cases and 81,899 controls, collected from the FinnGen R8 release (19). In the second and third MR experiments, we adopted the summary data of RA in Okada et al., which included 14,361 cases and 42,923 controls of European ancestry (http://plaza.umin.ac.jp/~yokada/datasource/files/GWASMetaResults/RA_GWASmeta_European_v2.txt.gz); the cases in this dataset fulfilled the 1987 criteria of the ACR for RA diagnosis or were diagnosed as RA by professional rheumatologists (20). The cases in our study from FinnGen and UK Biobank were the registry data with clinically meaningful disease endpoints by the clinical expert groups.

Instrument variants selection

Variants related to hypothyroidism and RA risk at p-value < 5×10^-8 were identified from GWAS summary data of European ancestry. To obtain the independent IVs, we used the 1,000 Genome linkage disequilibrium (LD) in European ancestry with r² larger than 0.01 in the 5,000-kb region as the threshold to apply data clumping. For the outcome variants, the minimum r² is 0.8 and the palindromic single-nucleotide polymorphisms (SNPs) are allowed. We also used the minor allele frequency (MAF) threshold 0.3 to infer palindromic SNPs. F statistic was employed to evaluate the impact of IVs (21), the formula is R²×(N − 2)/(1 − R²), where R² is the cumulative explained variance of selected SNPs in exposure that used (2×EAF×(1 − EAF)×beta²)/[(2×EAF×(1 − EAF) ×beta²) + (2×EAF×(1 − EAF)×N×SE(beta)²)], where N is the sample size of research, EAF is the effect allele frequency, beta is the estimated genetic effect, and SE (beta) is the standard error of the beta. The strong IVs were chosen when F > 10.

MR analysis

The three key assumptions (22) of IVs that are vital to the valid IVs for causal inference in MR analysis should be satisfied (1): the selected IVs should be directly associated with exposure; (2) the selected IVs should not be associated with confounders; (3) the selected IVs must have no effects on the outcome other than through the exposure (no horizontal pleiotropy exists).

Five complementary methods of TSMR were conducted on GWAS summary data to estimate the causal relationship that included MR-Egger, weighted median, inverse variance weighted (IVW), simple mode, and weighted mode. The IVW method, which uses the ratio estimate of each variant in a fixed effect meta‐analysis to estimate the causal effect (23), is strictly based on the three assumptions. Under the INSIDE assumption (instrument strength independent of the direct effects) (22), the MR-Egger method conducts a weighted linear regression rather than setting the intercept to be zero in IVW, and the intercept in MR-Egger can be used to estimate the horizontal average pleiotropic effect (24). If up to 50% of genetic variants are invalid, the total weight of the instrument comes from the median of the weighted ratio estimates of valid variants, then this method is called the weighted median method (25). We also incorporated mode-based methods that include simple mode and weighted mode. These methods estimate the causal effect of individual SNPs to form clusters (26). The simple mode takes the largest cluster of SNPs’ causal estimation, but the weighted mode assigns the weights to each SNP (27). If the beta values in summary data of exposure and outcome have significantly different distributions, the correct factor will be used.

The TwoSampleMR package (version 0.5.6) in R (version 4.0.2) was used to carry out the five MR methods.

Pleiotropy, heterogeneity, and sensitivity evaluation

We required the directions of all five methods to be consistent, and p-values less than 0.05 indicate significant results. The MR-Egger method was performed to return intercept values to test the horizontal pleiotropy. If directional pleiotropy was detected with the p-value < 0.05, we employed MR-PRESSO to evaluate the horizontal pleiotropy in the generated MR model and removed the outlier and then estimated the causal effect again. Cochran’s Q statistic of IVW was applied to check the effect of heterogeneity. Leave-one-out and single SNP analyses were also performed to identify if a single SNP is driving the causal estimates.

Functional instrument variant prioritization

We merged all IVs and used the b values in single SNPs analysis as the b score of each variant. We extracted variant call format (VCF) files of the merged IVs and adopted the Sei framework (28) (https://github.com/FunctionLab/sei-framework), which integrated 21,907 chromatin profiles to predict functional noncoding variants. We used the max absolute difference in 40 sequence classes of each variant as the Sei score. The top 10 variants with high b and Sei scores were the selected variants to reveal the potentially affected features in a more detail manner. The framework of functional analysis with details is described in Supplementary Figure 1. To visualize the regulatory features of selected variants in immune cells, we downloaded ReMap (https://remap.univ-amu.fr/storage/remap2022/hg38/MACS2/remap2022_all_macs2_hg38_v1_0.bed.gz) and CCCTC-binding factor (CTCF), estrogen receptor alpha gene (ESR1), androgen receptor (AR), and recombinant forkhead box protein A1 (FOXA1) chromatin immunoprecipitation sequencing (ChIP-seq) data (Supplementary Table 1). GWAS Catalog track was downloaded from the UCSC genome browser (https://genome.ucsc.edu/cgi-bin/hgTables?hgsid=1509267435_AXk036LRJbjcQTT5fpxOODCrVW6Y). PolyPhen-2 (29) was used to predict the function of variants falling into coding regions.

Results

RA and hypothyroidism, drug reimbursement

The results showed a significant causal relationship between hypothyroidism, drug reimbursements of hypothyroidism and risk of RA in the European ancestry through the IVW method (OR = 1.96; 95% CI 1.49, 2.58; p < 0.001). The result of the MR-Egger method (OR = 3.08; 95% CI 1.57, 6.04; p = 0.006) also showed a significant result.

The results of weighted median, simple mode, and weighted mode also showed a significant relationship between hypothyroidism, drug reimbursements of hypothyroidism and risk of RA (Figure 2A). MR-Egger (Cochran’s Q = 331.25; p < 0.001) and IVW methods (Cochran’s Q = 379.72; p < 0.001) were used to analyze the heterogeneity of outcomes, and there was heterogeneity among the selected IVs. Then, the IVW (multiplicative random effects) method was adopted to estimate the causal relationship, and the consequence was consistent (OR = 1.96; 95% CI 1.49, 2.58; p < 0.001). Though the heterogeneity existed among the selected IVs, it does not affect the results of IVW, and the relationship between hypothyroidism, drug reimbursements of hypothyroidism and risk of RA is still significant. MR-Egger was used to evaluate horizontal pleiotropy among the selected IVs; the results (β intercept = −0.09; SE = 0.07; p = 0.174) showed that no pleiotropy would affect the results (Figure 3). The MR-PRESSO method also showed significant results (Outlier-corrected Causal Estimate = 0.70; SE = 0.06; p < 0.001).

FIGURE 2

Figure 2 Summary view of the MR results. Summary of MR analysis results derived from the inverse variance weighted, MR-Egger, simple mode, weighted mode, and weighted mean methods. RA was used as the outcome. The hypothyroidism and drug reimbursements (A), non-cancer illness code self-reported: hypothyroidism/myxoedema (B), and levothyroxine sodium (C) as the exposures and significant associations were detected for these traits and RA risk.

FIGURE 3

Figure 3 SNP effect evaluation in TSMR between hypothyroidism, drug reimbursement, and RA. (A) MR test scatter plot of five methods. The x-axis is the SNP effect on hypothyroidism, drug reimbursement. The y-axis is the SNP effect on rheumatoid arthritis. (B) MR funnel plot of IVW and MR-Egger methods. (C) Forest plot of MR sensitivity analysis. All MR-Egger and IVM methods showed that MR effect sizes that are larger than 0 mean that hypothyroidism, drug reimbursement had a causal effect on RA. (D) Forest plot of MR leave-one-out sensitivity analysis. MR, Mendelian randomization; SNP, single-nucleotide polymorphism.

RA and self-reported hypothyroidism/myxoedema

To test the conclusion from hypothyroidism, drug reimbursements of hypothyroidism and risk of RA, we also applied TSMR on self-reported hypothyroidism/myxoedema and RA of European ancestry. Five methods all showed significant relationships between self-reported hypothyroidism/myxoedema and RA of European ancestry. The IVW result showed that self-reported hypothyroidism/myxoedema significantly increased the risk of RA (OR = 1.84; 95% CI 1.62, 2.10; p < 0.001). MR-Egger also showed a significant relationship between the exposure and outcome (OR = 1.64; 95% CI 1.12, 2.39; p = 0.013) (Figure 2B). The heterogeneity results evaluated by MR-Egger (Cochran’s Q = 84.10; p = 0.005) and IVW (Cochran’s Q = 84.10; p = 0.007) showed that there was heterogeneity among the selected IVs. MR-Egger was used to evaluate horizontal pleiotropy among the selected IVs, and the results (β intercept = 5×10^-4; SE = 0.01; p = 0.970) showed that no pleiotropy would affect the results. In addition, the results of single SNP analysis and leave-one-out sensitivity analysis indicated that no single SNP significantly leads to the causal results (Supplementary Figures 2C, D).

RA and levothyroxine sodium

Levothyroxine sodium is a medicine used to treat hypothyroidism (30). We also found in the trait of people who took levothyroxine sodium, which had significant causal effects on RA risk. Forty-two SNPs were selected as the IVs and we used MR-Egger to analyze the horizontal pleiotropy; we further found that there was pleiotropy among IVs (β intercept = −0.07; SE = 0.04; p = 0.060). We then adopted MR-PRESSO to analyze and filter the nine outlier SNPs, and used MR-Egger again to test the horizontal pleiotropy and eliminate the horizontal pleiotropy (β intercept = 0.02; SE = 0.02; p = 0.245). The causal relationship was evaluated with the selected IVs; the IVW method (OR = 1.47; 95% CI 1.36, 1.60; p < 0.001) and the MR-Egger method (OR = 1.40; 95% CI 1.08, 1.81; p = 0.001) also showed significant consequence (Figure 2C). MR-Egger (Cochran’s Q = 31.43; p < 0.001) and IVW methods (Cochran’s Q = 33.14; p < 0.001) were used to analyze the heterogeneity of outcomes, and there was heterogeneity among the selected IVs. The IVW (multiplicative random effects) method (OR = 1.49; 95% CI 1.37, 1.62; p < 0.001) also showed significant results. Moreover, the single SNP analysis and leave-one-out analysis showed that no single SNP affects the causal result significantly (Supplementary Figures 3C, D).

Functional analysis of IVs

We deciphered the relationship between hypothyroidism and RA, but MR results could not provide the potential functional mechanisms of these diseases. Previous studies showed that almost 90% of variants in GWAS fall into noncoding regions. Based on the third assumption of MR analysis, the selected IVs have to affect the outcome through exposure. Selected IVs of three experiments were merged into a 102-variant dataset; Chen et al. developed the Sei framework (28) to discover the regulatory code of human genetics data based on sequence information of traits and diseases. We adopted Sei framework to predict the 92 variants and regarded the max absolute fold change in the 40 sequence classes (Supplementary Table 2). To obtain the MR estimates using each of the selected IVs singly, we employed the Wald ratio to perform the single SNP MR analysis. We used the b values in the single SNP analysis as b score. Integrating b score and Sei score, we filtered the top 7 variants with absolute Sei score > 1 and b core > 0.5 (Figure 4). SNP rs3850765 has the largest b score and a larger fold change on the TF3 sequence class such as FOXA1, AR, and ESR1 (Supplementary Figure 4). The UCSC genome browser annotation suggested that this variant fell into the regulatory regions (Supplementary Figure 5). Another SNP rs4409785 has the largest fold change on CTCF and CTCF-cohesin and a larger b score in self-reported hypothyroidism/myxoedema and RA, and levothyroxine sodium and RA. This variant fell into the regulatory enrichment regions and the peak of CTCF in GM12878. Compared with the CD8-positive alpha-beta T cell, the variant locus is more activated in the activated CD8-positive alpha-beta T cell in CTCF-ChIP data. The chromatin interaction analysis by paired-end tag sequencing (ChIA-PET) of CTCF in GM12878 also showed this variant located at the anchor region (Figure 5). The UCSC genome browser annotations showed that this variant fell in regulatory regions (Supplementary Figure 6).

FIGURE 4

Figure 4 Top functional IVs. Heatmap of top 7 functional IVs. The Sei score is the fold change of each variant in 40 sequence classes. The b score is the b value in the sensitivity analysis of each variant. The SNP rs4409785 has the largest Sei score and b score > 0.5. The SNP rs3850765 has the largest b score and absolute Sei score > 1.

FIGURE 5

Figure 5 Functional annotations of SNP rs4409785. The top 1 track is the density track of ReMap 2022. The top 2 to 7 tracks are CTCF ChIP-seq bigwig peaks in immune cells. The top 8 to 11 tracks are ChIA-pet data in GM12878, T cell, activated T cell, and B cell. ChIP-seq, chromatin immunoprecipitation sequencing; ChIA-pet, chromatin interaction analysis by paired-end tag sequencing.

Discussion

In this study, two independent MR analyses were conducted by using European ancestry GWAS data to detect the relationship between hypothyroidism and RA risk, and found that hypothyroidism has a causal effect on RA risk. Moreover, because levothyroxine sodium is used to treat hypothyroidism, we also found that the genetic data of people who have hypothyroidism have a causal effect on RA risk. To validate the conclusion, we adopted larger datasets. Based on the FinnGen R8 release, the relationship of hypothyroidism, drug reimbursements between hypothyroidism and risk of RA in European ancestry is significant with the IVW method (OR = 1.98; 95% CI 1.34, 2.92; p < 0.001) (Supplementary Figure 7A). The results of Non-cancer illness code self-reported: hypothyroidism/myxoedema (OR = 2.13; 95% CI 1.52, 2.97; p < 0.001) and Treatment/medication code: levothyroxine sodium (OR = 24.92; 95% CI 7.05, 88.15; p < 0.001) with RA in the larger dataset also showed a significant relationship between these phenotype (Supplementary Figures 7B, C). The results of the other four methods also suggested a significant relationship, which validated the conclusion in the main datasets. To investigate whether the conclusion of hypothyroidism and RA reached from those with European ancestry also applies to Asians, a replication TSMR study using Asian GWAS data was also conducted. Though rs1049087 (F statistic = 30. 08) is the only significant lead SNP in the GWAS of hypothyroidism after filtering by conditions, we used the Wald ratio method but still obtained the same conclusion (OR = 2.12; 95% CI 1.79, 2.51; p < 0.001).

A reverse direction of hypothyroidism and RA TSMR analysis showed that there was no significant causal relationship between RA and hypothyroidism with the MR-Egger method. However, the IVW method showed a significant causal relationship between RA and hypothyroidism with no significant horizontal pleiotropic effect (Supplementary Table 3). Previous studies have reported that patients with RA are at high increased risk of thyroid dysfunctions (31). However, the earlier diagnosis does not mean which one starts earlier. Nisihara et al. (32) reported that rheumatic autoantibodies have been detected from 17.5% of autoimmune thyroid disease patients without clinical rheumatic disorder for 5 years. It suggested the need to closely monitor hypothyroidism patients so that RA can be diagnosed early.

MR is now widely used to assess the causal relationship through genetic markers between exposure and outcome (33). For the coding variants of IVs, we adopted PolyPhen-2 to predict the effect. SNP rs3775291 was predicted to be "likely damaged", which fell into the gene toll-like receptor 3 (TLR3) to change the amino acid leucine to phenylalanine. TLR3 is highly expressed in the thyrocytes of Hashimoto’s disease patients and may affect the innate immune response (34). Moreover, toll-like receptors play complex roles in the pathogenesis of RA; arthritis was ameliorated when interference targeted TLR3 in rats (35), which suggested the important role of TLR3 in both thyroid disease and RA. Though the genetic markers always come from GWAS data, the functional molecular mechanism remained elusive since 90% of disease-related GWAS SNPs fall into the noncoding regions. The complex regulatory code in noncoding regions made it difficult to decipher the key SNPs and mechanisms. We found that the SNP rs4409785 locus is more activated in the activated CD8-positive alpha-beta T cell in CTCF-ChIP data. CTCF transcription factor binding sites are often found next to the binding sites of thyroid hormone receptors (36). Studies in human samples suggested that the CD8-positive T cell plays an important role in the initiation and maintenance of the disease process (37). A new strategy was proposed here; we mined the causal effects deeper in IVs. In our study, we integrated the effects of IVs evaluated by MR analysis, noncoding variant predicted tool, and functional genomic annotations to pinpoint the potential functional noncoding variants in hypothyroidism and RA risk.

Although we obtained significant results from the TSMR methods, there were still several limitations that should be noted. First of all, the traits of hypothyroidism are derived from the public GWAS dataset with hypothyroidism, drug reimbursement, and self-reported hypothyroidism with the non-cancer illness code, rather than the diagnosis code; the more specific diagnosis code of the GWAS dataset of hypothyroidism will increase the statistical power. Second, the number of GWAS data of hypothyroidism, especially the case number, is small. The statistical power will be increased with GWAS data with a larger sample size. Third, it has been reported that women and elderly RA patients have an increased risk of developing hypothyroidism (31), but the GWAS summary data in our study were unable to adjust for gender and age because the information was not available. Finally, the result from MR reflects the change in RA risk due to a lifelong change in exposure hypothyroidism status, rather than a specific time in life (38). Therefore, the short-term effect size from MR analysis merits additional investigation and should not be regarded as a short-term intervention.

Conclusions

In summary, in this study, we provided strong evidence of hypothyroidism having a causal effect on RA risk. Importantly, we provided the potential critical SNPs in IVs and biological functional genomics annotations of these variants. Although translating findings into treatment required further validation, we still provided evidence that hypothyroidism patients need to be monitored for RA risk, which might help in the early detection of RA.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. All code is available on Github (https://github.com/shify-nxu/MR-RA). Further inquiries can be directed to the corresponding author.

Author contributions

YG: Methodology, software, data curation, formal analysis, visualization, writing—original draft, and writing—review and editing. Z-RF: Data curation, formal analysis, and writing—review and editing. F-YS: conceptualization, project administration, supervision, funding acquisition, resources, and writing—review and editing. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by the Key Research and Development Program of Ningxia (Special Talents) (grant numbers: 2022BSB03042 and 2022BSB03043), and the Natural Science Foundation of Tianjin Municipal Science and Technology Commission (grant number: 22JCQNJC00850).

Acknowledgments

We want to acknowledge the participants and investigators of FinnGen study and UK Biobank collaborators.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2023.1179656/full#supplementary-material

Abbreviations

GWAS, genome-wide association studies; RA, rheumatoid arthritis; TSMR, two-sample Mendelian randomization; IVs, instrument variants; OR, odds ratio; CI, confidence interval; SE, standard error; SNP, single-nucleotide polymorphism; MR, Mendelian randomization; LD, linkage disequilibrium; IVW, inverse variance weighted; CTCF, CCCTC-binding factor; ESR1, estrogen receptor alpha gene; AR, androgen receptor; FOXA1, recombinant forkhead box protein A1; MAF, minor allele frequency; VCF, variant call format; ChIP-seq, chromatin immunoprecipitation sequencing; ChIA-pet, chromatin interaction analysis by paired-end tag sequencing.

References

1. Smolen JS, Aletaha D, McInnes IB. Rheumatoid arthritis. Lancet (2016) 388:2023–38. doi: 10.1016/S0140-6736(16)30173-8

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ge X, Frank-Bertoncelj M, Klein K, McGovern A, Kuret T, Houtman M, et al. Functional genomics atlas of synovial fibroblasts defining rheumatoid arthritis heritability. Genome Biol (2021) 22:247. doi: 10.1186/s13059-021-02460-6

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Klareskog L, Padyukov L, Lorentzen J, Alfredsson L. Mechanisms of disease: genetic susceptibility and environmental triggers in the development of rheumatoid arthritis. Nat Clin Pract Rheumatol (2006) 2:425–33. doi: 10.1038/ncprheum0249

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Li Q, Wang B, Mu K, Zhang J, Yang Y, Yao W, et al. Increased risk of thyroid dysfunction among patients with rheumatoid arthritis. Front Endocrinol (Lausanne) (2018) 9:799. doi: 10.3389/fendo.2018.00799

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Conigliaro P, D’Antonio A, Pinto S, Chimenti MS, Triggianese P, Rotondi M, et al. Autoimmune thyroid disorders and rheumatoid arthritis: a bidirectional interplay. Autoimmun Rev (2020) 19:102529. doi: 10.1016/j.autrev.2020.102529

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Antonelli A, Ferrari SM, Corrado A, Di Domenicantonio A, Fallahi P. Autoimmune thyroid disorders. Autoimmun Rev (2015) 14:174–80. doi: 10.1016/j.autrev.2014.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Fallahi P, Ferrari SM, Ruffilli I, Elia G, Biricotti M, Vita R, et al. The association of other autoimmune diseases in patients with autoimmune thyroiditis: review of the literature and report of a large series of patients. Autoimmun Rev (2016) 15:1125–8. doi: 10.1016/j.autrev.2016.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Ferrari SM, Fallahi P, Ruffilli I, Elia G, Ragusa F, Benvenga S, et al. The association of other autoimmune diseases in patients with graves’ disease (with or without ophthalmopathy): review of the literature and report of a large series. Autoimmun Rev (2019) 18:287–92. doi: 10.1016/j.autrev.2018.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Gubbi S, Hannah-Shmouni F, Verbalis JG, Koch CA. Hypophysitis: an update on the novel forms, diagnosis and management of disorders of pituitary inflammation. Best Pract Res Clin Endocrinol Metab (2019) 33:101371. doi: 10.1016/j.beem.2019.101371

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Wu D, Xian W, Hong S, Liu B, Xiao H, Li Y. Graves’ disease and rheumatoid arthritis: a bidirectional mendelian randomization study. Front Endocrinol (Lausanne) (2021) 12:702482. doi: 10.3389/fendo.2021.702482

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Mahagna H, Caplan A, Watad A, Bragazzi NL, Sharif K, Tiosano S, et al. Rheumatoid arthritis and thyroid dysfunction: a cross-sectional study and a review of the literature. Best Pract Res Clin Rheumatol (2018) 32:683–91. doi: 10.1016/j.berh.2019.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

12. McCoy SS, Crowson CS, Gabriel SE, Matteson EL. Hypothyroidism as a risk factor for development of cardiovascular disease in patients with rheumatoid arthritis. J Rheumatol (2012) 39:954–8. doi: 10.3899/jrheum.111076

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Liu Y, Miao H, Lin S, Chen Z. Association between rheumatoid arthritis and thyroid dysfunction: a meta-analysis and systematic review. Front Endocrinol (Lausanne) (2022) 13:1015516. doi: 10.3389/fendo.2022.1015516

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Smith GD, Ebrahim S. “Mendelian randomization”: can genetic epidemiology contribute to understanding environmental determinants of disease? Int J Epidemiol (2003) 32:1–22. doi: 10.1093/ije/dyg070

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Farh KK-H, Marson A, Zhu J, Kleinewietfeld M, Housley WJ, Beik S, et al. Genetic and epigenetic fine mapping of causal autoimmune disease variants. Nature (2015) 518:337–43. doi: 10.1038/nature13835

PubMed Abstract | CrossRef Full Text | Google Scholar

16. FinnGen consortium. R5 release of FinnGen consortium genome-wide association analysis data (2021). Available at: https://finngen.gitbook.io/documentation/.

Google Scholar

17. Stahl EA, Raychaudhuri S, Remmers EF, Xie G, Eyre S, Thomson BP, et al. Genome-wide association study meta-analysis identifies seven new rheumatoid arthritis risk loci. Nat Genet (2010) 42:508–14. doi: 10.1038/ng.582

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Sakaue S, Kanai M, Tanigawa Y, Karjalainen J, Kurki M, Koshiba S, et al. A cross-population atlas of genetic associations for 220 human phenotypes. Nat Genet (2021) 53:1415–24. doi: 10.1038/s41588-021-00931-x

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Kurki MI, Karjalainen J, Palta P, Sipilä TP, Kristiansson K, Donner KM, et al. FinnGen provides genetic insights from a well-phenotyped isolated population. Nature (2023) 613:508–18. doi: 10.1038/s41586-022-05473-8

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Okada Y, Wu D, Trynka G, Raj T, Terao C, Ikari K, et al. Genetics of rheumatoid arthritis contributes to biology and drug discovery. Nature (2014) 506:376–81. doi: 10.1038/nature12873

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Papadimitriou N, Dimou N, Tsilidis KK, Banbury B, Martin RM, Lewis SJ, et al. Physical activity and risks of breast and colorectal cancer: a mendelian randomisation analysis. Nat Commun (2020) 11:597. doi: 10.1038/s41467-020-14389-8

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bowden J, Davey Smith G, Burgess S. Mendelian randomization with invalid instruments: effect estimation and bias detection through egger regression. Int J Epidemiol (2015) 44:512–25. doi: 10.1093/ije/dyv080

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Burgess S, Butterworth A, Thompson SG. Mendelian randomization analysis with multiple genetic variants using summarized data. Genet Epidemiol (2013) 37:658–65. doi: 10.1002/gepi.21758

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Burgess S, Thompson SG. Interpreting findings from mendelian randomization using the MR-egger method. Eur J Epidemiol (2017) 32:377–89. doi: 10.1007/s10654-017-0255-x

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Bowden J, Davey Smith G, Haycock PC, Burgess S. Consistent estimation in mendelian randomization with some invalid instruments using a weighted median estimator. Genet Epidemiol (2016) 40:304–14. doi: 10.1002/gepi.21965

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Hartwig FP, Davey Smith G, Bowden J. Robust inference in summary data mendelian randomization via the zero modal pleiotropy assumption. Int J Epidemiol (2017) 46:1985–98. doi: 10.1093/ije/dyx102

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Walker VM, Davies NM, Hemani G, Zheng J, Haycock PC, Gaunt TR, et al. Using the MR-base platform to investigate risk factors and drug targets for thousands of phenotypes. Wellcome Open Res (2019) 4:113. doi: 10.12688/wellcomeopenres.15334.1

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Chen KM, Wong AK, Troyanskaya OG, Zhou J. A sequence-based global map of regulatory activity for deciphering human genetics. Nat Genet (2022) 54:940–9. doi: 10.1038/s41588-022-01102-2

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, et al. A method and server for predicting damaging missense mutations. Nat Methods (2010) 7:248–9. doi: 10.1038/nmeth0410-248

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Almandoz JP, Gharib H. Hypothyroidism: etiology, diagnosis, and management. Med Clin North Am (2012) 96:203–21. doi: 10.1016/j.mcna.2012.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Huang C-M, Sung F-C, Chen H-J, Lin C-C, Lin C-L, Huang P-H. Hypothyroidism risk associated with rheumatoid arthritis. Med (Baltimore) (2022) 101:e28487. doi: 10.1097/MD.0000000000028487

CrossRef Full Text | Google Scholar

32. Nisihara R, Pigosso YG, Prado N, Utiyama SRR, De Carvalho GA, Skare TL. Rheumatic disease autoantibodies in patients with autoimmune thyroid diseases. Med Princ Pract (2018) 27:332–6. doi: 10.1159/000490569

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Dong S-S, Zhang K, Guo Y, Ding J-M, Rong Y, Feng J-C, et al. Phenome-wide investigation of the causal associations between childhood BMI and adult trait outcomes: a two-sample mendelian randomization study. Genome Med (2021) 13:48. doi: 10.1186/s13073-021-00865-3

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Harii N, Lewis CJ, Vasko V, McCall K, Benavides-Peralta U, Sun X, et al. Thyrocytes express a functional toll-like receptor 3: overexpression can be induced by viral infection and reversed by phenylmethimazole and is associated with hashimoto’s autoimmune thyroiditis. Mol Endocrinol (2005) 19:1231–50. doi: 10.1210/me.2004-0100

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Meng L, Zhu W, Jiang C, He X, Hou W, Zheng F, et al. Toll-like receptor 3 upregulation in macrophages participates in the initiation and maintenance of pristane-induced arthritis in rats. Arthritis Res Ther (2010) 12:R103. doi: 10.1186/ar3034

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Lutz M. Thyroid hormone-regulated enhancer blocking: cooperation of CTCF and thyroid hormone receptor. EMBO J (2003) 22:1579–87. doi: 10.1093/emboj/cdg147

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Carvalheiro H, da Silva JAP, Souto-Carneiro MM. Potential roles for CD8+ T cells in rheumatoid arthritis. Autoimmun Rev (2013) 12:401–9. doi: 10.1016/j.autrev.2012.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Davies NM, Holmes MV, Davey Smith G. Reading mendelian randomisation studies: a guide, glossary, and checklist for clinicians. BMJ (2018) 362:k601. doi: 10.1136/bmj.k601

PubMed Abstract | CrossRef Full Text | Google Scholar