This study is designed to determine the causal effect of gut microbiota on hypothyroidism and further explore through which pathway the screened taxon of microbiota causes hypothyroidism. MR analyses (Burgess and Thompson, 2021) are used to provide evidence for the causal association between microbiota and steps in the clinical pathway to hypothyroidism.

Study A is designed to screen the potential taxa of the gut microbiota that are causally related to hypothyroidism (Figure 1A). In Study A, two families and genera of gut microbiota that can causally affect the occurrence of hypothyroidism are investigated. Since the rigorous correction for multiple comparisons might neglect potential taxa that are causally associated with hypothyroidism, multiple testing is not performed in Study A. Screening is done according to three criteria. The first criterion is significant at both family and genus levels, which can increase the robustness and reliability of the results of the selected microbiota. Although it might wrongly remove some potentially valuable microbiota, it can better make sequent results correct. The second criterion is better to be validated by multiple MR methods. Each MR method has its own strengths. With significant results of more MR methods, the sequent analyses become more reliable. The third criterion is, if there are probiotics and harmful microbiota screened at the same time, the probiotics will be chosen as the final candidate target. Because the supplement of a bacteria is easier and more economical than targeting a bacterium, the findings will have higher clinical values. The screened key taxon of the gut microbiota is then used in Study B to determine more specific mechanisms.

In Study B, the screened “key taxon” is used to investigate whether it has a causal effect on known causative factors of hypothyroidism with MR analyses (Figure 1B). Then, the former known causative factors are further investigated to identify a causally affected clinical pathway. By confirming the linear causal effect of “key taxon” at each step of the known clinical pathway, the application potential of the identified taxon of gut microbiota can be proven.

All data used in this study are summary-level data of genome-wide association studies (GWAS). The source of all data is summarized in Supplementary Table S1.

