The Potential Link between Gut Microbiota and Serum TRAb in Chinese Patients with Severe and Active Graves' Orbitopathy

Background and Objective A previous study reported alterations in the intestinal microbiota in patients with Graves' orbitopathy (GO). Thyrotropin receptor autoantibody (TRAb) stimulates orbital and periorbital tissues and plays a pivotal role in the development of GO. However, the association between gut microbiota and TRAb in GO patients has still remained elusive. In this study, we explored the relationships between gut microbiota and GO-related traits, in which we applied a metabolic-network-driven analysis to identify GO trait-related modules and extracted significant operational taxonomic units (OTUs). Methods In the present study, we profiled gut microbiota using 16S rRNA gene sequencing in 31 GO patients. We performed metabolic-network-driven analysis to investigate the association between gut microbiota and GO-related traits (e.g., TRAb, TGAb, and TPOAb) in the combination of microbial effects. Results Applying microbiome network analysis of cooccurrence patterns and analysis of topological properties, we found that s_Prevotella_copri and f_Prevotellaceae showed a significant correlation with TRAb. In particular, we applied the latent class model to explore the association between gut microbiota and GO-related traits in the combination of microbial effects. It was revealed that the subjects involved in the latent class model with the higher abundance of s_Prevotella_copri and g_Bacteroides had a higher TRAb level. Conclusions Our results revealed the potential relationships between gut microbiota and GO-related traits in the combination of microbial effects. This study may provide a new insight into the interaction between the intestinal microbiota and TRAb-associated immune responses in GO patients.


Introduction
Graves' orbitopathy (GO) is an autoimmune disease, commonly associated with Graves' disease (GD). It influences appearance, visual acuity, and even quality of life of patients [1][2][3]. To date, the pathogenesis of GO has not been fully understood. Gut microbiota influences various autoimmune diseases, such as type 1 diabetes (T1D) [4] and systemic lupus erythematosus (SLE) [5]. A recently conducted study demonstrated that gut microbiota is associated with some thyroid diseases, including Hashimoto's thyroiditis (HT) and GD [6]. Increased intestinal permeability and infiltration of intraepithelial lymphocytes have been previously found in patients with HT [7]. Intestinal symbiotic microorganisms may influence extraintestinal immune responses inducing loss of tolerance to self-antigens, including thyroglobulin that underlies HT [8]. Our previous study revealed alterations in the intestinal microbiota in patients with severe and active GO. For instance, community diversity was significantly reduced in patients with GO. At the phylum level, the proportion of Bacteroidetes was notably increased, while at the genus and species levels, significant differences were observed [9]. e present study highlighted the role of microbiota in occurrence of GO.
To our knowledge, evaluating the autoimmune inflammation of GO patients into different stages, active and nonactive forms, is highly significant for the treatment of GO. An evidence showed that thyrotropin receptor autoantibody (TRAb) stimulates orbital and periorbital tissues and also plays a pivotal role in the development of GO; thus, detection of TRAb may be of clinical significance in active assessment of disease [10,11]. Elevated expression of the thyroid-stimulating hormone (TSH) receptor in orbital tissues supports the substantial role of TRAb in the pathogenesis of GO [10,12]. Recently, Kahaly et al. reported that high titers of TRAb are related to thyroid-associated orbitopathy in patients with GD [13]. However, to date, a limited number of studies explored the role of gut microbiota in GO patients, while none of the studies reported a potential association between gut microbiota and TRAb in such patients.
In the present research, to explore the relationships between gut microbiota and GO-related traits, e.g., TRAb and TGAb, we applied a metabolic-network-driven analysis to identify GO trait-related modules and extract important operational taxonomic units (OTUs). We identified some novel associations between gut microbiota and GO-associated traits. Our study provided a framework to better perceive the interactions of gut microbiota and extracted the important bacteria associated with TRAb.

