Alteration in gut microbiota is associated with dysregulation of cytokines and glucocorticoid therapy in systemic lupus erythematosus

ABSTRACT

A growing corpus of evidence implicates the involvement of the commensal microbiota and immune cytokines in the initiation and progression of systemic lupus erythematosus (SLE). Glucocorticoids have been widely used in the treatment of SLE patients, however, glucocorticoid treatment carries a higher risk of other diseases. Using the 16S rRNA technique, we investigated the differences between the gut microbiota associated with the immune cytokines of SLE and relevant glucocorticoid treatment in a female cohort of 20 healthy control subjects (HC), 17 subjects with SLE (SLE-G), and 20 SLE patients having undergone glucocorticoid treatment (SLE+G). We observed that the diversity and structure of the microbial community in SLE+G patients were significantly changed compared to that of SLE-G patients, whereas the gut microbial community of the SLE+G group showed a similarity with the HC group, which implicate that the shift in the gut microbiome could represent a return to homeostasis. Furthermore, the up-regulations of immune cytokines in SLE-G were identified as closely related to gut dysbiosis, which indicates that the overrepresented genera in SLE patients may play roles in regulating expression level of these immune cytokines. This associated analysis of gut microbiota, glucocorticoid therapy, and immune factors might provide novel and insightful clues revealing the pathogenesis of SLE patients.

KEYWORDS:

• Systemic lupus erythematosus (SLE)
• gut microbiota
• cytokines
• glucocorticoid
• autoimmune disease

Introduction

Systemic lupus erythematosus (SLE) is a prototypical systemic autoimmune disease with a significant disease burden across the world among different ethnic, racial, and age groups.¹ Its clinical manifestations are characterized by a higher frequency in women and combinations of symptoms, including skin rashes, glomerulonephritis, neurologic disorders, and severe vasculitis, which suggests that SLE may be a complex rather than a single pathogenic process.² However, the etiology and pathogenesis of SLE remain largely unknown.

The human immune system, which is supposed to defend the body against threats, can attack its own healthy tissues when facing certain disorders, leading to serious autoimmune diseases.³ An autoimmune disease is a kind of chronic systemic disease characterized by the excessive activation of auto-reactive T cells and B cells, which induces a great quantity of auto-antibodies and leads to extensive damage to multiple systems.⁴

Increasing studies have provided strong evidences that gut microbiota could stimulate the production of a variety of inflammatory cells and inflammatory cytokines through specific molecular patterns which are associated with autoimmune diseases.^3,5,6 Recent findings show that abnormalities in the gut microbiota can be detected in patients with SLE.⁷ Active SLE patients had significant dysbiosis in gut microbiota with reduced bacterial diversity and biased community constitutions.^8,9 Furthermore, Doua Azzouz et al. reported that SLE patients were found to have characteristic patterns of gut microbiome dysbiosis that directly paralleled disease activity, and suggested that gut commensals may contribute to the immune pathogenesis of lupus nephritis.¹⁰ Nevertheless, the detailed association among SLE, the immune system, and gut microbiota is still far from clear. Glucocorticoids are widely used as anti-inflammatory and immunosuppressive agents in many immune-mediated gastrointestinal diseases.¹¹ Recent study has shown that gut microbiota can affect the efficacy of glucocorticoids.¹² Moreover, glucocorticoid treatment carries a higher risk of opportunistic infections, iatrogenic osteoporosis, and avascular necrosis, as well as an increase in the risk of cardiovascular events and of adverse psychiatric effects, such as psychosis and manic episodes.¹³ Hence, it is important to find a new way to mitigate the side effects of glucocorticoid therapy.

In this study, we first aimed to investigate the difference in gut microbiota composition among SLE patients without glucocorticoid therapy compared, respectively, to healthy controls and SLE patients who had undergone glucocorticoid therapy. Second, to determine the significant genera associated with deregulated immune cytokines, we then detected the levels of 12 cytokines in three groups and performed the associated analysis for the differential abundance bacteria and the differential expressed immune cytokines between the two groups, respectively. Finally, we identified the specific genus correlative to the deregulation of cytokines leading to immune system disorder. This associated analysis of gut microbiota, glucocorticoid therapy, and immune factors might provide novel and insightful clues revealing the pathogenesis of SLE patients.

Results

Clinical characteristics of the recruited subjects

