Intestinal dysbiosis associated with systemic lupus erythematosus.

UNLABELLED: Systemic lupus erythematosus (SLE) is the prototypical systemic autoimmune disease in humans and is characterized by the presence of hyperactive immune cells and aberrant antibody responses to nuclear and cytoplasmic antigens, including characteristic anti-double-stranded DNA antibodies. We performed a cross-sectional study in order to determine if an SLE-associated gut dysbiosis exists in patients without active disease. A group of 20 SLE patients in remission, for which there was strict inclusion and exclusion criteria, was recruited, and we used an optimized Ion Torrent 16S rRNA gene-based analysis protocol to decipher the fecal microbial profiles of these patients and compare them with those of 20 age- and sex-matched healthy control subjects. We found diversity to be comparable based on Shannon's index. However, we saw a significantly lower Firmicutes/Bacteroidetes ratio in SLE individuals (median ratio, 1.97) than in healthy subjects (median ratio, 4.86; P < 0.002). A lower Firmicutes/Bacteroidetes ratio in SLE individuals was corroborated by quantitative PCR analysis. Notably, a decrease of some Firmicutes families was also detected. This dysbiosis is reflected, based on in silico functional inference, in an overrepresentation of oxidative phosphorylation and glycan utilization pathways in SLE patient microbiota. IMPORTANCE: Growing evidence suggests that the gut microbiota might impact symptoms and progression of some autoimmune diseases. However, how and why this microbial community influences SLE remains to be elucidated. This is the first report describing an SLE-associated intestinal dysbiosis, and it contributes to the understanding of the interplay between the intestinal microbiota and the host in autoimmune disorders.

INTRODUCTION

Metagenomic studies on gut microbiota burst onto the scientific scene during the last decade, due to the advent of next-generation sequencing techniques. In a very short period of time, microbiologists moved from the study of single, isolated, cultivable microorganisms, specifically, those able to grow under standard laboratory conditions, to the investigation of very complex microbial communities, mainly composed of uncultivable bacteria (1, 2). The first metagenomics reports enabled an overview of the complexity of our gut microbial communities (3, 4). Further studies focused on establishing the correlation between the human gut microbiome, the collective genomes of all microbes inhabiting the gut (5), and different physiological states, including those having an influence on health. Currently, we know that the gut microbiota might affect food and drug metabolism (6), influences human behavior (7), shifts during the course of pregnancy (8), displays age-associated changes (9–12), and possesses distinctive features depending on geographical location (12, 13), among other features. It is also becoming clear that there is a strong link between dietary patterns and the gut microbial profile (14, 15). Furthermore, some links have been established between some disorders (for example, obesity and metabolic syndrome) and an imbalance in the gut microbial ecology, also called dysbiosis (16–18). Remarkably, intestinal dysbiosis has also been associated with autoimmune diseases, such as rheumatoid arthritis, type 1 diabetes, and inflammatory bowel disease (IBD) (19–21).

Systemic lupus erythematous (SLE) is a prototypical autoimmune disease in humans that is characterized by the presence of hyperactive immune cells and aberrant antibody responses to nuclear and cytoplasmic antigens. Genetic, immunological, hormonal, and environmental factors contribute to disease susceptibility (22), and its prevalence varies greatly depending on the population under study, although a prevalence of 2 to 5 cases per 10,000 inhabitants is reportedly considered normal (23). Among the environmental factors, growing evidence suggests that molecular mimicry as a result of viral infection may contribute to the development of lupus (24). Also, some reports have highlighted intestinal infections that may ameliorate SLE symptoms (25), and a marked difference in the specificity of antibodies to bacterial DNA in healthy people and SLE patients has been indicated (26). In fact, there is early evidence of a different abundance of cultivable intestinal bacteria in SLE (27). Remarkably, it has been suggested that novel SLE biomarkers can be potentially found in the human microbiota (28). However, a study of the potential dysbiosis associated with SLE had not been tackled until now. Therefore, in this report we took advantage of next-generation sequencing techniques to explore the potential interplay of the human microbiome and SLE. We have proven, for the first time, that there is a gut microbial dysbiosis associated with SLE.

RESULTS AND DISCUSSION

Despite all the scientific knowledge generated in the last few years, and although few studies published so far support the dysbiosis theory as a key factor promoting chronic inflammation in autoimmune diseases (29–32), there is no scientific work that has taken advantage of next-generation sequencing techniques to explore the potential interplay of the human microbiome and SLE, the prototypical autoimmune disease in humans. Therefore, we designed our work with the aim of answering if there is an SLE-associated intestinal dysbiosis and, if so, which microbial population groups are related to the dysbiosis.

