Alterations and Mechanism of Gut Microbiota in Graves’ Disease and Hashimoto’s Thyroiditis

Graves’ disease (GD) and Hashimoto’s thyroiditis (HT) are the most common autoimmune thyroid disease (AITD) types (Fallahi et al. 2019; Rayman 2019; Knezevic et al. 2020). AITD is a set of organ-specific autoimmune diseases with similar genetic and immunological features. The pathogenesis of AITD depends on multiple factors, but the exact mechanism is still unclear. It is generally believed that genetic susceptibility, environmental and survival factors (gender difference), stress, and other factors have important roles (Ajjan and Weetman 2015; Yoo and Chung 2016; Banga and Schott 2018). Moreover, recent evidence has suggested that the gut microbiota is closely associated with some immune-related diseases, including type 1 diabetes mellitus (T1DM) (Kugelberg 2017), rheumatoid arthritis (RA) (Lynch and Pedersen 2016; Horta-Baas et al. 2017), multiple sclerosis, Graves’ ophthalmopathy (Covelli and Ludgate 2017; Shi et al. 2019), HT, and inflammatory bowel disease (Masetti et al. 2018; Zhao et al. 2018; Kozhieva et al. 2019). The gut microbiota has a crucial role in the metabolism, absorption, immune function, and defense mechanism against pathogens (Pickard et al. 2017; Azad et al. 2018; Reddel et al. 2019). Nevertheless, the exact effect of gut microbiota on AITD, particularly HT and GD, is still not well defined.

It has been suggested that gut microbiota targets the TSH receptor (Knezevic 2020; Yao et al. 2020). The combination of microbe and thyroid autoantibody suggests that it may have a role in AITD (Kristensen 2016). Therefore, a deep understanding of the exact mechanism behind these changes and their relationship with AITD may help develop new prevention and treatment strategies. This study explored the alterations and putative activity of gut microbiota in GD and HT. Fecal samples from the GD, HT patients, and healthy people were collected and analyzed using 16s rRNA sequencing.

Study population. Twenty-seven GD patients, 27 HT patients, and 16 healthy people were included in the study. The demographic and clinical parameters of the subjects are summarized in Table I.

The gut microbiota abundance and diversity in the GD and HT groups were similar to those in the healthy groups, but the overall structure was different. To identify whether the GD and HT were associated with changes in microbiota diversity, we sequenced and analyzed fecal samples. Thirteen phyla, 23 classes, 43 orders, 75 families, 221 genera, 422 species, and 595 operational taxa (OTU) were found in the GD group; 12 phyla, 21 classes, 33 orders, 61 families, 201 genera, 394 species and 585 out in the HT group; and 13 phyla, 22 classes, 35 orders, 64 families, 180 genera, 322 species, and 436 OTU were in the control group, all of which had 97% similarity. According to OTU analysis results, the grade-abundance curves of the GD and HT patients and the healthy control group presented similar patterns (Fig. 1A and 1B). The results showed that the richness and diversity of gut microbiota in the healthy control group tended to be lower than those in the GD and HT, but the differences were not significant. According to the Sobs and Simpson index in PAN/Core species analysis, alpha diversity analysis, and a Shannon index and dilution curve where both species richness and uniformity are considered, the species abundance, total species, and core species obtained by sequencing were sufficient. Consequently, the sample sequencing quantity was considered satisfactory, indicating the results were convincing.

Fig. 1

The dilution curve analysis showed that the gut microbiota of the GD and HT patients had a similar species richness compared to the healthy group. A total of 686 OTUs were detected in all the samples, among which 389 were commonly shared among groups. Sixty-three, 61, and 21 unique OTUs were identified in the GD, HT, and healthy control samples.

Next, taxon-dependent analysis was performed using the Ribosome Database Project (RDP) classifier to describe gut microbiota composition in different groups. The HT group had the highest content of Proteobacteria and Actinomycetes, followed by the GD group and the healthy control group. Notably, the HT group contained a small number of Verrucomicrobiaceae. At the level of “family” and “genus”, the community composition of the GD and HT group were different compared with that of the healthy control group (Fig. 1C and 1D). Moreover, PlS-DA analysis using a binary Jaccard similarity algorithm showed that the overall microbial composition of the GD group and HT group was somewhat different from that of the healthy control group, but there was no significant difference between the GD group and HT group (Fig. 1E). In addition, ANOSIM showed that gut microbiota composition was significantly different between the GD group, HT group, and healthy control group, with an R-value of 0.2519 (Fig. 1F).

Fig. 1

Fig. 1

These results indicated that the levels of bacterial abundance and diversity in the gut microbiota of the GD and HT patients were similar to those of the healthy controls. In contrast, the overall structure of the gut microbiota of both patients and healthy controls were significantly different.