Seventeen of SLE patients without glucocorticoid therapy (SLE-G), 20 of SLE patients with glucocorticoid therapy (SLE+G) and 20 healthy controls (HC) were recruited in this study. All patients and controls are female. There were no significant differences in age, and body mass index (BMI) among the groups, whereas there were statistically significant differences in erythrocyte sedimentation rate (ESR) among the groups (Table S1). The detailed clinical characteristics including the dose, duration of glucocorticoid treatment of all subjects, and SLE disease activity index (SLEDAI) are shown in Table S2.

SLE-G patients displayed alterations of gut microbiota compared to HC and SLE+G, respectively

The alpha diversity of bacterial communities was evaluated according to Observed OTUs, Pielou’s evenness, Faith’s phylogenetic diversity, and Shannon’s diversity index (Figure 1a to c and Figure S1a). Generally, the Observed OTUs of HC and SLE-G showed a significant difference (P < .05), as did SLE-G and SLE+G (Figure 1a). However, HC and SLE+G showed no statistical difference (Figure 1a, P > .05). Similarly, we found a significant difference in Pielou’s evenness (Figure 1b) and Faith’s phylogenetic diversity (Figure 1c) between HC and SLE-G, but no significant difference in the Shannon index was shown between HC and SLE-G (Figure S1a; P > .05).

Figure 1. SLE patients displayed alternation of gut microbiota and glucocorticoid treatment changed gut microbiota composition of SLE patients.

(a), (b) and (c). Alpha diversity in SLE patients and healthy controls. Box plots depicted greater gut microbial diversity in SLE-G compared with other two groups according to the Observed OTUs, Pielou’s Evenness and Faith’s Phylogenetic Diversity indexes. (d). Phylogenetic distances between samples from three groups were determined by Unweighted UniFrac PCoA (Principal co-ordinates analysis) of the overall gut microbiota. (e). Venn diagram of the shared species number at the level of phylum among HC, SLE-G and SLE+G groups. (f). Venn diagram of the shared species number at the level of genus among HC, SLE-G and SLE+G groups. (g). The ratio of Firmicutes/Bacteroidetes of HC, SLE-G and SLE+G group (“p” means phylum). (h). ROC analysis of sensitivity and specificity of the model for prediction of HC and SLE-G patients. (i). ROC analysis of sensitivity and specificity of the model for prediction of SLE-G and SLE+G patients. (j). The top 30 most differential abundance genera across groups at the criteria of P value < .05 by Kruskal–Wallis test, the P value corrected by the Benjamini–Hochberg FDR (“g” means genus). (k). Cladogram using the LDA model results for the bacterial hierarchy. Differences were represented by the color of the most abundant class (red, SLE-G; green, SLE+G; yellow, insignificant). The diameter of each circle is proportional to the taxon’s abundance. Each ring represents the next lower taxonomic level. (l). LEfSe analysis indicated genera strikingly different in gut microbiota of SLE patients with/without hormone therapy.

Beta diversity was measured via Unweight UniFrac distances, using the first two principal coordinates to visualize the dissimilarity distances and variation between samples. The PCoA scatterplot revealed the clear clustering of gut bacterial communities between SLE-G and other two groups (Figure 1d). A one-way analysis of similarity (ANOSIM) also revealed that the significant separation among the three groups (R = 0.226, P = .001), while HC and SLE+G showed similar community (R = −0.004, P = .488). Furthermore, we compared the number of bacteria at the phylum and genus levels. Distinct differences in bacterial composition were observed in the gut microbiota of all subjects. At the phylum level, the results showed 11 common phyla among three groups, and SLE-G displayed 14 specific phyla (Figure 1e). At the genus level, the analysis also showed 64 overlapping genera among three groups, and SLE-G displayed 101 specific genera (Figure 1f). In short, a Venn diagram revealed that the gut bacteria of SLE-G were different from those of HC and SLE+G. At the phylum level, compared to HC and SLE+G group, the gut microbiota of SLE-G patients was characterized by an increased Bacteroidetes level (Figure S1b) and a lower ratio of Firmicutes/Bacteroidetes (Figure 1g). Interestingly, glucocorticoid therapy may relate to restoring the ratio of Firmicutes/Bacteroidetes by reducing the level of Bacteroidetes.