We defined the SLE population group by considering that there is a census of about 300 SLE patients in Asturias (from a total population of about 1,000,000 inhabitants). Thus, we were able to obtain a group of SLE patients from a well-defined geographical location to compare them with a similar group of healthy controls (HC), considering factors such as sex, age, medication (absence of antibiotic, steroid, and immunological treatments during the last 6 months), medical history (presenting a wide variety of clinical SLE manifestations), duration of the disease (2 to 24 years), and absence of flares of disease activity at the time of sampling (systemic lupus erythematosus disease activity index [SLEDAI] score of ≤8 at the time of sample collection). The group of SLE patients included individuals with a large variety of symptoms (Table 1), allowing us to establish correlations between the microbial profile and SLE, which are very likely independent of a specific pattern of symptoms. This variability in the phenotype of the disease is an intrinsic characteristic of SLE (22). We also selected patients with no active disease at the time of sampling, because the clinical manifestations of the disease in this population group are not biased by the pharmacological treatment necessary to treat SLE individuals during disease relapse. Furthermore, mean dietary intakes of energy, macronutrients, micronutrients, fiber, and phyto-compounds were recorded, both from patients with SLE and healthy subjects, and we found that there was no significant difference between the 2 groups (Table 2). Also, no significant difference was found between the 2 groups regarding lifestyle-related factors (smoking, alcohol consumption, physical activity, and use of vitamin and mineral supplements [data not shown]). This reduced the possibility that our analysis was affected by factors shown to have an influence on the gut microbial profile, such as age (9), diet (15), or phenolic compound intake (33, 34).

TABLE 1 

Demographic, clinical, and immunological features of SLE patients

Subject no.	Age (yrs)	Disease duration(yrs)	Anti-dsDNA titer (U/ml)[a]	Complement C3 (g/liter)[a]	Complement C4 (g/liter)[a]	Clinical and immunological features[b]
SLE1	43	2	0.3	0.93	0.2	MR, PH, HD, ANA
SLE2	68	3	0.7	0.96	0.18	MR, DL, PH, AR
SLE4	35	4	7.7	1.53	0.22	PH, OU, RD, ANA
SLE5	50	5	18	1.67	0.44	PH, OU, HD, ANA, anti-SSa
SLE6	35	3	48	0.81	0.13	MR, OU, AR, HD, ANA, anti-dsDNA, anti-SSa
SLE7	70	3	27	1.74	0.37	PH, OU, HD, ANA, anti-dsDNA
SLE11	54	24	99.1	1.43	0.28	MR, DL, PH, AR, HD, ANA, anti-dsDNA, anti-Sm
SLE12	58	6	13	0.84	0.22	DL, PH, AR, HD, ANA, anti-SSa, RF
SLE13	40	6	0.6	0.83	0.25	MR, OU, ANA
SLE14	40	12	4	0.92	0.18	AR, SE, RD, ANA, anti-SSb
SLE15	51	24	104	0.83	0.14	MR, DL, PH, AR, SE, HD, ANA, anti-dsDNA, anti-SSa, anti-SSb, anti-Sm, anti-RNP, anti-CLP
SLE16	54	24	45	1.76	0.3	PH, AR, ANA, anti-dsDNA, anti-SSa
SLE17	46	13	19	0.8	0.11	MR, DL, PH, OU, AR, ANA, anti-dsDNA, anti-SSa, RF
SLE18	43	12	4.1	1.04	0.16	DL, PH, OU, AR, SE, HD, ANA, anti-SSa, RF
SLE19	34	4	0	1.19	0.25	MR, PH, OU, ANA
SLE20	51	7	5.8	0.67	0.14	PH, OU, ANA, anti-SSa
SLE21	59	11	1.2	1.16	0.17	PH, ANA, anti-dsDNA, anti-SSa, anti-CLP
SLE22	64	11	4.4	1.17	0.25	MR, PH, AR, ANA, anti-dsDNA, anti-SSa
SLE24	46	14	0.4	1.08	0.4	MR, PH, HD, ANA, anti-SSa, RF
SLE26	46	20	38	0.89	0.18	MR, PH, OU, RD, HD, ANA, anti-dsDNA.

At the time of sampling.

Cumulatively registered. Abbreviations: ANA, antinuclear antibodies; anti-RNP, antiribonucleoprotein antibodies; anti-Sm, anti-Smith antigen antibodies; anti-CLP, anticardiolipin antibodies; RF, rheumatoid factor; AR, arthritis; DL, discoid lesions; HD, hematological disorder; MR, malar rash; OU, oral ulcers; PH, photosensitivity; RD, renal disorder; SE, serositis.

