Alterations of Gut Microbiota in Patients With Graves’ Disease

Graves’ disease (GD) is a systemic autoimmune disease characterized by hyperthyroidism. Evidence suggests that alterations to the gut microbiota may be involved in the development of autoimmune disorders. The aim of this study was to characterize the composition of gut microbiota in GD patients. Fecal samples were collected from 55 GD patients and 48 healthy controls. Using 16S rRNA gene amplification and sequencing, the overall bacterial richness and diversity were found to be similar between GD patients and healthy controls. However, principal coordinate analysis and partial least squares-discriminant analysis showed that the overall gut microbiota composition was significantly different (ANOSIM; p < 0.001). The linear discriminant analysis effect size revealed that Firmicutes phylum decreased in GD patients, with a corresponding increase in Bacteroidetes phylum compared to healthy controls. In addition, the families Prevotellaceae, and Veillonellaceae and the genus Prevotella_9 were closely associated with GD patients, while the families Lachnospiraceae and Ruminococcaceae and the genera Faecalibacterium, Lachnospira, and Lachnospiraceae NK4A136 were associated with healthy controls. Metagenomic profiles analysis yielded 22 statistically significant bacterial taxa: 18 taxa were increased and 4 taxa were decreased. Key bacterial taxa with different abundances between the two groups were strongly correlated with GD-associated clinical parameters using Spearman’s correlation analysis. Importantly, the discriminant model based on predominant microbiota could effectively distinguish GD patients from healthy controls (AUC = 0.825). Thus, the gut microbiota composition between GD patients and healthy controls is significantly difference, indicating that gut microbiota may play a role in the pathogenesis of GD. Further studies are needed to fully elucidate the role of gut microbiota in the development of GD.


INTRODUCTION
Graves' disease (GD) is an autoimmune disorder in which stimulatory autoantibodies activate thyroid stimulating hormone receptors (TSHR), leading to excessive amounts of thyroid hormones. It is the most common cause of hyperthyroidism, with typical symptoms including weight loss, fatigue, heat intolerance, tremor, and palpitations. GD can occur at all ages, but is more common in women than in men (Ehlers et al., 2019). The prevalence of hyperthyroidism in the United States is approximately 1.2% (0.5% overt and 0.7% subclinical) (Doubleday and Sippel, 2020). The prevalence rate of overt hyperthyroidism estimated by a population-based study in Taiwan is around 0.27-0.37%, with an annual incidence of 0.97 to 1.06 cases per 1,000 individuals (Kornelius et al., 2018). The pathogenesis of GD remains unclear, but it is hypothesized that infections, iodine intake, smoking, and mental stress play a role, and both genetic and environmental factors appear to be involved (Shukla et al., 2018).
The human gut microbiota harbors up to 100 trillion microbes. These microbes play a crucial role in the maintenance of human health by affecting nutrient digestion, immune response, metabolic homeostasis, and protection against pathogen colonization (Lin and Zhang, 2017;Pickard et al., 2017). An imbalance in the gut microbiota or dysbiosis can increase epithelial barrier dysfunction and lead to intestinal and systemic disorders. Recent studies have established a link between the alterations in the composition of gut microbiota in many disease states, including inflammatory bowel disease, Crohn's disease, obesity, and type 2 diabetes (Ortega et al., 2020;Szilagyi, 2020). In addition, the components and metabolites derived from the gut microbiota are key regulators of host immune responses. Alterations in the gut microbiota may result in an imbalance of T-cell subpopulations, triggering the release of proinflammatory and anti-inflammatory cytokines (Lee and Kim, 2017). Moreover, this can trigger the onset of several types of autoimmune diseases, including type 1 diabetes, rheumatoid arthritis, and systemic lupus erythematosus (Alkanani et al., 2015;Correa et al., 2017;Picchianti-Diamanti et al., 2018).
The thyroid gland is the primary endocrine gland in the human body. However, the relationship between the gut microbiota and GD remains unclear. The aim of this study was to characterize the changes in the gut microbiota of patients with GD using 16S rRNA gene sequencing, compare it with the gut microbiota of their healthy counterparts, and determine its possible correlation with the relevant clinical factors.