Study Subjects.
In the current study, the 16S rRNA gene sequence was used to reconstruct the taxonomic structure of gut microbial communities using the fecal samples of GO patients. is study was performed at the Department of Endocrinology, Beijing Tongren Hospital, Capital Medical University, Beijing, China. Between March 2017 and March 2018, 31 patients with severe and active GO with hyperthyroid were enrolled. All patients received only an antithyroid drug ( yrozol; Merck & Co., Inc., Kenilworth, NJ, USA). e diagnosis of GO was carried out according to the European Group on Graves' Orbitopathy (EUGOGO) guidelines [2]. e enrolled GO patients had not received any treatment for ocular discomfort. e active GO was defined by a clinical activity score (CAS) ≥3/7, and severe GO was defined by the NOSPECS score ≥IV. TRAb was measured using a commercially available electrochemiluminescence assay kit (Roche Diagnostics GmbH, Mannheim, Germany) based on the M22 monoclonal antibody, with a normal range <1.75 U/L. e exclusion criteria for the GO patients were patients' age <18 or >65 years, history of chronic diarrhea or constipation, history of gastrointestinal surgery, therapy of probiotics or antibiotics over the previous 4 weeks, use of hormonal medication (<3 months), severe disease (acute infections, diabetes, stroke, renal or hepatic dysfunction, cancer, or autoimmune diseases), pure vegetarian, etc. [9]. e characteristics of the study subjects are summarized in Table 1 GTGCCAGCMGCCGCGGTAA; 806R: GGACTACHVGG GTWTCTAAT) were amplified using specific primers [14]. e sequences were analyzed and those with ≥97% similarity were classified into the same OTUs. e representative sequence for each OTU was screened out for further annotation. e abundance information on the OTUs was normalized using a standard sequence number, corresponding to each sample with the least number of sequences as previously described [9].

Network Analysis
e composition of bacterial communities could be positively influenced, in addition to negative relationships between contributing microorganisms in human disease. us, intermicrobial relationships can be inferred from the cooccurrence network of taxa, and this network can be used to investigate ecological interactions between microbes [15]. In the present research, we identified cooccurrence networks of gut microbiota based on Spearman's correlation analysis using the relative abundance tables, and P values were adjusted for multiple testing using the Benjamini-Hochberg approach. ese networks kept those correlations with the adjusted P value <0.1, and the absolute value of correlation coefficients was >0.3. By applying the modularity scores, we attempted to find out a dense subgraph. We applied the igraph package of R software (http://www.r-project.org) to implement the analysis. (WGCNA). Herein, the WGCNA was used to generate the network and identify network modules based on the OTU relative abundance. In WGCNA, we kept those OTUs that were found in at least 30% of the samples in order to ensure that the data were less sparse, as well as being more amenable to calculate correlation coefficients. In addition, as data dimensionality reduction is essential for the limited sample size, we referenced a previous study for preselecting bacteria. Eventually, the relative abundance of 51 OTUs was applied for WGCNA. Module preservation was assessed using the Z-summary as implemented in the WGCNA package of R software. We, in the present study, selected three topological properties (degree, betweenness, and closeness) of the network to extract important OTUs from the network. Accordingly, those OTUs with higher degree, betweenness, and closeness demonstrate that those are potential OTUs associated with GO.

Correlation between GO-Related Traits and Gut
Microbiota. We performed the correlation analysis between International Journal of Endocrinology the modules and GO-related traits, which included TRAb, TPOAb, TGAb, and CAS. In the module-trait correlation analysis, the module eigenOTU was defined as the first principle component of a module, which was used to calculate Pearson's correlation coefficient between a module and a GOrelated trait. Significance of the correlation was determined by an asymptotic P value. For those modules associated with GO-related traits, Pearson's correlation analysis was carried out to determine the association between each of the OTUs involved in the modules and GO-related traits.

Negative Binomial Regression Model for Further Investigation of Relationships between GO-Related Traits and
Gut Microbiota. As the identified bacteria showed a significant association with GO-related traits based on Pearson's correlation analysis, we further investigated the relationships between these bacteria and GO-related traits after adjusting for age and sex. Owing to the skewed and overdispersed distribution of the absolute abundance of bacterial taxa, the negative binomial regression model was used in the analysis. e absolute abundance of each microbiota was taken as a dependent variable, while TRAb, TPOAb, TGAb, and CAS were taken as independent variables, and age and sex were adjusting factors. P < 0.05 was considered statistically significant.