TABLE 2 

General characteristics and mean dietary intake of energy, macronutrients, fiber forms, and phyto-compounds in patients with SLE and healthy controls

Characteristic	SLE patients (n = 20)	Healthy controls (n = 20)
Female sex (%)	100	100
Age (yrs)	49.2 ± 10.7[b]	46.9 ± 8.6
BMI (kg/m2)	26.1 ± 5.3	25.2 ± 4.2
Energy (kcal/day)	2,173.1 ± 722.4	1,875.9 ± 332.8
Lipid (g/day)[a]	84.5 ± 41.0	85.4 ± 20.5
MUFA (g/day)[a]	35.3 ± 19.7	35.7 ± 7.6
PUFA (g/day)[a]	17.2 ± 9.7	17.5 ± 9.4
SFA (g/day)[a]	24.9 ± 14.1	25.0 ± 6.0
Protein (g/day)[a]	104.9 ± 27.6	100.6 ± 20.9
Carbohydrates (g/day)[a]	205.0 ± 75.6	203.5 ± 47.0
Dietary fiber (g/day)[a]	24.9 ± 10.4	25.3 ± 9.1
Insoluble fiber (g/day)[a]	16.0 ± 8.6	16.6 ± 7.5
Soluble fiber (g/day)[a]	2.9 ± 1.5	2.8 ± 1.1
Total isoflavones (mg/day)[a]	2.4 ± 2.4	2.5 ± 2.7
Total phenolics (mg/day)[a]	833.2 ± 527.3	916.3 ± 437.8

Model was adjusted for energy and BMI. PUFA, polyunsaturated fatty acids; MUFA, monounsaturated fatty acids; SFA, saturated fatty acids.

Values are means ± SD.

Our work is based on 16S rRNA gene-based data for fecal microbiota and the bioinformatic analysis of the results. In a previous work (35), we optimized protocols to study the human fecal microbial population by using an Ion Torrent PGM sequencing platform. This methodology was applied in the current study, and we obtained an average of 592,305 high-quality reads per fecal sample (see Table S1 in the supplemental material). Rarefaction curves obtained by plotting the Shannon, Chao1, and phylogenetic diversity indexes against the number of sequences (see Fig. S1 in the supplemental material) showed that a large part of the diversity of the samples was detected. The microbiota composition at the phylum and family levels was obtained (Fig. 1; see also Fig. S2 in the supplemental material). Remarkably, even considering the broad heterogeneity of the clinical manifestations of SLE in the individuals under study (Table 1), we obtained a particular type of microbiota for the SLE group. In this regard, the presence of anti–double-stranded DNA (dsDNA) antibodies and other clinical data were organized in a metadata file for all the microbiota profiles. A principal component analysis (PCoA) was performed with both metadata/microbiota profiles, using the variability of the 16S rRNA gene profiling at different taxonomic levels. Sample classification according to the metadata revealed no specific clustering of the samples or correlations with the different clinical features or anti-dsDNA antibodies (data not shown).

FIG 1 

Aggregate microbiota composition in fecal samples from control (HC) and lupus-affected (SLE) subjects at the phylum level (a) and family level (b). In panel b, only taxonomic groups representing >0.5% are shown.