Ethics Statement
This study's protocol was approved by the Institutional Review Board and Ethics Committee of Chang Gung Memorial Hospital (IRB 201701806B0C502). Written informed consent was obtained from all participants before the start of the study. All experimental procedures were performed in accordance with the approved guidelines.

Study Population and Sample Collection
A total of 55 GD patients were recruited from the Division of Endocrinology and Metabolism at Chang Gung Memorial Hospital between October 2017 and March 2020. These patients were previously diagnosed with GD and had an average follow-up of 45.33 months. All patients were initially treated with anti-thyroid drugs, including propylthiouracil (PTU), methimazole (MMI), or carbimazole (CBZ). GD was diagnosed on the basis of clinical presentations and laboratory data, including symptoms and signs of thyrotoxicosis, diffuse goiter, the presence of ophthalmopathy, elevated free T4 and suppressed TSH levels, and the presence of autoantibody against TSH receptor. In addition, 48 age-, sex-, and BMI-matched healthy subjects from a health screening center constituted the control group. The control group had normal values of free T4, TSH, and thyroid peroxidase (TPO) antibody, and no history of thyroid disease. Free T4 and TSH levels were measured using a one-step chemiluminescent immunoassay (CLIA) on the ADVIA Centaur XP system (Siemens Healthcare, Munich, Germany). TPO antibody and TSH receptor antibody were measured by chemiluminescent microparticle immunoassay (CMIA) (Architect, Abbott, Wiesbaden, Germany). Reference intervals were defined as follows: free T4 0.70-1.48 ng/dL; TSH 0.35-5.50 mIU/mL; anti-TPO <5.61 IU/mL; anti-TSHR <1.75. A questionnaire was used to assess dietary habits and health condition of all subjects. In each group, subjects who met any of the following criteria were excluded from this study: pregnancy; gastrointestinal disorders (chronic diarrhea or inflammatory bowel disease); gout; stroke; cancer; autoimmune diseases (type 1 diabetes, systemic lupus erythematosus, or rheumatoid arthritis); history of gastrointestinal surgery; use of antibiotics, probiotics, prebiotics, symbiotics (<2 months), hormonal medication or Chinese herbal medicine (<3 months); and pure vegetarian. The fecal samples were collected in a clean container, aliquoted, and all aliquots were immediately frozen and stored at -80°C.

16S rRNA Gene Sequencing
DNA was extracted using the QIAamp Fast DNA Stool Mini Kit (Qiagen, Germany) following the manufacturer's protocol. Amplification of the 16S V3-V4 region was performed using the universal primers reported previously (Klindworth et al., 2013). Samples were amplified using 2× KAPA HiFi HotStart Ready Mix (KAPA Biosystems, MA, USA), and indexed using the Nextera XT 96 Index Kit (Illumina, CA, USA). Lastly, the library was sequenced on an Illumina MiSeq platform and 300 bp paired-end reads were generated.