Latent Class Analysis (LCA) to Classify GO Patients.
As the identified bacteria showed significant association with GO-related traits, we further used them to classify GO patients by applying LCA and indicate their contributions to GO-related traits. LCA is a subset of structural equation modeling and is used to find out groups or subtypes of cases in multivariate categorical data, in which these subtypes were previously called "latent classes" [16]. At present, the latent class models are rarely applied in the analysis of microbiome data, in spite of the evolutionary, temporal, and count structure that could be directly incorporated through such models [17]. In the current study, we tried to use LCA to analyze the microbiome data: the absolute abundance of the identified bacteria was divided into three levels (<33rd percentile, >33rd percentile and <66th percentile, and >66th percentile), and these three levels were encoded as 1, 2, and 3.
e LCA can detect the presence of the combination patterns for identified bacteria in latent classes. e general evaluation indexes of LCA are as follows: likelihood ratio (LR), Akaike information criterion (AIC), Bayesian information criterion (BIC), and entropy. e smaller the LR, the better the model fits the data. e smaller values of the AIC and BIC indicate the better fitting of the data. e greater the entropy, the better the model fits the data. Herein, we applied multiple evaluation indexes which outperformed the single index to determine the number of latent classes. To our knowledge, different indexes may have different evaluation results for the model; thus, the interpretation of the model is extremely important. We also compared differences in CAS, TRAb, TPOAb, and TGAb among the latent classes. We attempted to indicate whether gut microbiota contribute to the GO-related traits. e flowchart of our study is shown in Figure 1.

Cooccurrence Network of Gut Microbiota.
In the present research, Pearson's correlation analysis was used to quantify the cooccurrence network of gut microbiota. High correlation reflected the interactions between sources of bacteria and similarities in their responses to environmental conditions. By applying the modularity scores, we found the dense communities in graphs. At the phylum level, the number of positive correlations was 9, whereas the number of negative correlations was 2, and the number of modalities was 3 (Figure 2(a)). As illustrated in Figure 2(a), Fimicutes and Bacteroidetes, as the most predominant phyla in GO patients, are involved in the same module. At the class level, the number of positive correlations was 28, whereas the number of negative correlations was 2, and the number of modalities was 4 ( Figure 2(b)).

WGCNA.
In WGCNA, we set up deepSplit � 2 and minModuleSize � 10 as parameters for the dynamic tree cut algorithm. e soft thresholding of power can be used in network construction; thus, we selected power � 1, which caused the fitted R 2 value and the mean connectivity to be the highest (Supplementary Figure 1).
Four modules, namely, MEgrey, MEturquoise, MEblue, and MEbrown, were eventually selected. Meanwhile, MEgrey was taken as the reference module into account, which included 9 OTUs, and the number of OTUs in MEturquoise, MEblue, and MEbrown was 17, 13, and 12, respectively. MEgrey, MEturquoise, MEblue, and MEbrown were colored grey, turquoise, blue, and brown, respectively. e cluster dendrogram of four identified modules is displayed in Figure 3(a). Figure 3(b) depicts the heatmap of an adjacency matrix of eigenOTU, which was defined as the first principle component of a module. e network heatmap of all OTUs is shown in Figure 3(c).
After the edge weights were filtered according to r > 0.1, a simplified network was attained. e size of the network, the number of network edges, the average degree, the average path length, and the graph density were 38, 94, 4.947, 2.891, and 0.134, respectively. Hence, we extracted the top 10 OTUs with outstanding topological properties ( Table 2). ese OTUs had higher degree, betweenness, and closeness demonstrating that these OTUs link with further OTUs, and thus, potential OTUs were associated with GO. As shown in Table 2, the majority of these OTUs (90%) were contained in the turquoise module (MEturquoise). Module preservation was assessed using the Z-summary, and the results showed that conservation of a turquoise module was confirmed. e simplified network obtained from WGCNA is illustrated in Figure 4. In this network, the isolated nodes were removed. From Table 2 and Figure 4, we can see that the proportion of Bacteroidetes in the top ten OTUs is very high. 4 International Journal of Endocrinology

Correlation between Modules and GO-Related Traits.
Pearson's correlation analysis was performed between the identified modules and GO-related traits, in which the traits were CAS, TRAb, TPOAb, and TGAb. e first principle component of a module was used to calculate Pearson's correlation coefficient between a module and a GO-related trait. We noted that the turquoise module (MEturquoise), which contained the majority of the OTUs with outstanding topological properties, was negatively correlated with TRAb (r � − 0.36, P � 0.04) (Figure 5(a)). As depicted in Figure 5(a), each cell of the matrix contained a correlation coefficient between one OTU module and a GOrelated trait. Gut microbiota data of Chinese patients Figure 1: Flowchart of our work. Firstly, WGCNA was used to generate the network and to identify GO-related network modules based on the OTU relative abundance. ose gut microbiota with outstanding topological properties were extracted. Secondly, the module-trait association analysis was performed, and the Pearson correlation analysis was used to determine the association between each of the OTUs involved in the modules and GO-related traits. irdly, the relationships between the identified bacteria and GO-related traits were further investigated based on the negative binomial regression after adjusting for age and sex. Finally, the identified bacteria were used to classify GO patients by applying the latent class analysis.
(a) (b) Figure 2: Operational taxonomic unit-(OTU-) based cooccurrence network analyses. e cooccurrence network of gut microbiota was constructed based on the correlation analysis. By applying the modularity scores, the dense communities were identified in graphs at the (a) phylum level and (b) class level.