In silico analysis of the sequences highlighted the key findings of our work. A high-quality filtering approach was used in order to process the Ion Torrent-generated sequencing data (see Table S1 in the supplemental material); a total of 293,436 unfiltered operational taxonomic units (OTUs) were identified by using uclust for de novo OTU picking. Based on each of five alpha-diversity measures (Chao1, PD whole tree, observed species, Shannon, and Simpson indexes), patients and controls were not significantly different (data not shown). Notably, one of the main results was the identification of a clear dysbiosis between the two study groups which was characterized by a higher relative abundance of Bacteroidetes in the SLE group. Overall, we detected a significant decrease in the Firmicutes/Bacteroidetes ratio in SLE individuals: the microbiota of SLE patients, compared with controls, had an almost-2.5-fold-decreased ratio (Fig. 2A) (P < 0.002). Looking at the different individuals, a gradient from lower (SLE) to higher (HC) Firmicutes/Bacteroidetes ratios was observed (Fig. 2B). These 16S rRNA gene-based analyses were corroborated by quantitative PCR (qPCR) analysis. The levels (reported as the log cells/g, with the interquartile range [IQR] in parentheses) of total fecal bacteria were 10.62 (9.46 to 10.80) in SLE patients and 10.35 (10.07 to 10.59) in controls. Firmicutes levels in the SLE and control groups reached 9.69 (8.86 to 10.38) and 9.99 (9.68 to 10.31), respectively, while those of Bacteroidetes were 10.52 (9.56 to 10.83) and 9.89 (9.59 to 10.23), respectively. No statistically significant differences in these levels of microbial groups were found between SLE patients and controls. However, when the data were expressed as the relative percentages of Firmicutes and Bacteroidetes with respect to the total bacterial level, a significantly higher (P < 0.05) percentage of Bacteroidetes was observed in the SLE group. Moreover, when the Firmicutes/Bacteroidetes ratio was calculated, a statistically significant (P < 0.01) decrease in the SLE group with respect to the control group was found (0.94 [0.90 to 0.98] versus 1.01 [0.96 to 1.06], respectively). The differences in the ratios based on qPCR were less pronounced than the differences obtained with the 16S rRNA profiling, likely because the two techniques provide different kinds of information: a relative proportion of sequences (from the 16S rRNA gene-based analysis) versus an absolute quantification of sequences (via qPCR). Thus, the fact that we obtained clear evidence for a significantly lower Firmicutes/Bacteroidetes ratio in SLE patients when we used two different culture-independent techniques supports the soundness and reliability of our analyses. The phyla Bacteroidetes and Firmicutes include the most abundant components of the human gut microbiota (36). Dysbiosis between Firmicutes and Bacteroidetes in the human gut has been described in previous studies in association with some disorders. The ratio between Firmicutes and Bacteroidetes decreases in human type 2 diabetes compared with controls (37). Also, most studies of the microbiota in people with Crohn’s disease report a decrease in the abundance of Firmicutes and an increase in Bacteroidetes in association with the disease (38). An opposite situation is observed in obesity, in which the dysbiosis is characterized by an increase in the Firmicutes/Bacteroidetes ratio (17). Therefore, this specific microbial balance between the more abundant phyla in the human gut seems to be dependent on the physiology of the disorder. In relation to this, it has been reported that this equilibrium is susceptible to modification by shifts in the dietary pattern. Wu et al. (15) reported that Firmicutes levels were positively associated with a low-fat/high-fiber diet. Also, dietary interventions, including whole grain in the diet, increase the Firmicutes/Bacteroidetes ratio (39).

FIG 2 

(A) Box plot of Firmicutes/Bacteroidetes ratios (median ± IQR) in SLE patients versus healthy controls. (B) Percentages of 16S rRNA reads of Bacteroidetes (red bars) and Firmicutes (blue bars) in the DNA extracted from fecal samples of SLE patients (SLE codes) and healthy controls (HC codes). Ratios are significantly different (P < 0.002).

We prepared scatter plots (see Fig. S3A in the supplemental material), and they clearly highlighted a positive association between Firmicutes in healthy controls (P < 0.01) and Bacteroidetes in SLE patients (P < 0.001). This association was confirmed at lower taxonomic levels, and normalized abundances of the classes Bacteroidia (phylum Bacteroidetes) and Clostridia (phylum Firmicutes) and the orders Bacteroidales and Clostridiales differed significantly between the SLE patients and the healthy controls (see Fig. S3B). At the family level, Lachnospiraceae (P < 0.05) and Ruminococcaceae (P < 0.05) were positively associated with healthy controls (see Fig. S4 and S5 and Table S2 in the supplemental material).

Statistical differences between the two groups (HC and SLE) were calculated by a PERMANOVA test, with the distance data obtained after ordination analysis (PCoA) using PAST v 3.1 (40). In all cases, data from the relative taxa abundances were used, and distances were computed according to the Bray-Curtis similarity index. The groups of SLE and HC differed statistically whenever phylum- or family-level data were used (P < 0.01 or P < 0.02, respectively).