Next, using ROC analysis of R package (version 3.1.4, pROC package), we tested whether the ratio of Firmicutes to Bacteroidetes can be used as a biomarker to distinguish HC, SLE-G and SLE+G group, respectively. In result, ROC analysis revealed that 88.8% AUC between HC and SLE-G (Figure 1h), 88.5% AUC between SLE-G and SLE+G (Figure 1i), and 67.2% AUC between HC and SLE+G (Figure S1c), which implied that the ratio of Firmicutes to Bacteroidetes may define a biomarker to diagnose SLE-G compared to HC and SLE+G, respectively.

To compare the gut microbiota composition of the three groups, we conducted a statistical analysis of gut microbiota at the genus level using the Kruskal–Wallis test. As shown in Figure 1j, the top 30 genera with statistically significant differences among the three groups (P < .05, FDR<0.05) were selected to display. Meanwhile, HC and SLE+G shared some common genera, such as Gemmiger, Lactococcus, Bifidobacterium, Streptococcus, and Desulfovibrio, which were enriched both in HC and SLE+G individuals, but decreased in the SLE-G group. Interestingly, compared with other groups, Akkermansia and Lactobacillus had the highest abundance in SLE+G group (Figure 1j).

We also analyzed the relationship among the top 30 genera with a statistical difference between the two groups. Bacteroides, Coprococcus, Prevotella, Parabacteroides, Succinivibrio, and Bilophila in SLE-G patients were increased compared to those of the HC group, and Gemmiger, Dialister were relatively reduced (Figure S1d). These genera, such as Bacteroides, Oscillospira, Phascolarctobacterium, Succinivibrio, Prevotella, and Bilophila were elevated in the SLE-G compared to the SLE+G group (Figure S1e). In fact, there was a relatively similar abundance of bacteria in the HC and SLE+G groups (Figure S1f), suggesting that glucocorticoid therapy might change the ecological structure of gut microbiota in the SLE patients by unknown mechanisms.

Next, to identify the specific genera associated with SLE and specific genus biomarkers for SLE diagnosis, a cladogram was used to describe the structure of the gut microbiota and the predominant bacteria in the HC, SLE-G, and SLE+G groups (Figure 1k). The LEfSe revealed that 10 taxa were differentially abundant in SLE-G and SLE+G groups. In these taxa, four were shown in red color, six were shown in green color were identified as enriched within the SLE-G and SLE+G, respectively. Moreover, compared with SLE+G, abundance of Bifidobacterium, Streptococcus, and Enterobacteriales was lower, Bacteroidetes and Bilophila were enriched in the SLE-G samples (Figure 1l).

In short, the above analysis indicates that SLE patients have an increased number of bacteria at both the phylum and genus levels. Glucocorticoid treatment might reduce the bacterial diversity in SLE patients through unclear molecular patterns.

To check whether the SLEDAI (systemic lupus erythematosus disease activity index) is associated with change of gut microbiota, we divided all of SLE-G patients into SLEDAI<6 and SLEDAI≥6 group. As shown in Figure S3a to c, no differences in alpha diversity were observed between two groups. Similar to the alpha diversity, no statistical difference in microbiota community structure was observed between SLEDAI<6 and SLEDAI≥6 group (Figure S3d and e). The abundance of top 10 genera between the two groups was calculated, but no statistically significant genera were found (Figure S3f).

Aberrant network of gut microbiota occurred in SLE patients

To compare the gut microbiota ecology of the three groups, we further assessed the interactions among the genera in three groups. The control samples and SLE+G patients mainly featured two networks with scattered genera from both Firmicutes and Proteobacteria phyla (Figure 2a and c). The microbial community of SLE-G obviously featured a more complex network, including 76 genera from the Proteobacteria phylum compared to both the HC and SLE+G groups (Figure 2b). The SLE+G patients displayed a similar network as the HC patients, showing a strong positive correlation between the genera (Figure 2c). However, the correlation of the genera in SLE-G was much different from those of HC and SLE+G. To quantify such a difference, we counted the number of edges (connections) and the centrality of nodes (genera) in the three microbial networks. Despite having a few overlapped edges, 802 of the edges and 801 of the edges were specific to the SLE-G, respectively, compared to the HC and SLE+G groups (Figure S2a and S2c). The closeness and eigenvectors of the shared genera in SLE-G were also quite different from those of HC and SLE+G (Figure 2d). Impressively, although SLE+G and HC displayed many common genera, only one edge overlapped in both networks; moreover, 37 and 26 edges were specific to HC and SLE+G, respectively (Figure S2b). Taken together, the network of SLE patients displayed many more connections than those of HC and SLE+G patients due to the increase in genera diversity in SLE patients.