The abundance of gut microbiota in the GD and HT groups. The linear discriminant effect size (LEfSe) method was used to identify specific bacteria associated with GD and HT. The branching diagram showed fecal microflora and major bacterial structures in the healthy controls and the GD and HT patients and compared the most considerable taxonomic group differences between the two communities (Fig. 2A). The LEfSe analysis revealed 24 discriminant features of class (n = 3), order (n = 3), family (n = 4), and genus (n = 14) (linear discriminant analysis LDA > 3, p < 0.05) when comparing the GD group and the control group, and 13 discriminant features of class (n = 2), order (n = 3), family (n = 3) and genus (n = 5) when comparing the HT group and the control group (linear discriminant analysis LDA > 3, p < 0.05).

Fig. 2

The abundance of Negativicutes in healthy control samples and Proteobacteria and Erysipelotrichia in GD patient samples increased. Coriobacteriaceae and Erysipelotrichia were more abundant in HT patient samples than in other samples (Fig. 2A).

At the “phylum” level, the proportions of Cyanobacteria in the GD samples were higher than those in the healthy control samples, while the proportions of abnormal cocci and Cyanobacteria were lower (Fig. 2B). Moreover, the proportions of Cyanobacteria in the samples of the HT patients were higher than that of the healthy control group, while the proportions of abnormal Coccinobacteria and Cyanobacteria were lower (Fig. 2C).

Fig. 2

At the level of “family”, Lachnospiraceae, Alcaligenaceae, Christensenellaceae, and Erysipelotrichaceae were prevalent in the GD patient samples (Fig. 2D); Enterococcaceae, Erysipelotrichaceae, and Bacillobacteriaceae were abundant in the HT patient samples, while Peptostreptococcaceae, Bacillaceae, and Matophyaceae were high in the healthy control samples (Fig. 2E).

Fig. 2

At the level of “genera”, Prevotella_9, Ruminococcus_2, and Lachnospiraceae_NK4A136_group were higher in GD patient samples, while the proportion of Megamonas genus was abundant in the healthy controls (Fig. 2F). The contents of Ruminococcus_2 and Enterococcus in the HT patient samples were relatively high (Fig. 2G).

These data suggest a difference in the microbiome of GD and HT patients compared to the healthy control group. Although there was no significant change in bacterial diversity, the abnormal composition of fecal microflora indicated gut microbiota imbalance in the GD and HT patients.

Bacillus, Blautia, and Ornithinimicrobium can be used as potential markers to distinguish GD and HT patients from the healthy people. Next, a random forest analysis was performed to compare the GD group (Fig. 3A) or HT group (Fig. 3B) with a healthy control group. The results of random forest analysis showed that the areas under the verification information curve of the top three strains of Bacillus, Blautia, and Ornithinimicrobium in the GD patients and the top two strains of Bacillus and Ornithinimicrobium in the HT patients were 0.98842 and 1, respectively (Fig. 3C and 3D). It suggested that these bacteria may be valuable markers for distinguishing healthy patients from GD and HT patients and could be used as potential diagnostic markers of GD and HT.

Fig. 3

Fig. 3

Functional categories according to the COG and KEGG in different groups. We also predicted the functional categories according to the COG and KEGG in different groups; the results are shown in Fig. 4. According to the COG distribution in Fig. 4A and 4B, the GD and HT groups were highly enriched in carbohydrate transport and metabolism (function G) compared to the control group, while the amino acid transport and metabolism (function E) was lower than that of the healthy control group. The function of “translation, ribosome structure, and biogenesis” (function J) was highest in the GD group, followed by the healthy control group, being lowest in the HT group. As this revealed a significant difference between the disease group and the healthy control group, it should be a focus of future research. Interestingly, the “transcription” (function K) of GD and HT groups was significantly stronger than the healthy control group.

Fig. 4

According to the KEGG distribution in Fig. 4C, there were significant differences in purine metabolism, aminoacyl tRNA biosynthesis, cysteine, and methionine metabolism between the GD group and the healthy control group (all p < 0.05). Moreover, there were significant differences in ribosome and pyrimidine metabolism between the HT and healthy control groups (Fig. 4D). According to the KEGG-based results in Fig. 4C and 4D, the ABC transporter, responsible for the ATP transport pathway, was significantly more abundant in the disease group than in the healthy control group.

Fig. 4

According to the COG database, the enzyme “glycosyltransferase” was a specific enzyme in the GD group (Fig. 4E). Also, the resolving enzyme has been suggested as a specific enzyme for the HT group (Fig. 4F).