Unsupervised PCoA of the 16S rRNA sequence data identified the phyla Bacteroidetes and Firmicutes as the main gradients for SLE and healthy control groups, respectively (Fig. 3A). PCoA at the family level showed that the families Lachnospiraceae and Ruminococcaceae were located near the healthy controls (Fig. 3B). Interestingly, Lachnospiraceae was the most abundant family in the feces of both study groups (Fig. 1; see also Fig. S2 in the supplemental material). Depletion of Lachnospiraceae and Ruminococcaceae has been associated with Clostridium difficile infections and nosocomial diarrhea (41). Also, several studies showed a decrease of Lachnospiraceae in IBD patients, and this family has been suggested as a biomarker of disease activity (42–44). However, the functional consequences of the depletion of these bacteria in the previously mentioned intestinal diseases remain to be investigated. Although the amplicons of the 16S rRNA sequences could be relatively short to perform a totally reliable population structure analysis at the genus/species level, it is noteworthy that the relative abundance of sequences assigned to Bacteroides spp. were significantly higher in SLE samples (P < 0.02). In the control group, we also found a significantly higher (P < 0.05) relative abundance of sequences tentatively assigned to Desulfovibrio, the most common genus of sulfate-reducing bacteria in the human gut (45). It is worth mentioning that some authors have described how dysbiosis could affect mucosal barrier function and impair immunoregulatory mechanisms, leading to pathological effects in systemic immunity. Using animal models, a direct involvement of components of the microbiota in chronic intestinal inflammation (46) and the protective role of specific commensals in avoiding bacterial translocation (47) have been demonstrated. Finally, we should bear in mind that nonsteroidal anti-inflammatory drugs and the antimalarial treatment of SLE patients could have an influence on the observed dysbiosis.

FIG 3 

PCoA results for the 16S rRNA profiles at the phylum (A) and family (B) level. The presence or absence of SLE was further included as metadata. HC, closed circles; SLE, open triangles.