Figure 2. Co-occurrence network of microbiota in HC (a), SLE-G patients (b) and SLE+G patients (c).

The circle size represents the relative abundance, and the density of lines between circles represents the Spearman coefficient. Green links stand for positive interactions between nodes, and red links stand for negative interactions (“g” means genus; “p” means phylum). (d). The heatmap for the discrepancies of counted nodes in three co-occurrence networks.

SLE patients showed deregulation of cytokines

Cytokines serve as key mediators of inflammation and tissue damage in a variety of immune-mediated disorders.¹⁴ It is well known that commensal microbiota modulate cytokine responses.¹⁵ Next, to reveal the association between gut microbiota and immune cytokines in SLE patients, we further analyzed the expression levels of inflammatory cytokines in the blood of SLE-G compared to HC and SLE+G patients (Figure 3). In total, twelve cytokines displayed higher expression levels in the SLE-G and SLE+G groups compared to HC, respectively. We observed that seven cytokines, including IL-1β, IL-6, IL-8, IL-17, IL-21, IL-22 and TWEAK showed no significant changes between the SLE-G and SLE+G groups. However, the expressions of five cytokines, including IFN-γ, IL-2R, IL-10, IL-35 and TNF-α were significantly different in the HC, SLE-G, and SLE+G groups. The expressions of these cytokines were highest in SLE-G patients and lowest in HC, while a slight decrease in SLE+G further supported previous studies that glucocorticoid therapy restored immune-mediated disorders.¹⁴

Figure 3. Analysis of twelve immune cytokines in HC, SLE-G and SLE+G groups.

All of subjects’ serum was collected and run on ELISA. The differences are calculated by T test (*p  < .05, **P < .01, ***P < .005).

The differential genera between any two groups showed a correlation with differentially expressed cytokines

To find the potential bacteria associated with the immune disorder in SLE, we conducted a correlation analysis of the different genera and the different cytokines between any two groups. As shown in Figure 4a, Bacteroides, Succinivibrio, Bilophila, and Parabateroides, which were significantly increased genera in SLE-G relative to HC, were positively correlative with IL-17, TWEAK, IL-2R, IL-21, IL-35, IL-10, and IFN-γ, which implicated these bacterial genera as core bacteria which possibly play essential roles in stimulating an increase in these immune factors. Conversely, Dialister and Gemmiger, which were significantly reduced in SLE-G, showed a strong negative association with immune factors IL-17, IL-2R, and IL-35 which suggested that these genera have the potential to reduce inflammatory cytokines in SLE patients. Nevertheless, there was relatively little significant positive correlation between the different genera and immune factors in HC and SLE+G. But impressively, the elevated Bifidobacterium abundance in SLE+G relative to HC showed positive correlation with IL-17, IL-22 and IL-6 (Figure 4b).

Figure 4. Correlation analysis of these differential genera and cytokines between any two groups.

The abundance of differential species (top 30) enriched in HC and SLE-G (a), in HC and SLE+G (b), and in SLE-G and SLE+G (c) was analyzed for conversation with differential immune factors using Spearman’s correlation analysis. The correlation effect is indicated by a color gradient from green (negative correlation) to purple (positive correlation). Correlation coefficients and p-values (*p < .05, **P < .01, ***P < .001) are shown.

Similar to the correlation analysis of the differential genera between HC and SLE-G, Succinivibrio, Bilophila, Oscillospira, and Phascolarctobacterium which were significantly elevated in SLE-G relative to SLE+G, were positively correlated with IL-35, IL-2R, IL-10, and IFN-γ (Figure 4c).

Functional alteration for the differential abundance of genera between any two groups