Negative Binomial Regression Model to Find Out the Relationships between GO-Related Traits and Gut Microbiota.
As the five identified bacteria showed a significant association with GO-related traits, we additionally used the negative binomial regression model to investigate the relationships between these bacteria and GO-related traits (TRAb, CAS, and TGAb) after adjusting for age and sex. As presented in Table 3, s_Prevotella_copri, s_Bacteroides_stercoris, and f_Prevotellaceae were all correlated with TRAb and TGAb, while g_Bacteroides was correlated with CAS.

LCA to Classify GO Patients.
In this stage, the five identified bacteria (s_Prevotella_copri (Y1), s_Bacter-oides_stercoris (Y2), f_Prevotellaceae (Y3), g_Bacteroides (Y4, OTU_736), and g_Bacteroides (Y5, OTU_1112)), which showed the significant correlation with GO-related traits, were selected. In LCA, the absolute abundance of the five bacteria was divided into three levels (<33rd percentile, >33rd percentile and <66th percentile, and >66th percentile) and were then coded as 1, 2, and 3. We applied LCA to detect the presence of the bacteria in latent classes. After applying multiple evaluation indexes in LCA (Supplementary Table 2 and Supplementary Figure 2), four latent variables were set; besides, the results of analysis are shown in Table 4 and Figure 6.
We then compared differences in TRAb among the four latent classes. e results showed that the fourth latent class had a higher TRAb value than other three latent classes. e fourth latent class included GO patients with high abundance of s_Prevotella_copri and g_Bacteroides (OTU_736) (Figure 7).