Bioinformatics and Statistical Analysis
The overlapping and merging of paired-end reads were conducted in FLASH (v1.2.11) (Magoc and Salzberg, 2011). Filtering of the raw tags was performed under specific filtering conditions according to the QIIME (v1.9.1) (Caporaso et al., 2010). Chimeras were removed by aligning sequences to the "gold" reference database using UCHIME (Edgar et al., 2011). The resulting filtered sequences were clustered into operational taxonomic units (OTUs) using UPARSE (v7.0.1090) with a 97% similarity threshold. The representative sequences were then classified taxonomically using the ribosomal database project (rdp) classifier (v2.2) with default 0.8 as the confidence threshold. The OTUs were aligned using PyNAST (v1.2) and assigned to taxonomic groups against the Greengenes database (v13.8) or SIlVA database (v132) . The alpha diversity metrics were estimated in QIIME (v1.9.1), including the coverage percentage (Good's), Venn diagram, richness estimators (ACE and Chao1), and diversity indices (Shannon and Simpson). Rank-abundance and rarefaction curves were calculated at a level of 97% similarity of the OTUs. Beta diversity measures, such as principal coordinate analysis (PCoA) and partial least squares-discriminant analysis (PLS-DA), were calculated using QIIME (v1.9.1). The bacterial taxonomic distributions of the sample communities were visualized using R software. Venn diagrams were plotted using the VennDiagram package in R. Significant differences in the microbiota structure between the two groups were evaluated by nonparametric analysis of similarity (ANOSIM) using the vegan package in R. The specific characterization of the gut microbiota that differentiate among the groups was analyzed using a linear discriminant analysis (LDA) effect size (LEfSe) algorithm (Segata et al., 2011). The parameters set with default p-value, a = 0.05, and an LDA score of 4.0 with LEfS. The relative abundances of OTUs between the GD group and the healthy controls were compared using Statistical Analysis of Metagenomic Profiles (STAMP) software (v2.1.3), Welch's t-test (p < 0.05), and Benjamini-Hochberg False Discovery Rate (FDR) method (q < 0.05) as a multiple test correction (Parks et al., 2014). The correlation analysis between clinical parameters and different microbiota was performed using Spearman rank correlation, where a p-value < 0.05 was considered statistically significant. Random forest algorithm was used to evaluate the prediction accuracy of the gut microbiota of GD patients and healthy controls. The receiving operational curve (ROC) was analyzed and the area under curve (AUC) using ten-fold cross-validation was calculated. Anthropometric and clinical characteristics were analyzed using SPSS version 18.0 (SPSS Inc., Chicago, USA).

Study Population
Fifty-five subjects were enrolled into the GD group and 48 volunteers were enrolled into the healthy control group. The demographic details of the subjects are summarized in Table 1.

Diversity Indices and Operational Taxonomic Units (OTUs) Assays of the Gut Microbiota
A total of 11,718,443 usable sequences were obtained from 103 samples using the Illumina MiSeq platform. From these, 5,783,276 high-quality sequences were selected, with an average of 56,158 sequences per barcode sample. Specifically, 671 OTUs in the healthy controls and 684 OTUs in the GD group were delineated at a 97% similarity level. The diversity of the gut microbiota observed in each group is shown in Table 2, and detailed characteristics of each sample are shown in Supplementary Table S1. The values of Good's coverage for the two groups were more than 99%, indicating that the current sequencing depth was sufficient to saturate the bacterial diversity of the gut microbiota. The bacterial richness estimates according to ACE and Chao were not significantly different. Analysis of bacterial community richness and evenness with Shannon and Simpson indices showed that the microbiota diversity of GD group, as a whole, was similar to that of the healthy controls. Although both the Shannon and Simpson indices were higher in the GD group, the differences were not statistically significant (p > 0.05). Based on the results of the OTU analysis, the rankabundance curves showed a similar pattern of bacterial communities in the two groups ( Figure 1A). Rarefaction curve analysis also indicated that the species richness of the GD group exhibited a tendency similar to that of the healthy controls ( Figure 1B). To better characterize the shared richness within  The coverage percentage (Good's), richness estimators (ACE and Chao1), and diversity indices (Shannon and Simpson) were calculated using Good's method and the QIIME pipeline (v1.9.1), respectively. C.I., confidence interval.  the two groups, the unique and overlapping OTU data were presented in a Venn diagram. These diagrams showed that 589/ 766 OTUs were shared between the two groups. Thus, 95 and 82 unique OTUs were identified in the GD group and healthy controls, respectively (Supplementary Figure S1).