To predict whether glucocorticoid therapy induced microbial changes could modulate the metabolic function of gut microbiota, we conducted a functional metagenomics prediction based on 16S rRNA sequencing using the PICRUSt.¹⁶ The PCA of STAMP analysis based on KEGG module prediction indicated that SLE-G had a more scattered distribution was different from that of HC and SLE+G (Kruskal–Wallis H-test, P < .05; Multiple testing adjustment, FDR<0.05, Figure 5a). The first principal component (PC1) could explain 73.3% of the functional changes caused by glucocorticoid therapy. As shown in Figure 5b, the KEGG pathways in the SLE-G group were different from those of the SLE+G and HC groups. Thetwenty-seven modules elevated in SLE-G were involved in Citrate cycle (TCA cycle), Glycolysis/Gluconeogenesis, Fructose and mannose metabolism, Galactose metabolism, Other glycan degradation, MAPK signaling pathway, Glycosaminoglycan degradation, and Glycerophospholipid metabolism, and Pentose phosphate pathway. The overrepresentation of the polysaccharide degradation pathway in SLE-G may be related to the high abundance of Bacteroides, which has been shown to have a wide capacity for degrading polysaccharides in those samples.^8,17 We observed eight modules which were elevated in SLE+G, including ABC transporters, Amino acid metabolism, Bacterial motility proteins, Bacterial secretion system, Fatty acid metabolism, Lipid metabolism, Glycan biosynthesis and metabolism, and Phosphotransferase system (PTS). In all, the predication analysis indicated that gut dysbiosis may induce a disease-linked state through the interference of physiological metabolic functions.

Figure 5. Microbial functions annotation in HC, SLE-G and SLE+G patients.

(a). PCA analysis based on the relative abundance of KEGG function modules in fifty-seven samples of three groups. KEGG function modules between two groups (Kruskal–Wallis H-test, P < .05; Multiple testing adjustment, FDR<0.05) were considered as significantly different. (b). The average abundance of 35 KEGG modules differentially enriched in HC, SLE-G, and SLE+G. Twenty-seven modules enriched in SLE-G, and eight modules overrepresented in SLE+G. (c). The relationship between thirty-five differential modules and the 30 top altered genera (Figure 1j) in healthy controls and SLE patients is estimated by Spearman’s correlation analysis. Green color means negative correlation, on the contrary, purple is positive correlation. Correlation coefficients and p-values (*p < .05, **P < .001, ***P < .0001) are shown.

To investigate which genera are correlative with specific KEGG pathways, we predicated the relationships of KEGG pathways associated with the 30 significantly different genera presented in Figure 1j. Impressively, as shown in Figure 5c, Bacteroides, Parabacteroides, Alistipes, Fusobacterium and Adlercreutzia which were over-represented in the SLE-G group showed positive correlations with Other glycan degradation, Glycosaminoglycan degradation, Glycosphingolipid biosynthesis, MAPK signaling pathway, Protein digestion and absorption and Adipocytokine signaling pathway, which implicated these four genera may involve in pathogenesis of SLE through regulating these correlative signaling pathways.

Discussion

In summary, as shown in Figure 6, we identified that the relative abundances of the phylum Bacteroidetes, Firmicutes, and Proteobacteria were significantly increased in SLE-G compared to HC and SLE+G, while a decrease in the phylum Verrucomicrobia, Firmicutes, and Proteobacteria was shown in SLE-G. The increasing genera, such as Bacteroides, Bilophila, Parabacteroides, and Succinivibrio in SLE groups were positively correlated with inflammatory cytokines IL-17, IL-21, IL-2R, TWEAK, IL-35, IFN-γ, and IL-10. Meanwhile, the decreasing genera, such as Dialister and Gemmiger, were negatively correlated with IL-17, IL-2R, and IL-35. Glucocorticoid therapy may inhibit the above changes in cytokines while modulating the composition of gut microbiota. Furthermore, we speculated that the deregulation of immune cytokines might accelerate the disorder of the host cell-signaling pathway, such as the mitogen-activated protein kinase (MAPK) signaling pathway, glycosaminoglycan degradation, other glycan degradation, and so on, leading to the pathogenesis of SLE disease. Nevertheless, in what way dose glucocorticoid regulate the community of gut microbiota need to be further exploration. To our knowledge, this study is the first research evaluating the effect of glucocorticoid therapy for SLE patients on the gut microbiota composition correlative with cytokine changes.

Figure 6. The summary for the analysis of gut microbiota and cytokines in SLE patients with/without Glucocorticoid treatment.

Dashed arrows and Question Mark indicated this interaction need to be investigated.