According to the KEGG database, the “ATP-dependent RNA helicase DHX58” (EC 3.6.3.14) was the highest in the GD group, followed by the healthy control group, while it was the lowest in the HT group (Fig. 4G and 4H). “Glutamine synthase” (EC 6.3.1.2) (Fig. 4G) and “DNA-directed RNA polymerase B subunit” (EC 2.7.7.6) (Fig. 4H) were the specific enzymes for the GD group and HT group, respectively.

Fig. 4

Fig. 4

Ten different strains of the two groups were divided into three categories, as shown in Fig. 5. Table II shows the top ten predictions using the KEGG database for the abundance of these three categories. The metabolic pathway of the “ABC transporter” (responsible for ATP transport) existed in the prediction results of the three different strains, indicating that this metabolic pathway is highly correlated with the occurrence of GD and HT.

Fig. 5

Next, we conducted a Venn diagram analysis based on the differential strains in Fig. 3A and 3B. Ten differential strains were further divided into two groups and three categories, as shown in Fig. 5. Table II shows the prediction results according to the KEGG database of the top ten strains different in abundance. The “ABC transporter” pathway (responsible for ATP transport) was found in the predicted results of three different strains, suggesting that this pathway was strongly associated with the development of the GD and HT.

GEO database screening for the HT differential genes according to the KEGG database. We first cross-referenced the differentially expressed genes in the two chips and then annotated them with the KEGG database. Next, the annotation results were compared to the predicted results for different strains using KEGG. The four metabolic pathways common in the two predicted results were glutathione metabolism, arachidonic acid metabolism, purine metabolism, and pyrimidine metabolism. Therefore, we suggest that four strains unique to the HT (Ruminococcus_1, Flavonifractor, Moryella, and Anaerotruncus) may affect the occurrence and development of HT by regulating glutathione metabolism and arachidonic acid metabolism. The six common bacteria (Bacillus, Corynebacterium, Ornithinimicrobium, Brachybacterium, Nocardioides, and Ruminococcus_gnavus_group) can participate in the occurrence and development of the HT by regulating purine metabolism and pyrimidine metabolism.

GD and HT are two major representative diseases of AITD. A previous study suggested an association between gut microbiota imbalance and HT or GD (Virili et al. 2018; Yao 2020). However, thus far, no studies have reported a common imbalance of gut microbiota in GD and HT patients. Our study found multiple bacteria with similar change direction and shared metabolic pathways involved in the GD and HT patients.

In this study, genomic DNA was extracted from the GD, HT, and healthy subjects feces and analyzed using the 16S rRNA gene sequencing. We found that the abundance and diversity of gut microbiota in the GD and HT patients were similar to the healthy control group. However, we also discovered that Actino-bacteria and Proteobacteria contents were the highest in the HT group, followed by the GD group, while they were the lowest in the control group. A previous retrospective study showed the highest alteration in the abundance of Bacteroidetes, Proteobacteria, and Firmicutes between the systemic inflammatory disease group and the healthy group (Nam et al. 2013; Clemente et al. 2018; Faucher et al. 2020). Zhao et al. (2018) found a higher gut microbiota richness and diversity in HT patients with normal thyroid function. Firmicutes were the most abundant, while Bacteroides were less common in HT patients, consistent with our findings. Nevertheless, in this study, HT patients were all hypothyroidism patients, that are different from the study reported by Zhao et al. (2018). Furthermore, Zhou et al. (2014) showed the gut microbiota diversity in the GD patients; Bifidobacteria and Lacto-bacillus were significantly reduced, but Clostridium and Enterococcus were increased compared to the healthy groups. Other studies reported on Helicobacter pylori and Yersinia enterocolitica, mainly focusing on the relationship between H. pylori CagA and AITD (Köhling et al. 2017; Figura et al. 2019; Cuan-Baltazar and Soto-Vega 2020). Bassi et al. (2012) and Soveid et al. (2012) suggested that H. pylori is associated with GD, but not with HT (Bassi et al. 2014), while Wenzel et al. (1988) found that IgA and IgG anti-Yersinia antibodies were significantly increased in GD and HT. Moreover, Takuno et al. (1990) found that Y. enterocolitica was significantly correlated with GD, but not with HT. Effraimidis et al. (2011) reported no causal relationship between Y. enterocolitica infection and the AITD. Another study showed an increased number of Bacteroidetes and Proteobacteria and decreased levels of Firmicutes in GD patients, which is consistent with our data (Ebert 2010; Köhling et al. 2017). The same study suggested a higher abundance of Pasteurellaceae and Prevotella in GD patients compared to healthy people (Ebert 2010; Köhling et al. 2017); yet, this was not observed in our study. ANOSIM analysis showed significant differences in gut microbiota composition between the HT, GD, and the healthy group, which further indicated that the intestinal microecology of the HT and GD patients was unbalanced.