Differences in Gut Microbiota Composition Between GD and Healthy Controls
To further analyze the composition of the gut microbiota in the different experimental groups, taxonomy from phylum to genus level was assigned to each OTU (Supplementary Figure S2). As a result, the relative abundance of the phyla Bacteroidetes and Actinobacteria was found to increase, while Firmicutes decreased in the GD group compared to the healthy controls ( Figure 1C).
In addition, at the genus level, the abundances of Bacteroides and Prevotella_9 were found to be significantly higher, while Faecalibacterium and Lachnospiraceae_NK4A136_group were slightly lower in the GD group compared to the healthy controls ( Figure 1D). To measure similarities between the microbial communities, beta diversity analysis was performed using principal coordinate analysis (PCoA) and partial least squares-discriminant analysis (PLS-DA). Despite significant inter-individual variations, the gut microbiota from GD group were clearly different from those of the healthy controls ( Figure  1E). A clear separation between groups was also observed in PLS-DA ( Figure 1F). The results of analysis of similarities (ANOSIM) indicated that the structure of the gut microbiota significantly differed between the two groups (ANOSIM, r = 0.158, p < 0.001). This comparison revealed that the bacterial diversity of the gut microbiota in the GD group was similar to that observed in the healthy controls. However, the overall community structure was distinctive between the two sample groups.
To identify the specific bacterial taxa associated with GD patients, the gut microbiota in the healthy controls and GD patients was compared using the LEfSe method. A cladogram representing the structure of the gut microbiota and the predominant bacteria in the healthy controls and GD patients is shown in Figure 2A, in which the greatest differences in the taxa between the two communities are displayed. LEfSe analysis revealed 16 discriminative features, including 2 phyla, 3 classes, 3 orders, 4 families, and 4 genera (LDA score >4; Figure 2B). Bacteria belonging to the Bacteroidetes phylum were enriched in the GD patients, while those in Firmicutes were more abundant in the healthy controls. At the family level, the proportion of Prevotellaceae and Veillonellaceae was increased in GD patients, while that of Lachnospiraceae and Ruminococcaceae was decreased compared to that in healthy controls. Moreover, the level of Prevotella_9 was significantly higher in GD patients; however, the levels of Faecalibacterium, Lachnospira, and Lachnospiraceae_NK4A136_group were markedly increased in the healthy controls. Thus, these microbial features may be used as potential biomarkers for distinguishing GD patients from healthy controls. The changes in intestinal microbiota composition from the phylum to the species level were further explored by statistical analysis of metagenomic profiles (STAMP) analysis using Welch's t-test. STAMP analysis yielded 22 statistically significant bacterial taxa (p < 0.05, q < 0.05) ( Table  3). At the level of phyla, Bacteroidetes and Actinobacteria were significantly more abundant in the gut microbiota of GD patients than in healthy controls, whereas Firmicutes was less abundant. At the family level, Prevotellaceae, Erysipelotrichaceae, Tannerellaceae, Coriobacteriaceae, and Actinomycetaceae were prevalent in GD patients, while the family Lachnospiraceae was found to be significantly higher in the healthy controls. At the genus and species levels, the abundance of the genera Prevotella_9, Parabacteroides, and Collinsella, as well as those of the species Actinomyces_odontolyticus, was higher in the GD patients than in the healthy controls ( Figure 2C).

Correlation Between Clinical Parameters and the Gut Microbiota
To explore the association between the clinical parameters and gut microbiota, a Spearman's correlation analysis was conducted among five clinical parameters (TSH, TPOAb, FT4, age, and BMI) with 29 bacterial taxa, which were differentially abundant between the GD group and healthy controls in the LEfSe and STAMP analyses (Figure 2). The abundance levels of 16 GDenriched taxa were positively correlated with TPOAb but negatively correlated with TSH, while six taxa enriched in the healthy controls showed contrasting correlations (p < 0.05) ( Figure 3A). Among these, 12 out of 16 GD-enriched taxa were also positively correlated with FT4, while three of the six control-enriched taxa were also negatively correlated with FT4 (p < 0.05). However, there were no significant correlations between the gut microbiota and other environmental factors, including age and BMI (p > 0.05). A heatmap was generated based on the 15 most abundant taxa from the 29 different taxa. Hierarchical clustering showed a clear separation between the GD group and the healthy controls ( Figure 3B). Further, to evaluate the predictive power of these taxa in predicting the GD status, the random forest analysis was used and an AUC value of 0.825 was achieved ( Figure 3C), indicating that these microbial features may be used as potential biomarkers for distinguishing GD patients from healthy controls.