We observed that the microbial community diversity and structure in SLE+G patients were significantly changed compared to those of SLE-G patients, whereas the gut microbial community of SLE+G patients was similar to that of the HC group. Consistent with previous studies,⁸ we also observed a reduction in the ratio of Firmicutes to Bacteroidetes in the feces of SLE patients compared to healthy controls. Impressively, SLE patients undergoing glucocorticoid therapy had an increase in the ratio of Firmicutes to Bacteroidetes. The result suggested the ratio of Firmicutes to Bacteroidetes may be used as biomarkers to evaluate the efficacy of glucocorticoid therapy in SLE patients. Interestingly, we found that a group of core bacteria genera, including Lactococcus, Streptococcus, and Bifidobacterium, were enriched in the HC and SLE+G groups, but reduced in the SLE-G group. However, Azzouz et al. found expansions of R. gnavus primarily in patients with active Lupus nephritis that were not prevalent in the current report.¹⁰ Recent studies have indicated that geographical location of the hosts’ residence showed the strongest associations with microbiota variations,¹⁸ and many bacterial species including Ruminococcus are unevenly distributed across human individuals in the same geographic area and between populations in different locations which might imply that the gut microbiota of SLE patients in different regions may have different specific characteristics. However, a recent research suggests that R. gnavus strains in China may have some of the same conserved strain features as those in the US,¹⁹ and Liu et al reported some species of Ruminococcus were reassigned include in the genus Blautia,²⁰ due to 16S rRNA sequencing technique limitation in species identification, future study focusing on species by using matagenomics technique will assistant us to identify specific Ruminococcus species.

Bacteria inhabiting the human gut are critical for digestion of the plant-derived glycans that compose dietary fiber.²¹ In agreement with these results, a function prediction based on gut microbiota composition showed that SLE displayed an overrepresentation of activity-related glycan metabolism, while SLE+G showed a decreased activity. Moreover, genera Bacteroides were overrepresented in SLE, which showed a high correlation with the glycan metabolism, and were largely reduced after glucocorticoid treatment. Bacteroides have been experimentally demonstrated to show extensive glycan degradability.¹⁷ Taken together, these results suggested that glucocorticoid treatment may play a pivotal role in targeting glycan metabolism through reducing the abundance of some bacteria involved in glycan activity.

Many studies have shown abnormal cytokine production in the T cells of SLE patients.²² Consistent with these studies, we found that the expression levels of 12 cytokines were relatively high in the SLE-G group compared to the healthy control group. Interestingly, six cytokines showed a decreased expression in SLE+G compared to the SLE-G group, which suggests that glucocorticoids may play a potential role in reducing the immune cytokines of SLE patients. The associated analysis of different bacteria in the three groups suggested that immune responses were highly related to the variations in genera. SLE+G patients showed elevated levels of Lactobacillus and Bifidobacterium, which are among the most abundant and important symbiotic bacteria in the human gut microbiota and which produce short-chain fatty acids (SCFAs) via the fermentation of dietary fiber.^23,24 The anti-inflammatory effects of SCFAs have been well characterized at both the epithelial and immune cell levels.²⁵ Moreover, SCFAs also influence the intestinal barrier function through the stimulation of tight junctions and mucous production,²⁶ as well as promoting the differentiation of regulatory T cells in the colon, spleen, and lymph nodes, thereby inhibiting inflammation.²² Recently, a number of studies have revealed that the abnormal function of Tregs is relevant to the pathogenesis of SLE;²² however, the immunotherapeutic use of Tregs in SLE is still in the beginning stages.²⁷ Bifidobacterium is an important probiotic necessary for intestinal microbial homeostasis, gut barrier, and lipopolysaccharide (LPS) reduction.²⁸ It has been reported that Bifidobacterium can prevent the excessive activation of CD4⁺ lymphocytes, thus maintaining the balance of Treg/Th17/Th1 in SLE patients.²⁹ Based on the above information, we suggested that the increased level of bacteria producing SCFA in SLE patients, such as Bifidobacterium, could be a potential probiotic of treating autoimmune diseases through the therapeutic modulation of Tregs.

Taken together, our analysis indicates that glucocorticoid therapy has the potential ability to stabilize the intestinal microbe community of SLE patients through some molecular pattern, further decreasing the production of cytokines. The analysis also suggested that a group of specific bacteria could be used to assess the effect of glucocorticoid therapy. Many studies describe dysbiosis as a change in the microbiota that accompanies autoimmune illnesses, but little is known about whether these changes are a cause or consequence of an altered immune state.³⁰ Therefore, further animal experiments are needed to understand exactly how the candidate bacteria influence immune dysfunction in SLE patients.