The LEfSe analysis showed similar trends in bacteria in the GD and HT groups, where the most apparent changes included an increased abundance of Erysipelotrichia, Cyanobacteria, Ruminococcus_2, and decreased abundance of Bacillaceae and Megamonas. These data suggest that there may be a common gut microbiota disorder in GD and HT patients. A previous study (Kozhieva et al. 2019) demonstrated the increased abundance levels of Blautia, Roseburia, Ruminococcus_torques_group, Romboutsia, Dorea, Fusicatenibacter, and Eubacterium_hallii_group, while the Fecalibacterium, Bacteroides, Prevotella_9, and Lachnoclostridium decreased in HT patients. Other studies (Liu et al. 2020) suggested that Lachnospiraceae_incertae_sedis, Lactonifactor, Alistipes, and Subdoligranulum were more enriched in HT patients with euthyroidism, while Phascolarctobacterium was more abundant in those with hypothyroidism. Yan et al. (2020) found that the number of Bacilli, Lactobacillales, Prevotella, Megamonas, and Veillonella strains in GD patients increased, while Ruminococcus, Rikenellaceae, and Alistipes decreased compared with the healthy people. Our data are consistent with the above results, suggesting that the gut microbiota of GD and HT patients is in an unbalanced state. However, our study confirmed that there are species with the similar trend change in HT and GD patients, confirming that there may be a common imbalance of flora involved in the occurrence of GD and HT.

The results of random forest analysis indicated that the areas under the verification information curve of Bacillus, Blautia, and Ornithinimicrobium in the GD patients and Bacillus and Ornithinimicrobium in the HT patients were equal to 0.98842 and 1, respectively, suggesting that these strains may be used as biomarkers to distinguish healthy individuals from the GD and HT patients. In order to analyze the differences in microbial composition between the disease group and the healthy control group, we used the PlS-DA analysis based on a binary Jaccard to replace the traditional PCoA analysis, as the PlS-DA analysis adds grouping information. Through group supervision, this method can ignore the random differences within the groups and highlight the systematic differences between the groups, which is more illustrative than PCoA and other methods.

In this study, the COG database was used to predict the function of the HT, GD, and healthy control groups. We found that HT and GD patients were enriched in carbohydrate transport and metabolism compared to the control group but had lower amino acid transport and metabolism activity. According to the COG prediction of the function of different strains, four unique strains of the HT group and six common strains of the HT and GD groups were identified as “S” (Function unknown). It may indicate that there is still an unknown metabolic pathway in the development of such diseases, and it might be the focus of subsequent research. In this study, KEGG was used to predict the functions of differential bacteria unique to the GD and HT groups and the bacteria common to these two groups. Venn diagram analysis was simultaneously performed. We found that the “ABC transporter” metabolic pathway existed in the predicted results of three different strains, indicating that this metabolic pathway was highly correlated with the occurrence of GD and HT. Studies have shown that the levels of L-arginine, L-ornithine, lysine, and guanbutamine in the GD and HT patients are higher than those in the healthy group, while the levels of putrescine, 1,3-diaminopropylene, spermine, and N-acetylputrescine are lower than those in the healthy group (Song et al. 2019). Some polyamine metabolites were only different in the GD or HT patients compared with the healthy group. Spermidine proportions were significantly reduced in all the patients. This study confirmed that most metabolites of the GD and HT had similar patterns compared with the healthy group, suggesting a common pathophysiological basis or metabolic pathway.

This study has a few limitations. First, it was a single-center study with a relatively small sample size, which may lead to bias. Second, the specific mechanism of abnormal gut microbiota involved in GD and HT was not examined. Third, all the people were selected from inland North China in this study. Considering the differences in dietary structure and ethnicity, the higher iodine intake in coastal areas, thyroid function, and gut microbiota may differ. Thus, further in-depth multicenter studies with a large sample size should be carried out to confirm these findings.

Our study demonstrated that Bacillus, Blautia, and Ornithinimicrobium could be used as potential markers to distinguish the GD and HT from the healthy population and the “ABC transporter” metabolic pathway may be involved in the pathogenesis of the GD and HT. In our future work, we plan to construct an animal model of the GD and HT, after which a fecal bacteria transplantation or intervention of differential metabolites will be conducted. It aims to clarify further the regulatory role of gut microbiota in the occurrence and development of GD and HT, and to explore the specific mechanism of abnormal gut microbiota involved in GD and HT. Our study suggested that GD and HT patients have similar changes in gut microbiota and the same metabolic pathway usage, suggesting that there may be common changes in gut microbiota in AITD patients, which needs to be further studied.