DISCUSSION
In this study, the gut microbiota of GD patients was compared with those of age-and sex-matched healthy subjects using 16S rRNA sequencing technologies. Patients with GD were found to have an altered gut microbiota composition compared to the healthy controls. In fact, differentially abundant bacterial taxa responsible for the differences between the two groups were identified. In addition, Spearman correlation analysis confirmed that several taxa were significantly correlated with clinical parameters. Among these, the top 15 most abundant taxa discriminated GD patients with high accuracy (AUC = 0.825).
Fecal samples from GD patients were found to show a slight increase in microbial diversity and richness compared to those of the healthy controls, although the differences were not statistically significant ( Table 2). Previous studies have observed bacterial diversity in hyperthyroid and hypothyroid patients, which may be associated with bacterial overgrowth in the intestinal tract (Lauritano et al., 2007;Zhou et al., 2014). Recently, increased microbial diversity has also been reported in patients with Hashimoto's thyroiditis and thyroid carcinoma (Zhao et al., 2018;Feng et al., 2019).
PCoA and PLS-DA analyses demonstrated that the microbial compositions of GD patients and healthy controls were divided into two separate groups; the GD group exhibited a relatively high homology, indicating that GD patients may share common microbiota features (Figures 1E, F). In our study, GD patients had significantly higher levels of Bacteroidetes and Actinobacteria and lower levels of Firmicutes compared to the healthy controls ( Figure 2 and Table 3). These findings are consistent with those of previous studies in patients with GD and Graves' orbitopathy (Shi et al., 2019;Jiang et al., 2021). However, in other studies of patients with Hashimoto's thyroiditis or thyroid cancer, Firmicutes was found to be more abundant and Bacteroidetes less abundant (Zhao et al., 2018;Feng et al., 2019). The ratio of Firmicutes to Bacteroidetes has been considered as an index of the health of the gut microbiota and has been correlated with obesity (Mariat et al., 2009;Tilg and Kaser, 2011). However, a recent meta-analysis suggested that the way in which the gut microbiota modulates obesity is more complex than a simple imbalance in these commensal phyla.
In addition, our study showed that the fecal samples of the GD patients were rich in members of the species Actinomyces_odontolyticus, the genus Collinsella, Parabacteroides, and Prevotella_9, as well as of the family Veillonellaceae ( Figures  2B, C). Actinomyces odontolyticus is normally present in the oral cavity and gastrointestinal tract of healthy humans. This species has been shown to be involved in the formation of biofilms on teeth (Liljemark et al., 1993), and has been associated with several human diseases, including actinomycosis, periodontitis, and colorectal cancer (Ximenez-Fyvie et al., 1999;Riegert-Johnson et al., 2002;Yachida et al., 2019). Although its role in the pathogenesis of inflammatory bowel disease has not been demonstrated, it may change the enteric environment and immune factors and aggravate Taxonomic cladogram and linear discriminant analysis (LDA) scores show differences among taxa between GD and healthy controls. Only taxa meeting a significant LDA threshold value of >4 are shown (A, B). Differentially abundant taxa from the phylum to genus level were further analyzed by STAMP analysis using Welch's ttest (p < 0.05, q < 0.05) (C).
injuries caused by inflammation . The genus Collinsella within the family Coriobacteriaceae has been previously linked to rheumatoid arthritis and has showed a strong correlation with high levels of alpha-aminoadipic acid and asparagine, as well as the production of the proinflammatory cytokine IL-17A (Chen et al., 2016). Parabacteroides can be considered an antiinflammatory symbiont (Kverka et al., 2011;Hiippala et al., 2020), increased levels of which have been observed in patients with active Behcȩt's disease (Ye et al., 2018). Prevotella is a known producer of propionate, which can induce the differentiation of regulatory T cells and suppress Th17 polarization (Li et al., 2016). The increased abundance of Prevotella has been associated with GD, rheumatoid arthritis and systemic lupus erythematosus (De Aquino et al., 2014;Mendonca et al., 2019;Yan et al., 2020). It may also influence the efficacy of drug therapy for GD (Yan et al., 2020). Our study found significantly high levels of Prevotella_9 in GD patients, which is consistent with previous reports in GD patients (Ishaq et al., 2018), while decreased levels of Prevotella_9 have been found in HT patients (Zhao et al., 2018). Besides, our study showed that the abundance of the family Veillonellaceae in patients with GD was increased, which is consistent with the previous studies (Yan et al., 2020;Chen et al., 2021). The Veillonella genus is normal flora of the mouth, gastrointestinal tract, and vagina. Recent studies have shown that Veillonella species were enriched in the fecal microbiota of primary sclerosing cholangitis patients (Little et al., 2020) and may contribute to the immune system development during early childhood (Poppleton et al., 2017).
In addition, our data showed significantly lower levels of the genus Faecalibacterium, Lachnospiraceae NK4A136 group, and Lachnospira in GD patients compared to the healthy controls ( Figures 2B, C). Lachnospiraceae NK4A136 group and Lachnospira belong to the family Lachnospiraceae, which can digest carbohydrates to produce butyrate, and which has been reported to have potent anti-inflammatory effects (Flint et al., 2012). Previous studies have reported a decrease in the abundance of the genus Lachnospira in patients with chronic kidney disease and childhood asthma (Stiemsma et al., 2016;Lun et al., 2019). The relative abundance of Lachnospiraceae NK4A136 group has been shown to decrease in children with cystic fibrosis (Coffey et al., 2019). The decrease in Faecalibacterium has been linked with GD patients (Ishaq et al., 2018;Zhao et al., 2018;Cornejo-Pareja et al., 2020) and other autoimmune diseases like multiple sclerosis and Crohn's disease (Thorkildsen et al., 2013;Cantarel et al., 2015). A previous research showed that Faecalibacterium was a protective factor, owing to its negative correlation with thyroid stimulating immunoglobulin antibody related to GD patients (Cornejo-Pareja et al., 2020).
Moreover, in addition to investigating alterations of the gut microbiota, we further investigated the relationship between gut microbiota and clinical parameters. TSH, TPOAb, and FT4 were found to be highly correlated with microbiota composition, indicating a close relationship between the gut microbiota and GD ( Figure 3A). Hierarchical clustering based on the relative abundance of the top 15 taxa from 29 gut microbiota also showed a clear separation between the two groups ( Figure 3B), and an AUC value of 0.825 was achieved ( Figure 3C), further indicating that some gut microbiota may be associated with GD. The potential role of gut microbiota in the development of GD warrants further investigation.
In summary, we demonstrated that GD patients, compared to healthy controls, exhibited gut microbiota dysbiosis characterized by unchanged overall bacterial richness and diversity, but an altered taxonomic composition. Several specific differential bacterial taxa were found to be significantly correlated with GD. Nevertheless, this study has several limitations. First, this was a single-center, cross-sectional study involving a limited number of patients. Second, whether the altered gut microbiota is a cause or a consequence of GD development will need to be confirmed. Lastly, gut microbiota studies using 16S rRNA sequencing have yet to provide a detailed list of functional genes and metabolic pathways for the study of host-gut microbiome interactions. Therefore, further multicenter studies with larger samples using metagenomics and metabolomics integrative analyses are needed to improve our understanding of the interactions between gut microbes, their metabolites, and host cells.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/, PRJNA695309.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board and Ethics Committee of Chang Gung Memorial Hospital (IRB 201701806B0C502). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
S-CC, S-FL, and J-JL designed the experiments. S-FL and S-TC collected the samples and clinical data. P-YC and F-SL helped supervise the project. Y-MY performed the statistical analyses. S-CC and J-JL wrote the manuscript. All authors discussed the results and commented on the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by grants from the Chang Gung Memorial Hospital (CMRPG3F1693).

SUPPLEMENTARY MATERIAL
The