Glucocorticoid toxicity is well known³¹ to increase cardiovascular events, ophthalmic conditions including cataracts and glaucoma, and adverse psychiatric effects, such as psychosis and manic episodes.¹³ Thus, it remains an urgent challenge to create effective, life-prolonging strategies in the clinical treatment of SLE patients. Although our data indicated that glucocorticoid therapy appears to modulate the abundance of some genera close to the level of the HC group, we also observed that glucocorticoid therapy may modify the microecology of gut, as SLE+G and HC stil showed differences in aboudance of some genera. Inappropriate modification of gut microbiota composition may lead to an increased risk of side effects from glucocorticoids treatment. Notably, the limitations of our current study are that the samples from same patients before Glucocorticoid treatment were not collected, future study focusing on investigating the gut microbiota shift from the same patients before and after treatment would be more impactful. Recently, Zagarra-Ruiz et al. reported that the dietary resistant starch exerts beneficial effects in lupus-prone hosts through suppressing a pathobiont that promotes interferon pathways implicated in the pathogenesis of human autoimmunity.³² Taken together, our result further suggests that the interactions of diet, microbiota, and host are involved in the development of SLE. The associated analysis of gut microbiota, glucocorticoid therapy, and immune factors in our study might provide novel and insightful clues revealing the pathogenesis of SLE patients.

Material and methods

Ethics statement

The study was approved by the Medical Ethics Committee of The Second Affiliated Hospital of Xi’an Jiaotong University (Ethics NO, 2016055). All participants, or their legally authorized representatives, provided written informed consent upon enrollment.

Study population recruitment and fecal collection

Seventeen female patients without glucocorticoid therapy (SLE-G), 20 female patients with glucocorticoid therapy (SLE+G) and 20 healthy female controls (HC) were consecutively recruited from March 2017 to August 2017 in the Second Affiliated Hospital of Xi’an Jiaotong University. All patients were eligible according to the Systemic Lupus International Collaborating Clinics Classification Criteria for SLE (SLICC),³³ and all enrolled SLE-G patients were further classified according to SLE disease activity index (SLEDAI).³⁴ All patients who had taken antibiotics, microbial agents, or underwent intestinal surgery within 1 month before admission were excluded. SLE-G group included 13 new patients, who were not given any medication before treatment and 4 patients were treated with tripterygium wilfordii or hydroxychloroquine before treatment. When SLE+G enrolled, they were only given glucocorticoid combined with gastric mucosa protection, calcium supplementation and other complementary measures, but they were not given immunosuppressive therapy. And the stool samples of SLE+G were collected after 2 months of glucocorticoid treatment. The HC group had no gastrointestinal disease and had not taken antibiotics for three months prior to the start of the study. In addition, these three groups had no significant differences in age, body mass index (BMI), and other related factors, whereas there were statistically significant differences in erythrocyte sedimentation rate (ESR) among the groups (Table S1). The clinical characteristics of patients, including body mass index (BMI), erythrocyte sedimentation rate (ESR), SLE disease activity index (SLEDAI), disease duration, and patients’ nephritis information are shown in Table S2. The stool and blood samples of each participant were collected in the hospital and were immediately stored at −80°C until analysis.

DNA isolation and 16S rRNA gene sequencing