Discussion
In the present study, we applied a metabolic-network-driven analysis to explore the relationships between gut microbiota and GO-related traits and identified a number of novel associations between gut microbiota and serum TRAb. In our previous study, linear discriminant analysis effect size (LEfSe) analysis and random forest analysis showed that Prevotellaceae was one discriminative feature, which could distinguish GO patients from controls obviously [9]. In the present study, Prevotellaceae was identified as an important family of bacteria associated with TRAb. Additionally, we applied LCA to classify GO patients into four different subclasses based on their gut microbiota constitution and found that those GO patients with high abundance of s_Prevotella_copri and g_Bacteroides had a higher TRAb level.
In humans, approximately 30-400 trillion microorganisms colonize the human intestinal tract, and their composition depends on environmental and immunogenetic factors [18]. Alterations in bacterial function and diversity may contribute to the development of autoimmune diseases and infectious diseases [19]. A number of scholars demonstrated that intestinal dendritic cells and macrophages are hyperresponsive to pathogen-associated molecular patterns during equilibrium. When epithelial barrier breakdown occurs, the pattern recognition receptors, which are present in innate immune cells, detect gut microbiota through toll-like receptors, soluble retinoic acid-inducible gene I, NOD-like receptors, or melanoma differentiationassociated protein 5, inducing an inflammatory cascade and activation of adaptive immune responses [20,21]. e variations in gut microbiota composition have been described in patients suffering from autoimmune diseases, including GD and HT [22][23][24]. A previous research demonstrated that the gut is mainly inhabited by two phyla of bacteria in humans: Firmicutes and Bacteroidetes, and the latter was mostly dominated by Bacteroides and Prevotella genera [25]. Studies showed that Bacteroides bacteria play an important role in consuming carbohydrates, while the Prevotella species dominate fibers especially [26,27]. A number of scholars reported that Prevotella is associated with chronic inflammatory conditions, including rheumatoid arthritis (RA) and systemic T-cell activation in HIV-1 infection [28,29]. In addition, s_Prevotella_copri, the most abundant species in Prevotella, has been found to be Negative binomial regressions were performed on gut microbiota after adjusting for age and sex. * P < 0.05; * * P < 0.01. CAS: clinical activity score; TGAb: antithyroglobulin antibody; TRAb: thyrotropin receptor antibody.  International Journal of Endocrinology correlated with the development of rheumatoid arthritis (RA). It has been previously unveiled that patients with RA have differential reactivity of immunoglobulin G (IgG) or IgA immune responses with s_Prevotella_copri. e responses were either IgG antibodies to P. copri, suggestive of a systemic immune response, or IgA antibody responses, suggestive of a mucosal immune response [30]. To date, the correlation of s_Prevotella_copri with thyroid-associated diseases has remained elusive. In the present study, the level of s_Prevotella_copri was positively correlated with a higher serum level of TRAb, which may be related to active orbital inflammation. In the current research, we found that f_Prevotellaceae was correlated with levels of TRAb in GO patients. Next, we detected f_Prevotellaceae DNA in the gut and antibody in the serum in patients [31]. e study showed that patients with inflammatory bowel disease (IBD) may exhibit a concomitant increase in Bacteroidetes [32]. In the present study, patients with higher serum levels of TRAb had high abundance of s_Prevotella_copri and g_Bacteroides in LCA. As expressed previously, TRAb is an independent risk factor for GO and can predict severity and outcome of the disease [33]. e present research revealed the relationship  between bacteria and TRAb, and our findings may be helpful for the prediction of GO by intestinal microorganisms in the future study. e relationship between s_Prevotella_copri and g_Bacteroides in regulating TRAb-associated immune responses needs to be further explored and elucidated. e gut microbiome is a complex and metabolically active community that directly influences host phenotypes [34]. In fact, the structure of the gut microbiome is influenced by several factors, including interactions between its members. erefore, it is highly essential to understand the combination of microbial effects underlying their contributions to disease [35]. A general approach to infer interactions between bacteria in the gut microbiota is to quantify the cooccurrence of OTUs. To address this issue, network biology is an emerging field that represents biology as networks that capture the relations between the parts of a complex biological system, such as molecules, processes, organs, or even different organisms. Meanwhile, WGCNA, which is used to identify modules in gene coexpression networks, can facilitate identification of communities within microbial cooccurrence networks [34,36]. For instance, a study showed that certain metabolites strongly correlate with the microbial community structure, and some of these correlations are specific for the prediabetic state [34]. In the present study, we assessed the association between gut microbiota and GO-related traits in the combination of microbial effects.
e present study has a number of limitations: First, the limited sample size influenced the achieved results. To address this issue, we performed the data dimensionality reduction. In WGCNA, we kept those OTUs found in at least 30% of the samples in order to ensure that the data were less sparse, as well as being more amenable to calculate correlation coefficients; besides, we referenced a previous study for preselecting bacteria. In addition, we used module preservation analysis to estimate the robustness of the identified modules and analyzed the topological properties to identify key bacteria. Specially, we applied the correlation analysis, negative binomial regression, and LCA to further validate the identified bacteria. Although these layer-by-layer analyses achieved the goal of the data dimensionality reduction and avoided the identified bacteria to be false-positive, an enlarged sample size will be required in the future study to validate the potential relationship between these bacteria and GO. Second, whether the changes in the fecal microbiota might be a cause or a consequence of GO development should be further explored. e analysis of 16S rRNA gene fragments did not provide data related to the functional traits of the bacterial genera being present in the communities. In addition, Pearson and Spearman's correlations were limited to detect the cooccurrence of bacteria; thus, further ensemble approaches, such as SparCC, Kullback-Leibler divergence, and Bray-Curtis dissimilarity, will be used in the next study.
ird, our findings may be related to thyroid autoimmunity rather than GO; thus, further analyses about the gut bacteria between GO patients and patients with GD without GO should be performed. Finally, TRAbs were subdivided into thyroid-stimulating antibody (TSAb), thyroid-blocking antibody (TBAb), and neutral antibody (neutral Ab) according to their functional effects [10]. e assay did not distinguish between stimulating, blocking, and no functional effects on the thyroid gland.
In summary, our study suggested the potential relationships between the composition of the gut microbiota and TRAb in GO patients. It might provide a new insight into the interaction between the intestinal microbiota and TRAb-associated immune responses in GO patients.

Data Availability
e data that support the findings of this study are available from the corresponding author upon reasonable request.

Ethical Approval
is study was conducted with the approval from the Ethics Committee of Beijing Tongren Hospital, Capital Medical University. All procedures performed in this study were in accordance with the 1964 Helsinki Declaration and its later amendments.

Consent
Informed written consent was obtained from all individual participants included in this study.

Conflicts of Interest
e authors declare no conflicts of interest.

Authors' Contributions
Ting-Ting Shi and Lin Hua contributed equally to this paper.