Community differentiation of the cutaneous microbiota in psoriasis

Study population and subject specimen matching

Between June 2008 and September 2011, we obtained consent (using the model consent forms for the HMP demonstration projects) and enrolled a total of 199 subjects (75 patients with psoriasis and 124 healthy controls) with ethical approval from New York University School of Medicine Institutional Review Board (IRB #08-709). The subjects were recruited from the same geographic region (NYC) and same clinic at NYU. Among the patients with psoriasis, 57 (76%) had not been exposed to antibiotics or received treatment relevant to psoriasis for at least 1 month before skin samples were obtained. Among the healthy controls, 112 (90.3%) had not been exposed to antibiotics or received treatments relevant to psoriasis for at least 1 month before skin samples were obtained. Psoriasis subjects receiving antibiotics less than 1 month before enrollment were excluded from further analysis, only six (11.8%) of the remaining subjects had taken any antibiotics in the preceding year. When we reviewed the control subjects who actually were included, none had received antibiotics in the 12 months prior to sampling. A total of 54 (72%) patients with psoriasis were studied by swabbing of the affected (lesion) and unaffected (unaffected) sites (see section on Specimen collection for details).

For these subjects, we sought control subjects of the same gender and ethnicity, and of similar age (± 5 years), from whom a cutaneous specimen was obtained in a region proximate to the site of the psoriasis lesion. In total, we obtained matching specimens from 37 (29.8%) of the control subjects. One or more sites from each of these controls were matched to the lesions in the 54 subjects with psoriasis. A control subject could be matched to more than one patient, since we also matched for cutaneous site. However, each control cutaneous site was uniquely mapped to only one triplet, thus there was no duplication of specimens in the analysis. In each of 48 matched pairs, the 2 sites match, but in 6 sites we matched a back specimen with an abdomen, which are relatively similar in composition in healthy skin. The final analyses were performed on a set of 51 triplets, which had adequate depth of sequencing (>1,000 sequences per sample).

The resulting set of 51 triplets contained samples from sites that are characteristic of where typical psoriatic lesions occur in the general population. All of the sites were of the dry or sebaceous cutaneous microenvironments. We grouped the exact location of the specimen by proximity to other samples into four categories: body, head, upper extremities and lower extremities. Upper and lower extremities contained only samples of the dry cutaneous microenvironment, while all head samples and 8 out of 12 body trunk samples were characterized as sebaceous. A table describing the matching of psoriasis lesions to control sites and skin environment is provided in Additional file 1: Table S1.

Although psoriasis affects each gender equally, our final set of subjects consisted of 75% males. The bias towards men being sampled possibly can be explained by the fact that the medicines used in the study are often not used in women of childbearing age, thus limiting the enrollment of women.

Psoriasis diagnosis and characteristics of populations

Patients were diagnosed with chronic plaque psoriasis in a dermatology clinic, and psoriasis was clinically classified based on characteristic morphologic features of the individual skin lesions and their distribution on the body. For each patient, disease duration, percentage cutaneous involvement, Psoriasis Area and Severity Index (PASI) and physician global assessment (PGA) scores were recorded. Means of severity scores for the subjects were PASI: 8.7 (± 10.1 SD), PGA: 6.6 (± 6.9 SD), and body surface area (BSA): 9.4 (± 13.9 SD). The characteristics of the control and affected study populations are given in Additional file 1: Table S3.

DNA sequencing and upstream processing

Samples were prepared for amplification and sequencing at the J. Craig Venter Institute (JCVI) Joint Technology Center (JTC) using a protocol for 16S rRNA gene amplification and sequencing developed as part of the NIH Human Microbiome Project [18, 21]. Negative control experiments were performed, when we developed the extraction protocol with Qiagen kit [22]. In short, we used a reagent control that included all DNA extraction and polymerase chain reaction (PCR) reagents, including the sterile swab and the buffers, without the skin sample. This specimen was examined in parallel using the identical procedures as with the skin samples. After electrophoresis and ethidium bromide staining, preparations from these controls did not generate any visible bands. Negative control reactions were performed for every pool of amplicons to ensure no visible detection of amplicons on ethidium bromide stained agarose gels. The V1 to V3 region of the 16S rRNA gene was amplified using forward primer 5′-AGAGTTTGATCCTGGCTCAG-3′ attached to the Roche B adapter for 454-library construction and reverse primer 5′-CCGTCAATTCMTTTRAGT-3′ attached to the Roche A adapter and a 10-nt barcode (5′-A-adapter-N (10) + 16S primer-3′). The barcoded primer design was completed using a set of algorithms developed at the JCVI for these purposes [23, 24]. PCR reactions were completed as follows (per reaction): 2 μl of gDNA (approximately 2 to 10 ng/μl), 1× final concentration of Accuprime PCR Buffer II (Invitrogen, Carlsbad, CA, USA), 200 nmol forward and reverse primers, 0.75 U of Accuprime TaqDNA polymerase high fidelity (Invitrogen), and nuclease-free water to bring the final volume to 20 μl. PCR cycling conditions consisted of an initial denaturation of 2 min at 95°C, 30 cycles of 20 s at 95°C, and 30 s at 56°C followed by 5 min at 72°C. A high number of amplification cycles is standard for skin studies because of typically low bacterial load in these specimens [25]. A negative control (water blank) reaction was examined after 35 cycles, and determined to be negative for the amplicon. Samples were then quantified, cleaned, and sequenced on the Roche 454-FLX (454 Life Sciences, Branford, CT, USA) as described previously [21], and a read processing pipeline consisting of a set of modular scripts designed at the JCVI were employed for upstream processing, consisting of deconvolution, trimming, and quality filtering, as described previously [26]. We performed a parallel analysis of the V3 to V5 16S rRNA gene region, but because of amplification and sequencing depth issues there were only 21 available triplets at this locus. Therefore, we focused exclusively on the V1 to V3 dataset.

Downstream sequence processing and statistical analyses

After upstream processing and quality checking the passing sequences were analyzed using QIIME scripts [27]. We first clustered the sequences into 97% identity operational taxonomical units (OTUs) using the UCLUST program [28]. A representative sequence from each OTU cluster was used to assign taxonomy to the cluster using the RDP Classifier [29] executed at 80% bootstrap confidence cut-off. These representative sequences were further aligned using PyNAST [30] with the Greengenes core-set alignment template. We used the alignment to reconstruct an approximate phylogenetic tree using FASTTREE [31]. The obtained phylogenetic tree and abundance tables were used to calculate unweighted and weighted UniFrac β diversity indices [32]. The OTU absolute abundance table and UniFrac β diversity matrices were extracted from the pipeline for further analysis in the R statistical programming environment [33]. After processing the median sequencing depth per sample was 8,621 (IQR 5,013 to 11,412). The sequencing effort was statistically similar across clinical skin types, body sites and cutaneous microenvironment (Additional file 1: Figure S1).

Chimeras were checked with ChimeraSlayer [34]. In all, 2,700 of the total 34,123 OTUs (7.9%) identified in the study were marked as potentially chimeric. On average, the total relative abundance (fraction of total sequences) per sample of putatively chimeric sequences was 3% (± 2% SD). The abundance of suspected chimera was similar across clinical skin types, body sites and cutaneous microenvironment (P >0.05 Kruskal-Wallis analysis of variance (ANOVA)).

The rarefactions for richness and Shannon diversity indices were calculated using scripts based on the community ecology package vegan. Comparisons of intergroup and intragroup β diversity were performed using one-way ANOVA with the Tukey honestly significant difference (HSD) multiple comparison correction procedure.

We used the ade4 package [35] in R to perform Principal Coordinates Analysis (PCoA) on weighted Unifrac distances. To avoid negative eigenvalues in the analysis, we used the Cailliez method [36] to convert the weighted Unifrac distance matrix into a closest corresponding matrix with Euclidean properties, which was further used for PCoA.

Univariate testing was performed on OTU relative abundances, calculated by dividing the absolute abundances by the total sequence count per sample analyzed. Differential relative abundance of specific taxa and OTUs was calculated on highly abundant taxa (mean relative abundance >1%) using the Kruskal-Wallis test with FDR correction for multiple testing [37]. This approach is analogous to standard ANOVA in that the test is significant if any pair of relative abundances (control vs unaffected, control vs lesion, unaffected vs lesion) is different. Post hoc pairwise testing with additional multiple testing control can be utilized to determine which pair is different.

Multiple testing correction and compositional data issues

The fact that the relative abundances present a compositional constraint violates the independence assumption. The relevant nature of the independence violation is that the individual significance values for the univariate tests are now potentially positively correlated. An example of such correlations is a situation where a statistically significant increase in one taxon abundance between two conditions is accompanied by a balancing decrease in one or more other taxa to have the abundances sum to a constant (1). However, the effect of this independence violation on the validity of the univariate findings is only mild for the following reasons: (1) the compositional constraint does not remove any true positive association, it only inflates the false negative rates, (2) false negatives are then controlled by the FDR multiple testing correction procedure, which is designed to take into account positive correlations in P values; (3) the effect of compositional constraint is minimized by the fact that we only focus on highly abundant taxa. Therefore, we believe that this study admits false positives at a rate similar to other genomic analyses, and this allows for discovery of useful associations with the phenotypes, which may be of potential diagnostic value, but may need further validation.

We utilized univariate χ² tests to compare the prevalence of specific taxa among clinical skin types. Spearman correlation tests were used to find associations between severity scores and taxa abundance. The P values were adjusted for false discovery using the Benjamini-Hochberg procedure [37]. Receiver operating characteristic curves (ROC) were computed in R using the package ROCR [38] and significance of the classification signal as measured by the area under the ROC curve (AUC) was established by Mann–Whitney test.

To establish the presence and identity of cutaneotypes in our data, we utilized methodology identical to that previously used for gut-microbiota enterotype classification [39]. In short, we applied the partitioning around medoids (PAM) method [40] to the square root of the Jensen-Shannon divergence distances to compute optimal clustering with given numbers of clusters (2 through 20). The Calinski-Harabasz index [41] was used to establish the number of cutaneotypes to optimally cluster the data. Additional evidence for clustering was obtained using the gap statistic [42], and is described in supplementary materials.

Non-Euclidean multivariate analysis of variance (MANOVA) was used to analyze the association of microbiota with clinicodemographic variables [43]. This analysis utilizes a matrix of squares of arbitrarily computed pairwise distances in lieu of the covariance matrix to be decomposed into within and between group sums of squares. This decomposition is used to compute a pseudo-F-statistic, the significance of which may be established by permutation. Post hoc pairwise testing of significant multilevel factors was likewise performed by permutation. This analysis was performed using the adonis function from the R package vegan [44].

Psoriatic lesions trend to decreased taxonomic diversity

To assess the changes in α diversity related to psoriasis, we examined the diversity of cutaneous microbiota at the three clinical skin types, in terms of taxonomical richness and evenness at phylum, class, order, family and genus taxonomical level, as well as at the level of OTUs defined as sequences with 97% identity to each other. At all taxonomical levels, the richness index demonstrated a trend of decreasing diversity from controls to unaffected to psoriatic lesion skin (Figure 1A). We tested the significance of this trend at depth of 2,000 sequences per sample; however, this trend was not statistically significant (P >0.05). The shape of saturation curves suggests that while our data accurately reflect genus-level and higher-level composition of cutaneous microbiota, our data do not encompass all subgenera-level skin taxa. In terms of evenness as measured by the Shannon index (also known as non-normalized Shannon equitability index) (Figure 1B), saturation was reached faster at each taxonomic level. Lesion specimens had less evenness (as measured by the Shannon index) than the unaffected and control specimens; we observed statistically significant differences in evenness (P <0.05 at the phylum, class, order, family and genus levels, however the differences were not significant at the 97% identity OTU level, P >0.05). The specimens from all three clinical skin types share a large fraction of taxa at all taxonomical levels (Figure 1C). Essentially no taxa uniquely characterize any of these skin types. Decreased evenness in the lesion microbiota along with maintenance of approximately the same number of total taxa predict that one or more taxa should exhibit higher relative abundance in the lesion samples, compared to the other groups.

Psoriasis status is associated with relative abundance and presence of specific taxa

Next, we examined the relative abundance of the taxa that constitute the multivariate distribution of the cutaneous microbiota. Skyline plots demonstrated generally similar compositions of cutaneous microbiota in terms of phylum-level relative abundance (Figure 2A) for the three groups of specimens. We observed that three phyla: Proteobacteria, Firmicutes, and Actinobacteria dominate the skin microbial communities in all three skin types, consistent with observations of previous studies of psoriasis and healthy skin composition [2]. At the genus level (Figure 2B), lesional specimens show similar composition to other groups with no strongly apparent differences. We further examined the dominant taxa (those with mean relative abundance of ≥1%) at each taxonomic level to assess their association with psoriasis status (Table 1 and Additional file 1: Figure S2). Many taxa were found to be associated with psoriasis status, despite high variance in each clinical group. Surprisingly, all of the highly abundant taxa that were significantly different showed a decrease in relative abundance in lesion specimens.

We examined the abundances of the major skin genera with respect to psoriasis status. Each of the major taxa that typically are found on skin (Propionibacterium, Corynebacterium, Streptococcus, and Staphylococcus), were not significantly different between lesion, unaffected, and control. However, the combined relative abundance of these four genera was significantly (P <0.01) different across the specimen groups. Upon further examination Propionibacterium does not play an important role for distinguishing skin types. Combined relative abundance of just three genera (Corynebacterium, Streptococcus, and Staphylococcus) attained statistical significance (P <0.05). The mean combined relative abundance of these genera increases from control (mean (± SEM): 22.03% (± 2.1%)) to unaffected (22.9% (± 2.5%)) to lesion (33.8% (± 3.3%)) specimens. Pairwise post hoc testing revealed that the combined abundance of the three genera in control and unaffected microbiota was different from lesion (P <0.05). Likewise, the univariate classification signal, as measured by AUC, for each of the four skin-associated taxa was not significant, while the combined signal of all four and just three (without Propionibacterium) as well as the univariate signals of other named differentially abundant taxa were stronger, approaching diagnostically relevant values (AUC 0.65 to 0.81) (Additional file 1: Figure S3). To place these results in context, the widely used prostate-specific antigen (PSA) has an AUC <0.75, and for some versions <0.65 [45]. We have further analyzed this and an additional dataset for multivariate signatures of psoriasis and the results suggest that even stronger reproducible signals are possible [46].

Psoriasis lesions are characterized by greater intragroup variability

We evaluated the skin microbiota for the degree of sample-to-sample variability using β diversity, a distance-based ecological measure that allows for comparison of specimens grouped according to skin type. The β diversity analyses were performed, using subject site as the unit of independent analysis, however, for lack of appropriate comparison methodology these analyses do not account for the fact that some control subjects contributed multiple specimens (from different sites) to the matched triplets. We expect the intrasite differences to be greater than the differences between subjects sampled at a single site, which at least partially alleviates this apparent loss of independence. We chose to examine the β diversity in our data using weighted and unweighted Unifrac dissimilarity measures [32] based on their performance [47] in showing the differences among clinical groups. In terms of both measures, lesional specimens showed the highest intragroup diversity, followed by the unaffected, and control specimens (Figure 3). These findings indicate that control subjects are more similar to one another in terms of their cutaneous microbiome than the specimens from unaffected sites of psoriasis patients, and that each lesion harbors a more distinct community of microbes compared to other lesions. The intergroup unweighted Unifrac distances are smaller between lesion and control than between lesion and unaffected sites. Similarly, the unaffected skin is closer to the control than to lesion. These findings suggest that the microbiota of unaffected skin from psoriasis patients maintains many compositional characteristics of the skin of the control subjects, while simultaneously showing divergence from lesion skin on the contralateral part of the body. This difference is indicated by the substantial β diversity of the lesional microbiota. Using the weighted Unifrac analyses, trends are similar, with increasing β diversity for lesion > unaffected > control specimens (Figure 3B). Distances in relation to the lesion specimens are greater in relation to the control and unaffected specimens than are the comparisons of control and unaffected. From these findings, we conclude that psoriasis is associated with a systemic change in the cutaneous microbiota that is evident in both the lesion samples and, to a lesser degree, in the clinically unaffected skin.

The psoriatic microbiota is associated with a cutaneotype enriched for Firmicutes and Actinobacteria

For clinical and investigational purposes, it would be useful if the cutaneous microbiota could be stratified into distinct subtypes. We performed a preliminary analysis of the cutaneous microbial community data for evidence that there might be such stratification into distinct clusters, using the methodology previously employed to report clustering of the gut microbiota [39]. Based on the Calinski-Harabasz index [41], we show that the data can be optimally separated into two clusters, which we call cutaneotypes (Figure 4A). A representation of these clusters on the first PCoA plane is shown in Additional file 1: Figure S4. These cutaneotypes differ in terms of the relative abundance of major phyla (Figure 4B); specifically, cutaneotype 1 is dominated by Proteobacteria (Figure 4D), while communities of cutaneotype 2 have higher relative abundance of Actinobacteria (Figure 4E), and Firmicutes (Figure 4F). We notice that the distribution of Proteobacteria is bimodal, with peaks corresponding to the modes of the respective contributing distributions, each linked to a cutaneotype. This suggests that changes in Proteobacteria may be driving the clustering, while the changes in relative abundance of the other two phyla may be due to compositional effects. Although the actual value of Proteobacteria relative abundance may vary continuously over the domain (0% to 100%), the bimodality of the Proteobacteria distribution serves as an important additional piece of evidence for the presence of two distinct clusters. The peaks of the modes correspond to the modes of the underlying distributions. However, the exact delineation between cutaneotypes may be difficult to define because the two distributions contributing to the bimodal mixture have substantial overlap in their tails.

In addition to differing in terms of taxonomic composition, the cutaneotypes are associated with psoriasis status (P <0.01). Non-lesion specimens from affected individuals have approximately even assignment to cutaneotypes, while most of the control specimens belong to cutaneotype 1 and the lesional specimens are dominated by cutaneotype 2 (Figure 4C). The difference in cutaneotype composition between the lesion and control samples yields a high odds ratio (3.52 (95% CI 1.44 to 8.98)). The evidence for existence of these two distinct cutaneotypes in the data is further corroborated by analyses based on the gap statistic [42] (Additional file 1: Figure S5 and Supplemental Methods section). In total, these findings provide strong support of the existence of distinct types of skin microbiota, cutaneotypes that differ in prevalence in psoriasis. Cutaneotypes serve as the description of the internal structure of our dataset, which may or may not be reproduced in other data. We caution the reader against overgeneralization of the reported preliminary observations to other datasets. The utility of clustering analysis rests in part on its reproducibility and generalizability; future studies of the cutaneous microbiota should examine this issue.

Psoriasis status is the major source of variability in microbial communities

We next attempted to decompose the variability of the microbiota into major components to identify possible associations with psoriasis status. We used PCoA on weighted Unifrac distances to examine the association of the entire cutaneous microbiome with specimen type (Figure 5). The first two axes are clearly separated from the rest of the components (Figure 5A), therefore we chose to represent the cutaneous microbiome of our study in these two dimensions. Only approximately 24% of total variation in the dataset is explained by the first plane, indicating that most of the variation is captured in multiple other dimensions. The high variability of the composition of individual cutaneous microbiota is indicated by the observation that 11 and 84 axes are necessary to represent, >50% and >90%, respectively, of the variability present in the specimen comparisons (data not shown). We tested the first two principal coordinates for association with psoriasis status. The microbiota from the lesions is distinguishable from that of unaffected and control skin along the first principal axis (Figure 5B) (P <0.001, ANOVA). This means that the major axis of variability as determined by Unifrac-based principal coordinates analysis is aligned with psoriasis status. Samples from control and unaffected skin are located across the origin from the lesion samples on PC1, but not on PC2 (Figure 5C). The individual samples from the different specimen types occupy common areas on two-dimensional and three-dimensional plots, without clear clustering according to psoriasis status. These findings reiterate the increased diversity in the psoriasis lesions and establish psoriasis status as one of the main sources of variability in our datasets, despite the substantial remaining variation unaccounted for by the clinical data.

We did not find evidence that links subject gender, ethnicity, age, or family history of psoriasis with the variability of the cutaneous microbiota. However, affected body site and month of sample collection each were associated with microbiota composition (P <0.001). Interestingly, in our data the cutaneous microenvironment (dry vs sebaceous) is not an important factor in accounting for variability of the skin microbiota (P >0.05). We further examined the association of the microbiota with the collection date for possible biases. We re-evaluated the model by considering only lesion and unaffected specimens (which are collected from the same individual) and stratifying the analyses by the subject. Under the modified model, the month of collection was not significant, suggesting that the association we originally observed is entirely due to the high intersubject variability, which is captured in part by the collection date variable.

Correlation analysis of psoriasis severity

An important prospect of studying the microbiota in psoriatic lesions is to identify taxa that correlate with disease severity, which may lead to better understanding of microbial roles in disease progression. We examined the identified distinct cutaneotypes for association with disease severity, as measured by the PASI, BSA and PGA scores. The cutaneotypes were not significantly associated with any of these severity measures (Additional file 1: Figure S6). However, several specific major taxa present across many taxonomic levels are associated with severity (Table 5). We did not find any taxa that were simultaneously correlated with all three measures of severity. This suggests that there is a weak link between the clinical assessment of disease severity and microbial colonization of the lesions.

Longitudinal analysis of the cutaneous microbiota

The longitudinal basis of sample collection allowed us to assess the effects of anti-inflammatory therapies on the composition of the microbiota communities. The psoriasis patients showed an overall improvement in the clinical severity of the lesions during the course of the treatment (Additional file 1: Table S6). Because we were only able to obtain follow-up samples from 17 and 9 subjects at 12-week and 36-week timepoints, we used the control specimens to reduce variability in the α diversity estimates. Our initial examination of the α diversity in the longitudinal cohort demonstrated high variability in the measurements of richness and Shannon index. Therefore, we looked for a way to minimize the variability in the data by restricting our analysis to triplets and viewing them in the light of natural variability that is represented by the control subjects. To do so, for each triplet we subtracted the α diversity of the control specimen from those of the unaffected and lesion specimens. The utility of this approach is in that in our case it allowed for consistent longitudinal trends to emerge in otherwise highly variable data. Thus, the longitudinal α diversity results are presented in terms of α diversity relative to the controls in each triplet (Figure 6). Although no statistically significant difference was observed between lesion and unaffected groups or longitudinally within groups, we observed several consistent patterns. Richness decreased over time (and treatment) in both lesion and unaffected specimens at all taxonomical levels, except for 97% identity OTUs, where the lesion demonstrated similar decrease, while the unaffected cutaneous flora rebounded to baseline levels at week 36. We also observed an increase in evenness (Shannon index), followed by a later decline. Both observations are consistent with findings in the cross-sectional cohort and lead to the similar conclusion that the abundance of a taxon or a group of taxa was increased, leading to elimination of other taxa (decreasing richness) and lower abundances of the others (decline in evenness). These patterns are preliminary and will require a larger cohort to confirm or reject.

The weighted and unweighted Unifrac intragroup β diversity estimates were similar among the three groups at all sampling times (Figure 7). Significant intragroup differences only were observed at baseline (Figure 7A) in case of unweighted Unifrac. One explanation for the absence of robust statistical differences in α and β diversity between the skin types is insufficient power in the longitudinal study.

However, the longitudinal samples allowed us to confirm the observation that cutaneotype 2 is most prevalent in psoriasis subjects even in the course of treatment. We found that although the prevalence changed over time, the relative trend of increased prevalence of cutaneotype 2 in unaffected and lesion specimens was consistent at all three timepoints (Figure 8). As can be seen from examining just the controls, subjects may switch from cutaneotype 1 to 2 within the course of the experiment. The source for this dynamic behavior is unknown and cannot be determined from the data we have, but suggests that in the skin as with other body sites (for example, the vagina; see [3]), the major clustering types may shift for an individual. This may be an important observation. This serves as an initial evidence of utility of cutaneotype clustering, which needs to be further examined in the future. A possible limitation to this observation is that the apparent increase in the proportion of cutaneotype 2 samples may in part be due to decreasing number of triplets over time in the longitudinal cohort (17 at baseline and 12-week mark and 9 at 36-week mark). A larger longitudinal cohort will be needed to see if this is a real limitation.

We also examined the longitudinal stability of the observed differential abundance of skin-associated genera (Corynebacterium, Propionibacterium, Staphylococcus, and Streptococcus). We found that the total combined relative abundance of these taxa remained surprisingly stable in the control specimens throughout the course of the experiment, while their abundance in the lesion and unaffected samples varied over time and was consistently different from that of control (Figure 9). The lesion specimens contained a higher proportion of bacteria belonging to these skin genera, and that proportion increased slightly from 40.9% at baseline to 46.7% at week 36. The abundance of these bacterial genera in the unaffected skin (25.5%) was similar to the control (27.8%) at baseline, but increased dramatically after 12 weeks of treatment (52.3%) before declining (37.1%) after further 24 weeks of treatment. This suggests that normally stable microbiota are perturbed by the treatment. A larger cohort of subjects is needed to examine the clinical significance, strength, and stability of such perturbation.