In total, 100-mg stool samples were used to extract total bacteria genome DNA following the protocol of the DNA extraction kit (#DP328, Tiangen Company, Beijing, China). The concentration and purity of the extracted bacterial DNA were detected using the Qubit 2.0 Fluorometer (Thermo Scientific, USA). The V4 regions of 16S rRNA genes were amplified using the specific primer with the barcode. The 16S rRNA V4 specific primers are 515F GTGCCAGCMGCCGCGGTAA and 806R GGACTACHVGGGTWTCTAAT. Polymerase chain reaction (PCR) products were purified using the GeneJET Gel Extraction Kit (Thermo Scientific). Sequencing libraries were generated using the Illumina TruSeq DNA PCR-Free Library Preparation Kit (Illumina, USA) following manufacturer recommendations, and index codes were added. Sequencing was performed in the Illumina Hiseq platform (Novogene, China). All raw data have been deposited in GenBank. BioProject accession number is PRJNA626512.

16S rRNA-based gut community analysis

Paired-end fastq files were generated by Illumina Hiseq. The analysis was performed by the QIIME 2 (Quantitative Insights Into Microbial Ecology, https://docs.qiime2.org/2018.2/) software pipeline.Firstly, we demultiplex sequence by running the demux emp-paired commands. After demultiplexing reads, we performed quality control using DADA2, denoise the read using dada2 denoise-paired commands. The result of this step generated a feature table that contains counts (frequencies) of each unique sequence in each sample of the dataset, and a sequences data which maps feature identifiers to the sequences which they represent. Then, diversity analysis were available through the q2-diversity plugin. Alpha diversity was evaluated by four indices using QIIME 2, including Shannon’s diversity index, Observed operational taxonomic units (OTUs), Faith’s phylogenetic diversity, and Pielou’s evenness. Beta diversity was measured by calculating the Unweighted UniFrac distance. Then, beta diversity was visualized in PCoA plots by R package (GUniFrac, ape, and ggplot2 package). Comparisons using distance matrices were evaluated by the ANOSIM analysis. The differential abundance genera among groups were examined by the Kruskal–Wallis test, the comparison of any two groups was conducted by the Wilcoxon rank-sum test,and the P value was corrected by the Benjamini–Hochberg FDR procedure for multiple comparisons. To identify biomarkers between healthy individuals and SLE patients at the genus level, taxa summaries were reformatted and inputted into the linear discriminant analysis effect size (LEfSe) via the Huttenhower Lab Galaxy Server.The LEfSe analysis used the non-parametric factorial Kruskal–Wallis (KW) sum-rank test (α = 0.05) to detect features with significant differential abundances (using one-against-all comparisons), and the threshold of LDA score was 2.0.

Co-occurrence network analysis

For the co-occurrence network analysis, the bacterial correlations in the HC, SLE-G, and SLE+G samples were computed, respectively, by Spearman’s correlation coefficient analysis. The significantly correlated genera (FDR<0.05) were visualized by Cytoscape (v3.6.1). The similarities between any two network structures were measured by node closeness and shared correlations.³⁸ Closeness of the nodes was calculated by Cytoscape (v3.6.1) to measure node centralities in each network. The concordance of bacterial correlations in any two groups was inferred by counting the same bacterial correlations across sample types. The results were visualized by the R (version 3.4.3) VennDiagram and gplots package.

Immune cytokines assay for healthy controls and SLE patients

We collected all the participants’ blood serum and tested 12 immune cytokines, including IL-1ß, IL-6, IL-10, IL-35, TWEAK, IL-2R, IL-17, TNF-α, IL-8, IL-21, IFN-γ and IL-22 following the company’s kit procedures. Besides IL-35 (Eiaab Company ELISA kit) and IL-21 (ThermoFisher Company), all other assay kits are from R&D System. The immune cytokines of the participants in the three groups were statistically analyzed GraphPad Prism software (T-test), respectively.

Correlation analysis of gut microbiota and immune cytokines

To determine the association between gut microbiota and immune factors between any two groups, we constructed a correlation analysis according to Spearman’s correlation in R (version 3.4.3, Hmisc package).

Functional predictions for 16S rRNA sequencing data

The functional prediction was performed using the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt)³⁹ for the Kyoto Encyclopedia of Genes and Genomes (KEGG) biomarker prediction implemented in Galaxy (http://huttenhower.sph.harvard.edu/galaxy/).

Author contributions

XL and YX conceived and designed project. MG, HW, ZX, KC, QL and ZD did experiments; MG, HW, XL and ZL performed 16S rRNA analysis. All authors did data analyses and interpretations; XL, MG, HW, and YX prepared and finished the manuscript.

Disclosure of potential conflicts of interest

Authors declare no competing interests.

Acknowledgments

This work was supported by NSFC grant 81671983 and 81871628, Natural science funding BK20161572 from Jiangsu province and starting package from NJMU (Xingyin Liu). Young scientist funding from Jiangsu province (Jianwei Wang and Qisha Liu). NSFC grant 81472876 and 81630081 (Yumin Xia). Starting funding for the team of gut microbiota research in NJMU (Xingyin Liu, Faming Zhang, Yankai Xia and Chuan Su).

Supplementary materials

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [81472876 and 81630081]; National Natural Science Foundation of China [81671983 and 81871628].