The data of gut microbiota are available from the MiBioGen database (https://mibiogen.gcc.rug.nl/), which was provided by Kurilshikov et al. (2021). These comprehensive GWAS data of gut microbiota were collected from a pool of 18,340 individuals from 24 cohorts. The data of 165 taxa (including 35 families and 130 genera) are finally used in the present study. The details about the measurement methods can be checked in previously published studies by Kurilshikov et al. (2021). Supplementary Table S2 lists all the taxa of gut microbiota used in this study (in the column “exposure”).

Fecal samples were collected from 18,340 people, the majority of whom are Europeans with a median age of approximately 55 years. After fecal sample collection from each cohort, DNA purification kits were used to extract the DNA. The MiSeq technology was used for DNA sequencing in most cohorts, while the HiSeq2500 technology was used for one cohort. Due to differences in the 16S rRNA gene domains among the cohorts, the author directly mapped the reads to the SILVA ribosomal RNA database for taxonomic classification. A cutoff posterior probability of 0.8 was applied at each taxonomic level to assign reads to their respective taxa. This posterior cutoff probability was assessed independently for each taxonomic level. Only taxa present in over 10% of the samples within each cohort were left. For the comprehensive mapping of multi-cohort Quantitative Trait Loci, the study-wide criteria were an effective sample size of a minimum of 3,000 samples and presence in at least three cohorts. The single nucleotide polymorphism (SNP) panels of cohorts were mainly provided by Illumina and Affymetrix. After quality control, SNPs strongly associated with one of the different taxa (P < 0.0001), as well as their corresponding positions, alleles, regression results, and sample size, all of which were recorded and uploaded to the MiBioGen database used in this study.

The outcome variables are all collected from the IEU OpenGWAS project (https://gwas.mrcieu.ac.uk/), which was built by the University of Bristol, and contains summary-level GWAS data from multiple databases. For example, the outcome variable of Study A, hypothyroidism, is derived from UK Biobank but downloaded from IEU OpenGWAS project.

In the investigation of factors contributing to hypothyroidism, the majority of data are derived from FinnGen database, which is built by the University of Helsinki, since it divided hypothyroidism into different subtypes, such as autoimmune hypothyroidism and hypothyroidism due to medicaments and other exogenous factors. However, in the exploration of infectious variables related to subacute thyroiditis, there are no data about coxsackievirus and adenovirus in the FinnGen database. Instead, the summary-level GWAS data of coxsackievirus and adenovirus receptor provided by Sun et al. (2018) are used. Details can be checked in each source database (Supplementary Table S1).

The instrument variable (IV) selection is done according to the standardized processes of MR analyses using P-value threshold and linkage disequilibrium (LD) clumping. To obtain the adequate number of genetic variations as IVs and retain IVs strongly associated with each taxon of gut microbiota, the significant threshold is finally set at 1 × 10⁻⁵, which is also a common threshold in many high-quality MR studies of gut microbiota (Li et al., 2022; Song et al., 2023). Then, SNPs with LD are also removed using a threshold of r² < 0.001 and a maximum distance of 10,000 kb. Finally, 2,105 independent SNPs associated with 165 bacterial traits are identified as IVs.

This study adhered to the three core assumptions of MR, as outlined by Bowden and Vansteelandt (2011): (1) the IVs are strongly associated with the taxa of gut microbiota, (2) IVs are independent of confounders affecting outcome variables, and (3) the IVs affect the outcome variable solely through the exposure and not via alternative pathways. The first assumption is met in the selection of SNPs, while the last two are assessed by pleiotropy tests.

In this study, five MR analytical methods are used, namely, inverse variance-weighted (IVW), MR Egger, weighted median, weighted mode, and simple median. IVW is the method with the greatest statistical power but biased if there is horizontal pleiotropy. On the contrary, the MR Egger method is robust to horizontal pleiotropy. The weighted median method is more robust to outliers than IVW and useful when significant result is found in the heterogeneity test. The weighted mode and simple median methods are supplementary to the weighted median method. Therefore, the primary methods of MR analyses are IVW when no horizontal pleiotropy or heterogeneity is found, MR Egger when there is horizontal pleiotropy, and weighted median when there is heterogeneity but no horizontal pleiotropy (Burgess and Thompson, 2021). The Supplementary material contains more detailed introduction of these analytical methods of MR and the way and rationale for the method selection.

MR analysis relies on Mendel's laws of inheritance and instrumental variable estimation methods, which are designed to calculate causal effects in the presence of unobserved confounding (Sanderson et al., 2022). To address potential confounders, heterogeneity and pleiotropy have to be estimated and used for the selection of the robust MR methods. Heterogeneity refers to differences and outliers in SNP-specific causal estimates. Cochran's Q statistic is calculated to evaluate heterogeneity in this study, which is a weighted sum of the squared distances of the SNP-specific causal estimates from the overall IVW estimate. The term pleiotropy is sometimes used interchangeably with “horizontal pleiotropy,” while “mediation” is sometimes called “vertical pleiotropy.” Horizontal pleiotropy indicates that IVs can affect outcomes through covariate rather than only through the exposure, while vertical pleiotropy indicates exposure affect outcome via the covariate. Horizontal pleiotropy is detected by the MR Egger method, while horizontal pleiotropy can also be detected robustly.

In addition, sensitivity analyses with the leave-one-out method are performed to assess the reliability, stability, and robustness of the causal estimates, specifically investigating whether the results were disproportionately influenced by any individual SNP. The calculation of all analyses was done with the “TwoSampleMR” package (version 0.5.6) on R (version 4.2.2). The results of MR analyses are presented with odds ratio (OR) and confidence intervals (CIs). With a significance level below P of < 0.05, it can be concluded that there might be a causal relationship.

Since the MR analysis of clinical outcomes does not reveal the molecular mechanisms, we attempted to carry out enrichment analyses of the key taxon IVs, although these results might be only suggestive rather than confirmed.

First, the single nucleotide polymorphism database (dbSNP) provided by the National Institutes of Health of the United State was utilized to determine the related genes of IVs (https://www.ncbi.nlm.nih.gov/snp/). Then, the tool known as the database for annotation, visualization and integrated discovery (DAVID) provided by the National Institutes of Health of the United States was used to perform gene-enrichment analysis with annotation categories of diseases (from Genetic Association Disease Database) and gene ontology (Huang and Sherman, 2009; Sherman et al., 2022) (https://david.ncifcrf.gov/). The default threshold used is P < 0.1 and count >2.

We attempted to use the enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, but it produced no significant results. As an alternative, the transcription factor genes that interact with the products of IV-related genes were first investigated to indirectly explore the possible related molecular pathways. The interacting transcription factors were also investigated using the DAVID tool. Then, the enrichment analysis of KEGG pathways was performed again with the identified transcription factors to ascertain the possible molecular pathways involved (It is worth noting that these indirect results are suggestive and have much weaker power).

A total of 2,151 SNPs were selected as IVs for 35 families and 130 genera of gut microbiota, which is listed in Supplementary Table S2. Five MR methods mentioned before were performed for each taxon, and all results are displayed in Supplementary Table S3, while the results of tests of heterogeneity and horizontal pleiotropy are listed in Supplementary Table S4. The directions of each pair of taxon and outcome were further confirmed with the Steiger test to avoid reverse causation (Supplementary Table S6).

For easy check, the significant results of primary methods are summarized and displayed in Table 1, while the corresponding results of heterogeneity and horizontal pleiotropy are listed in Supplementary Table S5. The relationships of significant families and genera are shown in Figure 2. The figure shows that two families and two genera can causally inhibit hypothyroidism, while two families and six genera can facilitate hypothyroidism. Only Defluviitaleaceae (OR 0.043, 95%CI 0.005–0.363, P = 0.018)/Defluviitaleaceae_UCG_011 (OR 0.385, 95%CI 0.172–0.865, P = 0.021) and Bacteroidaceae (OR 3.363, 95%CI 1.312–8.617, P = 0.012)/Bacteroides (OR 3.363, 95%CI 1.312–8.617, P = 0.012) are significant at both the taxonomical levels.

Although the primary method results of both pairs are significant, more methods of Defluviitaleaceae also support significance. Meanwhile, a supplementary of probiotic bacteria is more economical and has higher application value than targeting harmful bacteria. Therefore, the Defluviitaleaceae family and Defluviitaleaceae_UCG_011 genus are selected as key taxa of gut microbiota for further investigation in Study B.

In this part, SNPs used as IVs is listed in Supplementary Table S7. All of the MR results in this part are summarized in Supplementary Table S8, while results of tests of heterogeneity and horizontal pleiotropy are listed in Supplementary Table S9. The directions of causation between Defluviitaleaceae and outcomes are checked with Steiger test, and there is no reverse causation (Supplementary Table S11). The MR results of outcomes, which contain significant results are listed in Table 2 for easy check, and corresponding results of heterogeneity and pleiotropy are in Supplementary Table S10.

Firstly, variables related to hypothyroidism in FinnGen database are selected as outcome variables to understand the causal effect of Defluviitaleaceae. These factors include “thyroiditis,” “congenital iodine-deficiency syndrome,” “autoimmune hypothyroidism,” “postinfectious hypothyroidism,” “hypothyroidism due to medicaments and other exogenous,” and “other/unspecified hypothyroidism.” The results of thyroiditis are significant for both of Defluviitaleaceae family (OR 0.174, 95%CI 0.046–0.653, P = 0.029) and Defluviitaleaceae_UCG_011 genus (OR 0.139, 95%CI 0.029–0.664, P = 0.043) (Figures 3A, B), while no significance is found in other outcome variables.

Then, MR analyses are performed between Defluviitaleaceae and three subtypes of thyroiditis in FinnGen database, which are “subacute thyroiditis” (family: OR 0.028, 95%CI 0.004–0.213, P = 0.007; genus: OR 0.018, 95%CI 0.002–0.194, P = 0.013) (Figures 3A, B), “autoimmune thyroiditis” (family: OR 0.982, 95%CI 0.459–2.104, P = 0.964; genus: OR 0.916, 95%CI 0.391–2.146, P = 0.840), and “unspecified thyroiditis” (family: OR 1.401, 95%CI 0.632–3.103, P = 0.406; genus: OR 2.254, 95%CI 0.926–5.487, P = 0.073) (Supplementary Table S8).

Therefore, the association of Defluviitaleaceae with mumps (family: OR 0.811, 95%CI 0.460–1.429, P = 0.469; genus: OR 0.791, 95%CI 0.419–1.492, P = 0.468) (Supplementary Table S8) and influenza (family: OR 0.818, 95%CI 0.676–0.989, P = 0.038; genus: OR 0.792, 95%CI 0.644–0.974, P = 0.027) (Figures 3A, B) are further explored, which are the common infections causing subacute thyroiditis. Because there are no summary-level GWAS data of coxsackievirus and adenovirus, the key proteins, coxsackievirus and adenovirus receptor, are analyzed as an alternative (family: OR 0.99,4 95%CI 0.818–1.209, P = 0.956; genus: OR 1.045, 95%CI 0.847–1.289, P = 0.679) (Supplementary Table S8). Although there is no causality found with coxsackievirus and adenovirus receptor in this study, it needs more direct evidence as proof.

Finally, MR analyses are performed between Defluviitaleaceae and antibody level of different subtypes of influenza. A significant increase of H3N2 IgG levels with Defluviitaleaceae family (OR 1.934, 95%CI 1.123–3.332, P = 0.017) was observed, while the result with Defluviitaleaceae_UCG_011 genus is not significant (OR 1.675, 95%CI 0.953–2.943, P = 0.073) but still suggests an increase of the H3N2 antibody level (Figures 3A, B). Both the family and genus fail to find significance with H1N1 antibody level (family: OR 1.260, 95%CI 0.705–2.255, P = 0.435; genus: OR 1.210, 95%CI 0.651–2.249, P = 0.547) (Supplementary Table S8).

In sum, the results of the present study suggest that Defluviitaleaceae might be able to increase anti-influenza H3N2 IgG levels to decrease the occurrence of influenza. The decrease of influenza by Defluviitaleaceae prevents subacute thyroiditis and hypothyroidism (Figure 3C). Figure 4 shows the significant findings in the form of scatter plots.

As supplementary, forest plots (showing robustness of MR results) and leave-one-out plots (showing any outliers in IVs) of variables are displayed in Supplementary Figures S1, S2.

There are 13 different SNPs as IVs of Defluviitaleaceae, of which 9 of them have corresponding genes, which are listed in Table 3A. Since only nine genes were used as input to perform enrichment analysis, the results with P < 0.1 were limited, and all of them are listed in Table 3B. The results show involvement in metabolic and immune fields, specifically with neutrophils, respiratory function tests, celiac disease, lipoproteins, and echocardiography results. In GO analysis, these genes are associated with cell adhesion. Although it is difficult to observe the association between cell adhesion and influenza, thyroiditis, and hypothyroidism, the enrichment results indicating immune system and neutrophils count are consistent with the finding of this study.

Since the nine SNP-related genes failed to find the involved pathway in KEGG pathway enrichment analysis, the transcription factors interacting with them were identified instead, which are listed in Supplementary Table S12A. The KEGG pathways of these transcription factors may give us some hints about the involved molecular pathways. The results can be observed in Supplementary Table S12B. It is worth noting that these indirect results are suggestive and have much weaker power. Since it is the indirect results of interacting transcription factors, the ranges get larger and are involved with the FoxO and JAK-STAT signaling pathway, immunity (e.g., Th17 cell differentiation), viral infection (e.g., measles and viral hepatitis), hormones (e.g., prolactin), metabolism, and cancer. The immunity, viral infection, and hormones are exactly consistent with the findings of the former part of this study.