In our study, we also determined the tentative metagenomes from phylogenetically associated reference genomes. Our aim was to highlight metabolic pathways and shifts associated with the SLE population compared with controls. We inferred the functionality of the different putative metagenomes by using PICRUSt software, which allows the prediction of metabolic pathways from the 16S rRNA reads (48). A functional analysis using the data obtained from the KEGG pathways at level 3 allowed us to detect certain processes potentially associated with either healthy controls or SLE. KEGG levels are the different hierarchical subdivisions in which the functions of a biological system (cell, organism, or ecosystem) are arranged according to information organized in the KEGG database (http://www.genome.jp/kegg/kegg1a.html). Pathways displaying a difference in mean proportions between healthy controls and SLE groups of at least 0.1% are represented in Fig. S6 in the supplemental material. Some glycan degradation pathways are slightly overrepresented in the microbiota of SLE patients, probably due to the higher abundance of Bacteroidetes in these samples. Bacteroidetes, and specifically the main genus of this phylum, Bacteroides, have been shown to display broad glycan-degrading abilities (49). The same occurs with lipopolysaccharide biosynthesis proteins, which is in direct relation to the higher abundance of Bacteroidetes in the SLE samples. Remarkably, oxidative phosphorylation processes seem to be associated with SLE patients. This finding indicates that some bacteria able to perform oxidative phosphorylation may be better adapted to the intestinal ecosystem of individuals with SLE, and this could be related to the imbalance in the oxidative stress environment at the intestinal level that has been linked to some autoimmune diseases (48).

In summary, understanding and potentially manipulating immune responses through the action of intestinal microbiota comprise one of the most active fields in probiotic and prebiotic research. Most likely, the dysbiosis defined in this work is the consequence of the altered immune function of SLE patients. At present, the treatment of SLE patients is exclusively performed with drugs, and our findings could indicate that challenging the immune system with the bacteria depleted in SLE, or substances promoting their growth, could influence SLE physiology. To the best of our knowledge, experimental papers about the relationship of the human gut microbiome and SLE have not been published, and this is the first report describing an SLE-associated intestinal dysbiosis. Thus, our results establish the basis to delve deeper into the understanding of the relationship between the human gut microbiota and autoimmune diseases.

MATERIALS AND METHODS

Ethics statement.

Ethics approval for this study (reference code AGL2010-14952; grant title “Towards a better understanding of gut microbiota functionality in some immune disorders”) was obtained from the Bioethics Committee of CSIC (Consejo Superior de Investigaciones Científicas) and from the Regional Ethics Committee for Clinical Research (Servicio de Salud del Principado de Asturias) in compliance with the Declaration of Helsinki. All determinations were performed with fully informed written consent from all participants involved in the study.

Study subjects.

The study sample comprised 20 patients with SLE (SLE codes) and 20 healthy controls (HC codes). SLE patients were selected from the updated Asturian Register of Lupus (Asociación Lúpicos de Asturias, Oviedo, Spain). All of them fulfilled at least four of the American College of Rheumatology criteria for SLE (50). The individuals recruited for the SLE group were in remission and had not had any immunotherapy or corticoid treatment during the previous months, since those are the treatments that could have the strongest influence on the physiology of the patients. All patients were women of Caucasian origin, 49.2 ± 10.7 years old (mean ± standard deviation [SD]), and had no active disease at the time of sampling (SLEDAI score of ≤8). Information on clinical manifestations was obtained by reviewing clinical records (Table 1). Patients were also asked precise questions regarding treatment received during the previous 6 months. Only those individuals who had not used antibiotics, glucocorticoids, immunosuppressive drugs, monoclonal antibodies, or other immunotherapies were recruited for the study. Eighteen patients were receiving antimalarial treatment, and all of them were regular consumers of nonsteroidal anti-inflammatory drugs. Twenty age-matched healthy women (46.9 ± 8.6 years old) from the same population were recruited as controls.

(i) Nutritional assessment.

Variables of macro- and micronutrient intakes were collected by means of a semiquantitative food frequency questionnaire that included 160 items. During a personal interview, carried out within 7 days after the collection of the fecal sample, subjects were asked item by item whether they usually ate the food and, if so, how much they ate. For this purpose, three different serving sizes of each cooked food were presented in pictures to the participants so that they could choose from up to 7 serving sizes (from “less than the small one” to “more than the large one”). For some of the foods consumed, amounts were recorded in household units, by volume, or by measuring with a ruler. Special attention was paid to cooking practices and number and amount of ingredients used in each recipe, as well as questions concerning menu preparation (e.g., type of oil or milk used).

Food intake was analyzed for energy and macro- and micronutrient contents by using the nutrient food composition tables developed by CESNID (51). Total, soluble, and insoluble fiber intake was completed based on Marlett food composition tables (52), and polyphenol content was calculated from the U.S. Department of Agriculture (USDA) Nutrient Database (53).

(ii) Anthropometric measures.

The body mass index (BMI) was calculated from the following formula: weight/(height)2 (in kg/m2). Height was registered by using a stadiometer with an accuracy of ±1 mm (Año-Sayol, Barcelona, Spain). Subjects were barefoot, in an upright position, and with the head positioned in the Frankfort horizontal plane. Weight was measured on a scale with an accuracy of ±100 g (Seca, Hamburg, Germany).

(iii) Lifestyle-related factors.

During the interview, other factors associated with the lifestyle of the subject were registered. Smoking habit, physical activity, alcohol consumption, and supplements use were included in the questionnaire. Smoking status was categorized as nonsmoker (including exsmokers and occasional smokers) or current smoker. Those subjects who reported that they never exercised were categorized as physically inactive. “Regular alcohol consumer” refers to those subjects who declared a regular consumption of alcoholic drinks. Also, the use of vitamin and mineral supplements during the last month was self-reported.

Fecal sample collection and DNA extraction.

Fresh fecal material (between 10 and 50 g per person) was collected in a sterile container and immediately manipulated and homogenized within a maximum of 3 h from defecation. During the waiting period, from defecation to homogenization, samples were kept at 4°C. Thirty milliliters of RNAlater solution (Applied Biosystems, Foster City, CA) was added to 10 g of sample, and the mixture was homogenized in a sterile bag, using a stomacher apparatus (IUL Instruments, Barcelona, Spain) with three cycles at high speed, 1 min per cycle. Homogenized samples were then stored at −80°C until use.

For DNA extraction, samples were thawed and the QIAamp DNA stool minikit (Qiagen Ltd., Strasse, Germany) was used as previously described (35).

16S rRNA gene amplification.

Partial 16S rRNA gene sequences were amplified from extracted DNA by using the primer pair Probio_Uni (5′-CCTACGGGRSGCAGCAG-3')/Probio_Rev (5′-ATTACCGCGGCTGCT-3′) (35), which targets the V3 region of the 16S rRNA gene sequence. The PCR conditions used were 5 min at 95°C, 35 cycles of 30 s at 94°C, 30 s at 55°C, and 90 s at 72°C, followed by 10 min at 72°C. Amplification was carried out using a Verity thermocycler (Applied Biosystems). The integrity of the PCR amplicons was analyzed by electrophoresis on an Experion workstation (Bio-Rad, Hertfordshire, United Kingdom).

Ion Torrent PGM sequencing of 16S rRNA gene-based amplicons.

The PCR products derived from amplification of specific 16S rRNA gene hypervariable regions were purified by electrophoretic separation on a 1.5% agarose gel and the use of a Wizard SV Gen PCR cleanup system (Promega, Madison, WI), followed by a further purification step involving Agencourt AMPure XP DNA purification beads (Beckman Coulter Genomics GmbH, Bernried, Germany) in order to remove primer dimers. The DNA concentration of the amplified sequence library was estimated through use of the Experion system (Bio-Rad). From the concentration and the average size of each amplicon library, the amount of DNA fragments per microliter was calculated and libraries for each run were diluted to 3E9 DNA molecules prior to clonal amplification. Emulsion PCR was carried out using the Ion OneTouch 200 template kit v2 DL (Life Technologies, Guilford, CA) according to the manufacturer’s instructions. Sequencing of the amplicon libraries was carried out on 316 chips by using the Ion Torrent PGM system and employing the Ion Sequencing 200 kit (Life Technologies) according to the supplier’s instructions at the DNA sequencing facility, GenProbio s.r.l. After sequencing, the individual sequence reads were filtered with the PGM software to remove low-quality and polyclonal sequences. Sequences matching the PGM 3′ adaptor were also automatically trimmed. All PGM quality-approved, trimmed, and filtered data were exported as SFF files.

Sequence-based microbiota analysis.

The SFF files were processed using QIIME (54). Quality control retained sequences had lengths between 150 and 200 bp, a mean sequence quality score of >25, and with truncation of a sequence at the first base if a low-quality rolling 10-bp window was found. Presence of homopolymers of >7 bp and sequences with mismatched primers were omitted. In order to calculate downstream diversity measures (alpha and beta diversity indices, Unifrac analysis), 16S rRNA OTUs were defined at ≥97% sequence homology by using uclust (55). Chimeric sequences were removed using ChimeraSlayer (56). Furthermore, OTUs that included less than 10 sequences were filtered using QIIME (54). All reads were classified to the lowest possible taxonomic rank by using QIIME and a reference data set from the Ribosomal Database Project (57). The sequence data features of all the samples are included in Table S1 in the supplemental material.

Different alpha diversity indexes were calculated using QIIME and information from the OTU tables using the alpha_diversity.py script. The different diversity metrics were set by passing the option –s to the script. The following indexes were calculated for every sample and compared between groups by using a two-sided Student’s t test: Chao1, PD whole tree, observed species, Shannon, and Simpson.

Analysis by qPCR.

Quantification of total fecal bacteria, Firmicutes, and Bacteroidetes by qPCR was performed by using previously described primers and conditions (58–60). Analyses were done in duplicate in a 7500 Fast real-time PCR system (Applied Biosystems) using Sybr green PCR master mix (Applied Biosystems). Standard curves were made with pure cultures, grown under anaerobic conditions at 37°C, of Escherichia coli LMG 2092 in Gifu anaerobic medium (GAM; Nissui Pharmaceutical Co., Tokyo, Japan), Faecalibacterium prausnitzii DMSZ 17677 in RCM formula (Oxoid Ltd., Basingstoke, Hampshire, United Kingdom), and Bacteroides thetaiotaomicron DSMZ 2079 in GAM.

Functional inference.

The functionality of the different metagenomes, grouped by disease status (healthy control versus SLE), was predicted using the software PICRUSt 0.9.1 (http://picrust.github.io), which has been explained in detail elsewhere (48). In short, this software allows the prediction of functional pathways from the 16S rRNA reads. First, a collection of closed-reference OTUs was obtained from the filtered reads by using QIIME v 1.7.0 (54) and by querying the data against the IMG/GG reference collection (GreenGenes database, May 2013 version; http://greengenes.lbl.gov). Reverse-strand matching was enabled during the query, and OTUs were picked at a 97% identity. A BIOM-formatted table (biological observation matrix [61]) was obtained with the pick_closed_reference_otus.py script. This table, containing the relative abundances of the different reference OTUs in all the metagenomes, was normalized based on the predicted 16S rRNA copy number by using the script normalize_by_copy_number.py. Final functional predictions, inferred from the metagenomes, were created with the script predict_metagenomes.py. When necessary, tab-delimited tables were obtained with the script convert_biom.py.

Analysis of predicted metagenomes.

PICRUSt and QIIME provide a number of scripts that can be useful for analyzing both 16S rRNA gene relative abundances and the predicted metabolic data. Predicted metagenomic contents were collapsed at KEGG pathway level 3 (http://www.genome.jp/kegg/pathway.html) with the categorize_by_function.py script, and the data were analyzed statistically by using STAMP v 2.0.0 (62). STAMP allows data filtering and the application of different statistical tests and corrections, including PCoA. It also generates different graphics, including box plots, error plots, and scatter plots. Data of the KEGG pathway distributions were plotted by using the script summarize_taxa_through_plots.py. Associations of different taxonomic categories to SLE were statistically analyzed with the script otu_category_significance.py.

Statistical analyses.

Statistical analysis was performed using IBM-SPSS version 19.0 (IBM SPSS, Inc., Chicago, IL). For descriptive purposes, in Table 2 the mean values are presented as means ± SD on untransformed variables. Differences between SLE patients and controls were compared by using a multivariate linear model, including energy intake and BMI as covariates.

Individuals were ordered according to their sequence data composition by principal component analysis using the taxonomic data at the phylum and family levels. Patterns were extracted using all the variations from the taxonomic data via an indirect method as a model and SLE as metadata. To analyze the associations between inferred metabolic pathways and SLE, metabolic pathways with very low abundance levels (<0.001 in 50% of the samples) were excluded from all analyses. Association of KEGG pathways to SLE were identified by running two-sided Welch tests on every pair of means. This test is a variation of Student’s t test and is used when equal variance cannot be assumed in both groups. Confidence intervals (95%) were obtained by inverting the Welch’s tests. The false discovery rate (FDR) correction (48, 63) was finally applied in all cases, and significant differences between healthy controls and SLE patients were only considered when below a P value of 0.05 or a q value below 0.2. P and q values at the phylum and family levels are included in Table S2 in the supplemental material. In the particular case of the family Desulfovibrionaceae (with P < 0.05), further statistical analysis was carried out for sequences tentatively assigned to the genus Desulfovibirio by using a one-sided t test.

In relation to qPCR results, not all the bacterial groups showed normal distribution; therefore, differences in bacterial levels between groups of individuals were analyzed using a nonparametric test (Mann-Whitney U test).

Nucleotide sequence accession number.

The raw sequences reported in this article have been deposited in the NCBI Short Read Archive (SRA; study accession number SRP028162).

Rarefaction curves generated for 16S rRNA gene sequences obtained from fecal samples of control subjects (HC samples; gray lines) and lupus-affected subjects (SLE samples; red lines). (a) Rarefaction curves derived using the Shannon index. (b) Rarefaction curves derived using the Chao1 index. (C) Number of phylotypes identified versus the number of sequences per sample. Download

Figure S1, DOC file, 0.2 MB

Interindividual variation in fecal samples of SLE subjects (a) and control subjects (b) in the proportion of the major microbial phyla/families. The aggregate microbiota composition of control (HC codes) and lupus-affected (SLE codes) subjects at the levels of phylum and family are indicated at the bottom of each inset. Only taxonomic groups representing >1% are shown in the keys. Download

Figure S2, DOC file, 0.1 MB

(A) Scatter plot at the phylum level in which the relative abundances (means) of Firmicutes and Bacteroidetes in healthy controls (x axis) and SLE patients (y axis) are represented. (B) The normalized abundances of the classes Bacteroidia and Clostridia and the orders Bacteroidales and Clostridiales differed significantly between the lupus (SLE) and the healthy control (HC) 16S rRNA sequence data, as represented in the error plots. Tables with relative abundances and sequence metadata were elaborated using QIIME and analyzed using STAMP v 2.0.0. Two-sided Welch’s tests were run on every pair of means. The axes represent mean proportions of the taxa in HC (x axis) or SLE (y axis) groups. The vertical and horizontal lines in each taxa represent data spread (2nd and 98th percentiles) in the HC group (horizontal lines) and in the SLE group (vertical lines). The statistical fit of the points to the gray dashed y = x line is represented by the coefficient of determination in the upper left part of the plot (R2). The y = x line represents the boundary delimiting taxa that are enriched in the HC group (blue dots) or in the SLE group (yellow dots). Download

Figure S3, DOC file, 0.3 MB

Scatter plot at the family level in which the relative abundances (means) of the bacterial families in healthy controls (x axis) and lupus patients (y axis) are represented. Families showing the greatest relative abundances and divergence from the diagonal are labeled with arrows. Download

Figure S4, DOC file, 0.2 MB

Summaries of the communities present in the different samples analyzed in this study. Plots were obtained at four taxonomic levels: phylum (A), class (B), order (C), and family (D), and the most abundant taxons are labeled with arrows. Data were obtained using QIME v 1.7.0 with the scripts pick_de_novo_otus.py and summarize_taxa_through_plots.py. Download

Figure S5, DOC file, 0.5 MB

KEGG pathways (hierarchical level 3) for which normalized abundances differed between SLE patients and healthy controls. Functional metagenomes were predicted using PICRUSt scripts, and only pathways reaching statistically significance are shown. On the y axis, the KEGG pathways are indicated. Mean proportions for each metadata group (healthy controls versus SLE) are shown together with the confidence intervals of the data. Variations due to every single metagenome on mean proportions are shown by horizontal error bars. Download

Figure S6, DOCX file, 0.2 MB

Sequence data features

Table S1, PDF file, 0.04 MB.

Statistically significant differences between healthy controls and SLE patients at the phylum and family levels

Table S2, PDF file, 0.4 MB.